Set the PETSc constants on the DS the integrals use, so expression values reach the kernels (#695)

src/underworld3/cython/petsc_maths.pyx

@@ -1,5 +1,6 @@
 from typing import Union
 import sympy
+import numpy as np
 
 import underworld3
 import underworld3.timing as timing
@@ -15,6 +16,31 @@ cdef extern from "petsc.h" nogil:
     PetscErrorCode DMPlexComputeCellwiseIntegralFEM( PetscDM, PetscVec, PetscVec, void* )
 
 
+def _pack_manifest(manifest):
+    """The current values of the JIT constants manifest as a contiguous array."""
+    from underworld3.utilities._jitextension import _pack_constants
+    if not manifest:
+        return None
+    return np.ascontiguousarray(_pack_constants(manifest), dtype=np.float64)
+
+
+cdef _set_ds_constants(PetscDS ds, manifest):
+    """Hand the current UWexpression values to the DS the integral kernel reads.
+
+    The JIT routes every ``uw.function.expression`` in the integrand to PETSc's
+    constants array (the same mechanism the solvers use, so a changed value does
+    not recompile). A DS that never receives the values hands the kernel zeros:
+    a viscosity, a time or any other expression in an integrand silently
+    integrated to nothing (found on the cylinder drag, 2026-09-05).
+    """
+    cdef double[::1] vals
+    values = _pack_manifest(manifest)
+    if values is None or len(values) == 0:
+        return
+    vals = values
+    CHKERRQ(PetscDSSetConstants(ds, len(values), <const PetscScalar*>&vals[0]))
+
+
 def dm_force_coordinate_field(dm):
     """Force coordinate field creation and strip boundary labels from the
     coordinate DM.  Must be called after createCoordinateSpace and after
@@ -122,8 +148,9 @@ class Integral:
         cdef DS ds = self.dm.getDS()
         cdef PetscScalar val_array[256]
 
-        # Now set callback...
+        # Now set callback (and the current constant values the kernel reads)...
         ierr = PetscDSSetObjective(ds.ds, 0, ext.fns_residual[0]); CHKERRQ(ierr)
+        _set_ds_constants(ds.ds, _getext_result.constants_manifest)
         ierr = DMPlexComputeIntegralFEM(dm.dm, cgvec.vec, &(val_array[0]), NULL); CHKERRQ(ierr)
 
         self.dm.restoreGlobalVec(a_global)
@@ -290,8 +317,9 @@ class CellWiseIntegral:
         elif isinstance(self.fn, sympy.vector.Dyadic):
             raise RuntimeError("Integral evaluation for Dyadic integrands not supported.")
 
-        cdef PtrContainer ext = getext(self.mesh, JITCallbackSet(residual=(self.fn,)),
-                                       self.mesh.vars.values()).ptrobj
+        _getext_result = getext(self.mesh, JITCallbackSet(residual=(self.fn,)),
+                                self.mesh.vars.values())
+        cdef PtrContainer ext = _getext_result.ptrobj
 
         # Pull out vec for variables, and go ahead with the integral
         self.mesh.update_lvec()
@@ -316,6 +344,7 @@ class CellWiseIntegral:
         cdef DM dm = self.mesh.dm
         cdef DS ds = self.mesh.dm.getDS()
         CHKERRQ( PetscDSSetObjective(ds.ds, 0, ext.fns_residual[0]) )
+        _set_ds_constants(ds.ds, _getext_result.constants_manifest)
 
         # DMPlexComputeCellwiseIntegralFEM writes Nf scalars per cell into a
         # flat [cell*Nf + field] layout when the output vector carries no
@@ -461,6 +490,11 @@ class BdIntegral:
         cdef PetscDMLabel sandbox_label = NULL
         CHKERRQ(DMGetLabel(sandbox_dm, boundary_bytes, &sandbox_label))
 
+        # The sandbox has its own DS (DMCreateDS): the constants go there.
+        cdef PetscDS sandbox_ds = NULL
+        CHKERRQ(DMGetDS(sandbox_dm, &sandbox_ds))
+        _set_ds_constants(sandbox_ds, _getext_result.constants_manifest)
+
         # Output value
         cdef PetscScalar result = 0.0
 

tests/test_0503_integral_expression_constants.py

@@ -0,0 +1,64 @@
+"""A `uw.function.expression` inside an integrand must reach the integral kernel.
+
+The JIT routes every expression constant to PETSc's constants array (so that a
+changed value does not recompile). The integral classes compiled through that
+path but never set the values on the DS they integrate with, so the kernel read
+zeros: any integrand carrying a viscosity, a time, or any other expression
+integrated to nothing, and a fresh Integral returned the same zero from the
+cache. Found on the cylinder drag (the viscous traction vanished), 2026-09-05.
+"""
+import numpy as np
+import pytest
+import sympy
+import underworld3 as uw
+
+pytestmark = [pytest.mark.level_1, pytest.mark.tier_a]
+
+
+@pytest.fixture(scope="module")
+def setup():
+    mesh = uw.meshing.UnstructuredSimplexBox(minCoords=(0.0, 0.0), maxCoords=(1.0, 1.0),
+                                             cellSize=0.25, regular=True, qdegree=3)
+    x, y = mesh.X
+    T = uw.discretisation.MeshVariable("T_c", mesh, 1, degree=2)
+    T.array[:, 0, 0] = uw.function.evaluate(y ** 2, T.coords).reshape(-1)     # dT/dy = 2y
+    return mesh, x, y, T
+
+
+def test_volume_integral_carries_the_expression_value(setup):
+    mesh, x, y, T = setup
+    c = uw.function.expression(r"c_{v}", 2.0, "probe constant")
+    integral = uw.maths.Integral(mesh, c * T.sym[0].diff(y))               # 2 * int 2y = 2
+    assert np.isclose(integral.evaluate(), 2.0, rtol=1e-8)
+    c.sym = 3.0                                                            # a changed value, no recompile
+    assert np.isclose(integral.evaluate(), 3.0, rtol=1e-8)
+    assert np.isclose(uw.maths.Integral(mesh, c * T.sym[0].diff(y)).evaluate(), 3.0, rtol=1e-8)
+
+
+def test_boundary_integral_carries_the_expression_value(setup):
+    mesh, x, y, T = setup
+    c = uw.function.expression(r"c_{b}", 2.0, "probe constant")
+    integral = uw.maths.BdIntegral(mesh, c * T.sym[0].diff(y), "Top")       # 2 * 2 * length 1
+    assert np.isclose(integral.evaluate(), 4.0, rtol=1e-8)
+    c.sym = 0.5
+    assert np.isclose(integral.evaluate(), 1.0, rtol=1e-8)
+
+
+def test_cellwise_integral_carries_the_expression_value(setup):
+    mesh, x, y, T = setup
+    c = uw.function.expression(r"c_{c}", 2.0, "probe constant")
+    cells = uw.maths.CellWiseIntegral(mesh, c * T.sym[0].diff(y)).evaluate()
+    assert np.isclose(np.asarray(cells).sum(), 2.0, rtol=1e-8)
+
+
+def test_constitutive_flux_in_a_boundary_integral(setup):
+    """The case that found it: the viscous traction on a wall."""
+    mesh, x, y, T = setup
+    v = uw.discretisation.MeshVariable("U_c", mesh, 2, degree=2)
+    p = uw.discretisation.MeshVariable("P_c", mesh, 1, degree=1)
+    stokes = uw.systems.Stokes(mesh, v, p)
+    stokes.constitutive_model = uw.constitutive_models.ViscousFlowModel
+    stokes.constitutive_model.Parameters.shear_viscosity_0 = 2.0
+    v.array[:, 0, :] = uw.function.evaluate(sympy.Matrix([[y ** 2, 0.0]]), v.coords).reshape(-1, 2)
+    sigma_xy = stokes.constitutive_model.flux[0, 1]                        # 2 eta (du/dy)/2 = 2y * 2 / ... = eta * 2y
+    assert np.isclose(uw.maths.BdIntegral(mesh, sigma_xy, "Top").evaluate(), 4.0, rtol=1e-8)