From 88807c26db09688efdbcf185ffe5efde343f0566 Mon Sep 17 00:00:00 2001 From: lmoresi Date: Thu, 7 May 2026 10:31:01 +1000 Subject: [PATCH 1/2] Revert "Merge pull request #172 from underworldcode/bugfix/integral-evaluate-linear-growth-171" This reverts commit fe05c632d84c78930f76690ec18579a836f79ae6, reversing changes made to 259958bca6907d37c20b2680aa935af2e4c8beb4. --- src/underworld3/cython/petsc_maths.pyx | 35 ++++++++------------------ 1 file changed, 10 insertions(+), 25 deletions(-) diff --git a/src/underworld3/cython/petsc_maths.pyx b/src/underworld3/cython/petsc_maths.pyx index 465bc215..29e26a2a 100644 --- a/src/underworld3/cython/petsc_maths.pyx +++ b/src/underworld3/cython/petsc_maths.pyx @@ -102,46 +102,31 @@ class Integral: raise RuntimeError("Integral evaluation for Dyadic integrands not supported.") - # BUGFIX(#171): clone the mesh DM and create a fresh DS per call, - # rather than mutating the shared mesh DS. The previous code path - # was: - # ds = self.mesh.dm.getDS() # shared DS - # PetscDSSetObjective(ds, 0, fn) # mutate it - # DMPlexComputeIntegralFEM(...) - # PetscDSSetObjective on the shared DS makes - # DMPlexComputeIntegralFEM cost grow linearly with repeated - # evaluate() calls in long simulations — observed in convection - # diagnostics where t_vol grew from 5 s to 833 s over ~3000 calls. - # The sister class CellWiseIntegral already takes the - # clone-DM + createDS path; this brings Integral in line. - # Explicit destroy of the clone at the end keeps long runs bounded. - mesh = self.mesh + self.dm = self.mesh.dm # .clone() + mesh=self.mesh _getext_result = getext(self.mesh, JITCallbackSet(residual=(self.fn,)), self.mesh.vars.values(), verbose=verbose) cdef PtrContainer ext = _getext_result.ptrobj # Pull out vec for variables, and go ahead with the integral + self.mesh.update_lvec() - a_global = self.mesh.dm.getGlobalVec() - self.mesh.dm.localToGlobal(self.mesh.lvec, a_global) + a_global = self.dm.getGlobalVec() + self.dm.localToGlobal(self.mesh.lvec, a_global) cdef Vec cgvec cgvec = a_global - 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() + cdef DM dm = self.dm + cdef DS ds = self.dm.getDS() cdef PetscScalar val_array[256] + # Now set callback... ierr = PetscDSSetObjective(ds.ds, 0, ext.fns_residual[0]); CHKERRQ(ierr) - ierr = DMPlexComputeIntegralFEM(dmc.dm, cgvec.vec, &(val_array[0]), NULL); CHKERRQ(ierr) + ierr = DMPlexComputeIntegralFEM(dm.dm, cgvec.vec, &(val_array[0]), NULL); CHKERRQ(ierr) - self.mesh.dm.restoreGlobalVec(a_global) - dmc.destroy() + self.dm.restoreGlobalVec(a_global) # We're making an assumption here that PetscScalar is same as double. # Need to check where this may not be the case. From d65fe9a22668c12d52b17642ac544508b480c3a4 Mon Sep 17 00:00:00 2001 From: lmoresi Date: Thu, 7 May 2026 10:54:06 +1000 Subject: [PATCH 2/2] CellWiseIntegral: fix self.dim typo + xfail regression tests for the field-cloning bug MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit While reverting PR #172 (which broke Integral.evaluate), audited the sister CellWiseIntegral class — the PR's premise was that CellWiseIntegral "already takes the right path". It does not: 1. CellWiseIntegral.evaluate() referenced self.dim (no such attribute), so any call to the class crashed with AttributeError before exhibiting the deeper bug. This commit fixes it to self.mesh.dim. 2. With the typo fixed, CellWiseIntegral exhibits the same field-cloning bug that broke Integral in PR #172: cloning mesh.dm and attaching a fresh single-field discretisation, but then integrating against mesh.dm.getGlobalVec() (packed for the multi-field layout), produces a ~2x over-count on the unit square. CellWiseIntegral has zero pre-existing test coverage. This commit adds two regression tests marked xfail with a clear reason; they will start passing once CellWiseIntegral is rewritten to integrate against mesh.dm + getDS() directly (the pre-PR-172 Integral pattern). An inline TODO(BUG) comment in the source points at the same write-up. Underworld development team with AI support from Claude Code --- src/underworld3/cython/petsc_maths.pyx | 11 ++++++-- tests/test_0501_integrals.py | 39 ++++++++++++++++++++++++++ 2 files changed, 47 insertions(+), 3 deletions(-) diff --git a/src/underworld3/cython/petsc_maths.pyx b/src/underworld3/cython/petsc_maths.pyx index 29e26a2a..4061bf29 100644 --- a/src/underworld3/cython/petsc_maths.pyx +++ b/src/underworld3/cython/petsc_maths.pyx @@ -300,10 +300,15 @@ class CellWiseIntegral: cdef Vec cgvec cgvec = a_global - ## Does this need to be consistent with everything else ? - + # 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.dim, 1, False, -1) + cdef FE fec = FE().createDefault(self.mesh.dim, 1, False, -1) dmc.setField(0, fec) dmc.createDS() diff --git a/tests/test_0501_integrals.py b/tests/test_0501_integrals.py index 512bc6a0..309e22ca 100644 --- a/tests/test_0501_integrals.py +++ b/tests/test_0501_integrals.py @@ -166,3 +166,42 @@ def test_integrate_swarmvar_deriv_00(): assert abs(value + 2) < 0.0001 return + + +# 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") +def test_cellwise_integrate_constants(): + + calculator = uw.maths.CellWiseIntegral(mesh, fn=1.0) + cell_values = calculator.evaluate() + + # Sum of per-cell volumes equals the total mesh volume (unit square = 1.0). + assert abs(np.sum(cell_values) - 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]) + + calculator = uw.maths.CellWiseIntegral(mesh, fn=s_soln.sym[0]) + cell_values = calculator.evaluate() + + # 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 + + return