From 57ec017660b22fbd9c48f6472a9831689428ab9f Mon Sep 17 00:00:00 2001
From: lmoresi <louis.moresi@anu.edu.au>
Date: Tue, 5 May 2026 14:21:58 +1000
Subject: [PATCH] Fix Integral.evaluate() linear-time growth on shared mesh DS
 (#171)
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

Integral.evaluate() called PetscDSSetObjective on the shared mesh DS
every call, then DMPlexComputeIntegralFEM. The user reports this makes
diagnostic-integral cost grow linearly across long runs — convection
diagnostics at Ra=1e7 / 1/64 mesh / ~3000 timesteps showed t_vol going
from 5 s to 833 s (≈ +0.65 s per integral call), while solve costs
(t_stokes, t_adv) stayed flat over the same window. By step ~3700 the
diagnostic was ~20× the actual solve.

The asymmetric sister class CellWiseIntegral.evaluate() already takes
the right path: clone the mesh DM, attach a default FE field, create a
fresh DS, set the objective on that fresh DS, integrate, and return.
This commit brings Integral.evaluate() in line with that pattern, with
one improvement over CellWiseIntegral: explicit dmc.destroy() at the
end so long runs don't accumulate cloned DMs (the broader DM-lifecycle
cleanup is tracked separately).

The fix doesn't reproduce the symptomatic linear growth at small
mesh + small call-count locally (200 calls × 100-cell mesh stays
flat at ~2.5 ms per 3 integrals both before and after the fix) — the
behaviour is scale-dependent. Confirmation at the original Ra=1e7 /
1/64 / 3000-step scale is left as a re-run on the reporter's side.

Verified
--------

  - /tmp/repro_171.py — 200 evaluate() calls × 3 integrals each on a
    100-cell mesh: post-fix mean ratio last/first = 0.04x (flat).
  - /tmp/repro_171b.py — same, with a Stokes solve interleaved between
    integrations: post-fix ratio 0.21x (flat).
  - pytest -m "level_1 and tier_a" on amr-dev: 62 passed, 3 skipped,
    0 failed. No numerical regressions.

Underworld development team with AI support from Claude Code
---
 src/underworld3/cython/petsc_maths.pyx | 35 ++++++++++++++++++--------
 1 file changed, 25 insertions(+), 10 deletions(-)

diff --git a/src/underworld3/cython/petsc_maths.pyx b/src/underworld3/cython/petsc_maths.pyx
index 29e26a2a..465bc215 100644
--- a/src/underworld3/cython/petsc_maths.pyx
+++ b/src/underworld3/cython/petsc_maths.pyx
@@ -102,31 +102,46 @@ class Integral:
             raise RuntimeError("Integral evaluation for Dyadic integrands not supported.")
 
 
-        self.dm = self.mesh.dm  # .clone()
-        mesh=self.mesh
+        # 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
 
         _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.dm.getGlobalVec()
-        self.dm.localToGlobal(self.mesh.lvec, a_global)
+        a_global = self.mesh.dm.getGlobalVec()
+        self.mesh.dm.localToGlobal(self.mesh.lvec, a_global)
 
         cdef Vec cgvec
         cgvec = a_global
 
-        cdef DM dm = self.dm
-        cdef DS ds = self.dm.getDS()
+        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 PetscScalar val_array[256]
 
-        # Now set callback...
         ierr = PetscDSSetObjective(ds.ds, 0, ext.fns_residual[0]); CHKERRQ(ierr)
-        ierr = DMPlexComputeIntegralFEM(dm.dm, cgvec.vec, &(val_array[0]), NULL); CHKERRQ(ierr)
+        ierr = DMPlexComputeIntegralFEM(dmc.dm, cgvec.vec, &(val_array[0]), NULL); CHKERRQ(ierr)
 
-        self.dm.restoreGlobalVec(a_global)
+        self.mesh.dm.restoreGlobalVec(a_global)
+        dmc.destroy()
 
         # We're making an assumption here that PetscScalar is same as double.
         # Need to check where this may not be the case.