uw.maths.Integral, BdIntegral and CellWiseIntegral read every expression constant as zero: the kernels are compiled with PETSc constants that are never set

[lmoresi]

What happens

Any uw.function.expression (UWexpression) in the integrand of uw.maths.Integral, uw.maths.BdIntegral or uw.maths.CellWiseIntegral contributes nothing: the integral of c * f is zero for any expression c, whatever its value, and a fresh Integral built after the value changed returns the same zero.

c = uw.function.expression(r"c", 2.0, "probe")
uw.maths.Integral(mesh, c * T.sym[0].diff(y)).evaluate()        # 0.0, expected 2.0
uw.maths.Integral(mesh, 2.0 * T.sym[0].diff(y)).evaluate()      # 2.0
uw.maths.BdIntegral(mesh, stokes.constitutive_model.flux[0, 1], "Top").evaluate()   # 0.0: eta is an expression

The constitutive viscosity is an expression, so every boundary traction formed from solver.stress or constitutive_model.flux integrates to its pressure part only. We found it on the DFG cylinder benchmark, where the drag came out 23 to 28% low at two mesh sizes and did not move with the SUPG weights: the whole deficit was the viscous drag, which the traction integral had silently dropped. The reaction (consistent-flux) form of the force, which uses a volume Integral of the same stress, was wrong in the same way.

Cause

The JIT (getext) routes every expression constant to PETSc's constants array (_extract_constants / PetscDSSetConstants), so that a changed value does not recompile. The solvers call PetscDSSetConstants before every solve (_update_constants). The three integral classes compile through the same getext but never set the constants on the DS they integrate with, so the kernel reads zeros. BdIntegral integrates on a sandbox DM with its own DS (DMCreateDS in UW_DMCreateBdIntegralSandbox), so the values have to go to that DS.

uw.function.evaluate is unaffected (it goes through a different path).

Fix

petsc_maths.pyx: pack the manifest (_pack_constants) and call PetscDSSetConstants on the DS each class uses, right after PetscDSSetObjective (volume, cell-wise) and on the sandbox DS (boundary). Regression test tests/test_0503_integral_expression_constants.py (volume, boundary, cell-wise, value update, and the constitutive-flux traction that found it). On branch feature/navier-stokes-supg (PR #688); it is a bug fix that belongs on development and should be extracted there.

Consequences worth checking

Any published number that came from an Integral or BdIntegral with an expression in the integrand: boundary heat flux with a diffusivity expression, tractions and dynamic topography formed from a constitutive flux, dissipation integrals with a viscosity expression. Integrands built from plain floats and mesh variables were never affected.

🤖 Generated with Claude Code

https://claude.ai/code/session_018T2VHUGaZiQVJ95qQ4DiSL