From 6d6ededde2c132dac5e1881183cf7ced9705d109 Mon Sep 17 00:00:00 2001 From: lmoresi Date: Sat, 5 Sep 2026 18:47:05 -0700 Subject: [PATCH] Set the PETSc constants on the DS the integrals use, so expression values reach the kernels (#695) uw.maths.Integral, BdIntegral and CellWiseIntegral compile their integrands through the same JIT as the solvers, which routes every uw.function.expression to PETSc's constants array, but none of them ever called PetscDSSetConstants: the kernels read zeros, so any integrand with a viscosity, a time or another expression in it integrated to nothing, and a fresh Integral returned the same zero from the cache. Found on the DFG cylinder drag, where the viscous traction (eta is an expression) vanished and the drag read 23 to 28% low on two meshes without moving with the SUPG weights. Each class now packs the manifest and sets the constants right after the objective; the boundary integral sets them on its sandbox DS, which has its own discrete system. Regression test test_0503 covers the three classes, a changed value without recompilation, and the constitutive-flux traction that found it. Underworld development team with AI support from Claude Code Co-Authored-By: Claude Fable 5.1 Claude-Session: https://claude.ai/code/session_018T2VHUGaZiQVJ95qQ4DiSL (cherry picked from commit 5131039e27cd87218a2af64bd62656b04565754d) --- src/underworld3/cython/petsc_maths.pyx | 40 +++++++++++- ...test_0503_integral_expression_constants.py | 64 +++++++++++++++++++ 2 files changed, 101 insertions(+), 3 deletions(-) create mode 100644 tests/test_0503_integral_expression_constants.py diff --git a/src/underworld3/cython/petsc_maths.pyx b/src/underworld3/cython/petsc_maths.pyx index 5be87ac74..a4445d780 100644 --- a/src/underworld3/cython/petsc_maths.pyx +++ b/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), &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 diff --git a/tests/test_0503_integral_expression_constants.py b/tests/test_0503_integral_expression_constants.py new file mode 100644 index 000000000..64ef38ed6 --- /dev/null +++ b/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)