Merge pull request #174 from underworldcode/bugfix/integrals-revert

Revert PR #172 (integrals regression) + xfail CellWiseIntegral parity tests

Author: lmoresi

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.
@@ -315,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()
 

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