Volume integral evaluation time grows linearly across long runs (~80x in 3000 calls)

[lmoresi]

Symptom

In a long Rayleigh–Bénard run at Ra=1e7 (1/64 mesh, ~3100 timesteps), per-step solver costs (Stokes + advection–diffusion) stayed flat while the volume-integral diagnostic time (t_vol — used for Vrms, mean_T, Vmax) grew linearly with the number of integral evaluations.

Quartile of run	mean t_vol	first sample	last sample
Q1 (steps 645–1425, 157 diag. calls)	44 s	5 s	84 s
Q2 (1430–2210)	149 s	84 s	418 s
Q3 (2215–2990)	267 s	194 s	833 s
Q4 (2995–3770)	348 s	296 s	576 s

Linear fit over the diagnostic samples: +0.13 s per simulation step (≈ +0.65 s per integral call). By step ~3700 the diagnostic step costs ~410 s — roughly 20× the actual solve (t_stokes + t_adv ≈ 21 s, also flat).

For comparison, t_stokes and t_adv over the same 3134 steps had slopes of +0.0001 s/step and +0.00009 s/step respectively (≈ 3% drift, well inside noise). So the growth is specific to the volume-integral path.

(Plot: t_vol vs step alongside t_stokes and t_adv, with the linear fit, can be attached.)

Reproduction

• Script: standard convection workflow (docs/examples/workflows/convection.py) running _compute_diagnostics(...) — calls three Integral.evaluate() invocations per diagnostic step (Vrms, mean_T, Vmax).
• The Integral objects themselves are cached (one object per integrand, reused — see convection_config.py:378–407), so this is not JIT recompilation: the JIT cache hits every call.

Suspect codepath

Integral.evaluate() in src/underworld3/cython/petsc_maths.pyx:

self.dm = self.mesh.dm  # .clone()                           # line 89
...
cdef DS ds = self.dm.getDS()                                  # line 105
ierr = PetscDSSetObjective(ds.ds, 0, ext.fns_residual[0])     # line 109
ierr = DMPlexComputeIntegralFEM(dm.dm, cgvec.vec, ...)        # line 110

Every evaluate() call:

1. takes a reference to the shared mesh DS (no clone, no clear),
2. registers a new objective callback via PetscDSSetObjective(ds.ds, 0, fn),
3. computes the integral.

The DS is never cleared between calls. Hypothesis (the obvious one): the integration function pointers are being accumulated on the shared DS rather than replaced — every diagnostic call adds another callback that gets walked at integration time, producing the linear growth observed.

Supporting evidence inside the same file

The sister class CellWiseIntegral.evaluate() (lines 245–302) takes the opposite approach for the same operation:

dmc = self.mesh.dm.clone()    # line 287
dmc.setField(0, fec)
dmc.createDS()                # line 290
...
PetscDSSetObjective(...)      # on the fresh DS

It clones the DM and creates a fresh DS for every integration. The asymmetry suggests that the accumulation issue was at least implicitly known when CellWiseIntegral was written — Integral.evaluate() inherits the bug.

The commented-out # .clone() on line 89 also looks like a hint that someone tried cloning here and reverted.

Proposed fix (to investigate, not yet implemented)

Either:

1. Clear the DS objective before re-setting: before line 109, call something equivalent to PetscDSSetObjective(ds.ds, 0, NULL) — or whatever PETSc exposes for resetting per-field objectives.
2. Match CellWiseIntegral: clone the DM and createDS() per call (heavier, but guaranteed isolation).

Option (1) is preferred if PETSc supports it cleanly, since cloning the DM allocates per call.

Impact

• The convection example becomes diagnostic-bound after a few thousand steps. For longer scientific runs (the 3134-step run above ran for ~24 wall-hours; ~85% of the wall-time was in t_vol after step 2200), this is a hard ceiling on what the workflow can do unattended.
• Affects any user calling mesh.integrate(...) (or Integral.evaluate()) repeatedly across a long run — diagnostics, time-series, mass-conservation checks, convergence monitors.
• Likely related to the broader DM lifecycle / object-accumulation pattern noted elsewhere in the codebase.

Discovered while

Validating the SLCN advection–diffusion fix and the Rayleigh–Bénard workflow tooling — see the long Ra=1e7 1/64 cubic baseline run in /Users/lmoresi/+Simulations/RayleighBenard/Ra1e7_1over64/.

Underworld development team with AI support from Claude Code

[lmoresi]

Closing — the linear-growth symptom does not reproduce on current `development`.

Investigation outcome

Ran the full reproducer documented at `~/.claude/plans/repro-issue-171-integral-growth.md` (mirrors the original 1/64 Ra=1e7 setup: `UnstructuredSimplexBox`, P2-P1 Stokes, P3 T, qdegree=3, three cached `Integral` objects), at cellsize 1/32 with 100 calls per integrand.

	t_first	t_steady	ratio	linear slope	verdict
unpatched dev	27.7 ms	16.2 ms	0.59x	flat (-0.08%/call)	NOT reproduced
with #172 (clone-DM fix)	12.0 ms	5.1 ms	0.43x	flat (-0.20%/call)	NOT reproduced

The original predicted trajectory at 1/32 was 1.3 s → 17 s (13× ratio); actual unpatched today is 28 ms → 16 ms — ~1000× lower than predicted-with-bug, and entirely flat. Some intermediate landing between the original report and now killed the pathology — most likely candidates are PR #155 (swarm-routed point eval), PR #160 (SNES + DM-hierarchy lifecycle leak), or PR #161 (ETD lifecycle changes), but bisecting which one is out of scope.

What landed regardless

PR #172 (merged) replaced `Integral.evaluate()`'s shared-mesh-DS path with a per-call clone-DM + `createDS` (matching the working `CellWiseIntegral` pattern). On current dev this gives a clean ~3× per-call speedup for volume integrals (16 ms → 5 ms), plus structural improvements:

• Eliminates the only remaining shared-DS-mutation site in `Integral.evaluate()` — closes the door on this pathology recurring as other lifecycle code shifts.
• Explicit `dmc.destroy()` so long runs don't accumulate cloned DMs.
• Closes the asymmetry with `CellWiseIntegral` that originally pointed at this site.

Reproducer (embedded for future regressions)

If volume-integral cost ever starts climbing during a long run again, save this as `/tmp/repro_issue_171.py` and run with `pixi run -e default python -u /tmp/repro_issue_171.py` (defaults: cellsize 1/32, 100 calls; takes ~30 s). The script ends with a machine-readable `# REPRODUCED` / `# NOT reproduced` verdict and a CSV per-call timing trail. Suspect line for any return: `src/underworld3/cython/petsc_maths.pyx` — was line 109 (`PetscDSSetObjective` on the shared mesh DS) before #172.

```python
"""
Reproducer for underworld3 issue #171.

uw.maths.Integral.evaluate() per-call wall time grew linearly across thousands
of calls in a long convection run. This script exercises the integral path in
isolation: build a faithful mesh + variables, populate v and T with non-trivial
fields, then call evaluate() on three cached Integral objects in a tight loop
with per-call timing.

Run:
pixi run -e default python -u /tmp/repro_issue_171.py
pixi run -e default python -u /tmp/repro_issue_171.py --cellsize 0.015625 --n-calls 50

Args:
--cellsize Mesh cell size. 1/32 (default) is fast; 1/64 matches the original.
--qdegree PETSc quadrature degree (default 3, matches the original).
--T-degree Temperature polynomial degree (default 3, matches the original).
--n-calls Number of evaluate() iterations per integrand (default 100).
--solve-once Also assemble + solve Stokes once before the loop, to seed
any DS state the solvers might leave behind. Off by default —
Integral DS is independent of solver DS, so solving should
not be needed for the bug to appear.

Output (stdout, CSV-friendly):
n,t_vol,t_meanT,t_v2,t_total # per call
# linear fit + reproduced/not-reproduced verdict at the end.
"""
import argparse
import time

import numpy as np
import sympy

import underworld3 as uw

def main():
ap = argparse.ArgumentParser()
ap.add_argument("--cellsize", type=float, default=1.0 / 32)
ap.add_argument("--qdegree", type=int, default=3)
ap.add_argument("--T-degree", type=int, default=3)
ap.add_argument("--n-calls", type=int, default=100)
ap.add_argument("--solve-once", action="store_true")
args = ap.parse_args()

print(
    f\"# repro_issue_171: cellsize={args.cellsize:g} (1/{1/args.cellsize:.0f}), \"
    f\"qdegree={args.qdegree}, T_degree={args.T_degree}, \"
    f\"n_calls={args.n_calls}, solve_once={args.solve_once}\",
    flush=True,
)

mesh = uw.meshing.UnstructuredSimplexBox(
    minCoords=(0.0, 0.0),
    maxCoords=(1.0, 1.0),
    cellSize=args.cellsize,
    regular=False,
    qdegree=args.qdegree,
)
v = uw.discretisation.MeshVariable(\"U\", mesh, mesh.dim, degree=2)
p = uw.discretisation.MeshVariable(\"P\", mesh, 1, degree=1)
T = uw.discretisation.MeshVariable(\"T\", mesh, 1, degree=args.T_degree)

print(
    f\"# DOFs: v={v.coords.shape[0]}, p={p.coords.shape[0]}, T={T.coords.shape[0]}\",
    flush=True,
)

sx, sy = sympy.symbols(\"sx sy\")
T_expr = 1.0 - sy + 0.1 * sympy.sin(sympy.pi * sx) * sympy.exp(
    -((sy - 0.5) ** 2) / 0.05
)
vx_expr = sympy.sin(sympy.pi * sx) * sympy.cos(sympy.pi * sy)
vy_expr = -sympy.cos(sympy.pi * sx) * sympy.sin(sympy.pi * sy)
T_fn = sympy.lambdify((sx, sy), T_expr, \"numpy\")
vx_fn = sympy.lambdify((sx, sy), vx_expr, \"numpy\")
vy_fn = sympy.lambdify((sx, sy), vy_expr, \"numpy\")

Tc = T.coords
T.array[:, 0, 0] = T_fn(Tc[:, 0], Tc[:, 1])
Vc = v.coords
v.array[:, 0, 0] = vx_fn(Vc[:, 0], Vc[:, 1])
v.array[:, 0, 1] = vy_fn(Vc[:, 0], Vc[:, 1])

if args.solve_once:
    stokes = uw.systems.Stokes(mesh, velocityField=v, pressureField=p)
    stokes.constitutive_model = uw.constitutive_models.ViscousFlowModel
    stokes.constitutive_model.Parameters.shear_viscosity_0 = 1.0
    stokes.bodyforce = sympy.Matrix([0.0, 1e7 * T.sym[0]])
    for face in (\"Top\", \"Bottom\", \"Left\", \"Right\"):
        stokes.add_dirichlet_bc((0.0, 0.0), face)
    print(\"# solving Stokes once...\", flush=True)
    t0 = time.time()
    stokes.solve()
    print(f\"# Stokes solve: {time.time() - t0:.2f}s\", flush=True)

integrator_vol = uw.maths.Integral(mesh, 1)
integrator_meanT = uw.maths.Integral(mesh, T.sym[0])
integrator_v2 = uw.maths.Integral(mesh, v.sym.dot(v.sym))

print(\"# warmup (JIT)...\", flush=True)
t0 = time.time()
_ = integrator_vol.evaluate()
_ = integrator_meanT.evaluate()
_ = integrator_v2.evaluate()
print(f\"# warmup total: {time.time() - t0:.3f}s\", flush=True)

print(f\"# loop: {args.n_calls} calls per integrand\", flush=True)
print(\"n,t_vol,t_meanT,t_v2,t_total\", flush=True)
rows = []
for n in range(1, args.n_calls + 1):
    t0 = time.time(); _ = integrator_vol.evaluate();   t_vol = time.time() - t0
    t0 = time.time(); _ = integrator_meanT.evaluate(); t_mT  = time.time() - t0
    t0 = time.time(); _ = integrator_v2.evaluate();    t_v2  = time.time() - t0
    t_tot = t_vol + t_mT + t_v2
    rows.append((n, t_vol, t_mT, t_v2, t_tot))
    print(f\"{n},{t_vol:.4f},{t_mT:.4f},{t_v2:.4f},{t_tot:.4f}\", flush=True)

a = np.array(rows)
ns, tt = a[:, 0], a[:, 4]
slope, icpt = np.polyfit(ns, tt, 1)
ratio = tt[-1] / tt[0] if tt[0] > 0 else float(\"inf\")
print()
print(f\"# linear fit:  t_total = {icpt:.4f}s + ({slope:+.5f} s/call) * n\")
print(f\"# t_1   = {tt[0]:.4f}s\")
print(f\"# t_n   = {tt[-1]:.4f}s   (n={int(ns[-1])})\")
print(f\"# ratio = {ratio:.2f}x\")
pct = (slope / icpt * 100) if icpt > 0 else 0.0
if pct > 1.0:
    print(f\"# REPRODUCED: slope = {pct:.1f}% of baseline per call.\")
else:
    print(f\"# NOT reproduced: slope = {pct:.2f}% of baseline per call (within noise).\")

if name == "main":
main()
```

Predicted trajectory (for future regressions)

The original 1/64 long-run data fit `t_n ≈ 5 s + 0.65 × n_call`. Linearly scaling by mesh DOF count gives:

cellsize	predicted t₁ if bug returns	predicted t₁₀₀ if bug returns	actual today (no bug)
1/16	~0.3 s	~4.3 s	<0.05 s
1/32	~1.3 s	~17 s	~0.02 s
1/64	~5 s	~70 s	(untested today)

Bug present → final-vs-first ratio ≈ 13× at 1/32; clear monotone climb visible by call ~10. Bug fixed/absent → roughly flat, ratio ≈ 1×, slope under ~1% of baseline per call.

Underworld development team with AI support from Claude Code