Add MPI-complete direct scalar boundary-flux integrals

docs/api/solvers.md

@@ -22,6 +22,27 @@ Viscoelastic extension of the Stokes solver.
 
 ## Scalar Equations
 
+### Integrated Boundary Flux
+
+After solving a continuous scalar problem, call collectively on all ranks:
+
+```python
+total_flux = poisson.boundary_flux_integral("Top")
+```
+
+This sums consistent nodal volume reactions on the requested essential
+boundary. It has the same raw CBF sign as `boundary_flux()`, without mass
+recovery, a temporary flux field, or boundary quadrature. Boundary membership
+is propagated through the PETSc point SF so ranks sharing a boundary node
+include their partial reactions even when they hold no labelled facet.
+
+The return value is an integral: no mean removal, area division or Nusselt
+normalization is applied. Divide by the boundary area and the appropriate
+reference conductive flux when a normalized diagnostic is required. At
+intersections of driven walls, a shared nodal reaction mixes contributions
+from both walls; this method does not separate those by facet. Use
+`boundary_flux()` or `boundary_flux_field()` for pointwise values instead.
+
 ### SNES_Poisson
 
 ```{eval-rst}

src/underworld3/cython/petsc_generic_snes_solvers.pyx

@@ -3220,6 +3220,19 @@ class SolverBaseClass(uw_object):
         return _bff(self, boundary, field, mass=mass, remove_mean=remove_mean,
                     scale=scale, normal=normal)
 
+    def boundary_flux_integral(self, boundary):
+        r"""Integrated scalar CBF flux through ``boundary``.
+
+        This is the direct integral diagnostic for quantities such as Nusselt
+        numbers. It sums the consistent scalar nodal reactions collectively,
+        avoiding pointwise de-smearing and a temporary flux MeshVariable. Use
+        :meth:`boundary_flux` or :meth:`boundary_flux_field` when nodal values
+        are required.
+        """
+        from underworld3.utilities.boundary_flux import boundary_flux_integral as _bfi
+        return _bfi(self, boundary)
+
+
 ## Specific to dimensionality
 
 

src/underworld3/utilities/boundary_flux.py

@@ -783,6 +783,57 @@ def boundary_flux(solver, boundary, mass="auto", remove_mean=False, normal=None)
     return xs, (np.column_stack(cols) if nodes else np.zeros((0, ncomp)))
 
 
+def boundary_flux_integral(solver, boundary):
+    r"""Return the integrated scalar flux through ``boundary``.
+
+    For a scalar solver, summing the consistent nodal reactions on the queried
+    boundary gives :math:`\int_\Gamma F\cdot\hat n\,d\Gamma` directly because
+    the boundary basis is a partition of unity. The raw reactions are partial
+    on partition-cut nodes, so the final sum is collective across ranks.
+
+    This path is intended for integral diagnostics such as Nusselt numbers. It
+    avoids pointwise boundary-mass recovery and a temporary MeshVariable. Use
+    :meth:`boundary_flux` or :meth:`boundary_flux_field` when nodal values are
+    required.
+
+    Call collectively after solving a continuous scalar problem with an
+    essential boundary condition on the queried boundary. The sign is the
+    raw CBF residual sign, identical to ``boundary_flux``; no mean is removed
+    and no area normalization is applied. For a mean flux divide by the
+    boundary area. At intersections of driven boundaries a nodal reaction
+    includes both contributions, so this is not a facet-separated flux there.
+    """
+    dm = solver.dm
+    if dm.getLocalSection().getFieldComponents(0) != 1:
+        raise ValueError(
+            "boundary_flux_integral requires a scalar solver field; use "
+            "boundary_flux(..., normal=...) for vector traction."
+        )
+    nodes, lsec, _csec, _cvec, _v0, _v1, _edge_nodes = _boundary_field_nodes(
+        solver, boundary, field_id=0
+    )
+    # A rank may share a boundary node without holding a labelled facet.
+    # Propagate membership through the point SF before summing raw reactions.
+    boundary_points = np.zeros(dm.getChart()[1], dtype=np.int32)
+    for point, _slot, _coordinate in nodes:
+        boundary_points[point] = 1
+    if dm.comm.size > 1:
+        sf = dm.getPointSF()
+        roots = boundary_points.copy()
+        sf.reduceBegin(MPI.INT32_T, boundary_points, roots, MPI.MAX)
+        sf.reduceEnd(MPI.INT32_T, boundary_points, roots, MPI.MAX)
+        sf.bcastBegin(MPI.INT32_T, roots, boundary_points, MPI.MAX)
+        sf.bcastEnd(MPI.INT32_T, roots, boundary_points, MPI.MAX)
+        np.maximum(boundary_points, roots, out=boundary_points)
+    ra = np.asarray(solver._assemble_volume_reaction()).ravel()
+    local_integral = sum(
+        float(np.sum(ra[lsec.getFieldOffset(point, 0):
+                        lsec.getFieldOffset(point, 0) + lsec.getFieldDof(point, 0)]))
+        for point in np.flatnonzero(boundary_points)
+    )
+    return float(dm.comm.tompi4py().allreduce(local_integral, op=MPI.SUM))
+
+
 def write_boundary_scalar_field(solver, field, value_by_key, dim):
     """Write ``value_by_key`` (coordinate-key → scalar) onto a scalar MeshVariable
     ``field`` at the matching nodes; interior nodes untouched. Returns ``field``.

tests/parallel/test_1065_boundary_flux_parallel.py

@@ -21,6 +21,8 @@
 
 # SERIAL reference: BdIntegral of the flux field over Bottom. `python <thisfile>`.
 GOLDEN_BDFLUX = -1.731543e-01
+ANALYTIC_DIRECT_INTEGRAL = -2.0 / np.sinh(np.pi)
+GOLDEN_DIRECT_INTEGRAL = -1.731790673330021e-01
 
 
 def _flux_diagnostics(res=48):
@@ -44,6 +46,7 @@ def _flux_diagnostics(res=48):
     xs, flux = poisson.boundary_flux("Bottom")
     poisson.boundary_flux_field("Bottom", q)
     bd_q = float(uw.maths.BdIntegral(mesh=mesh, fn=q.sym[0], boundary="Bottom").evaluate())
+    direct_integral = poisson.boundary_flux_integral("Bottom")
 
     # gather + dedup for a whole-boundary relL2 vs analytic (on rank 0, then bcast)
     comm = uw.mpi.comm
@@ -61,17 +64,26 @@ def _flux_diagnostics(res=48):
         c = np.dot(F, q_an) / (np.linalg.norm(F) * np.linalg.norm(q_an))
         F = F if c >= 0 else -F
         relL2 = float(np.linalg.norm(F - q_an) / np.linalg.norm(q_an))
-    return bd_q, comm.bcast(relL2, root=0)
+    return bd_q, direct_integral, comm.bcast(relL2, root=0)
 
 
 def test_boundary_flux_partition_independent():
     """boundary_flux reproduces the serial reference at np=2 and np=4 (flux boundary cut
     at np=4): both the collective BdIntegral of the flux field and the whole-boundary
     accuracy vs analytic."""
-    bd_q, relL2 = _flux_diagnostics(res=48)
+    bd_q, direct_integral, relL2 = _flux_diagnostics(res=48)
     assert np.isclose(bd_q, GOLDEN_BDFLUX, rtol=1e-5, atol=0), (
         f"BdIntegral flux differs serial vs np={uw.mpi.size}: {GOLDEN_BDFLUX} vs {bd_q}")
     assert relL2 < 0.01, f"heat flux relL2 vs analytic {relL2:.4f} too large at np={uw.mpi.size}"
+    assert np.isclose(
+        direct_integral, GOLDEN_DIRECT_INTEGRAL, rtol=1.0e-10, atol=0.0
+    ), (
+        "Direct reaction integral differs from the serial reference at "
+        f"np={uw.mpi.size}: {GOLDEN_DIRECT_INTEGRAL} vs {direct_integral}"
+    )
+    assert np.isclose(
+        direct_integral, ANALYTIC_DIRECT_INTEGRAL, rtol=1.0e-7, atol=0.0
+    )
 
 
 def _uniform_flux_3d_error(degree, mass):
@@ -133,6 +145,6 @@ def test_boundary_flux_degree3_partition_independent():
 
 
 if __name__ == "__main__":
-    _b, _r = _flux_diagnostics()
+    _b, _i, _r = _flux_diagnostics()
     if uw.mpi.rank == 0:
         print(f"DIAG_FLUX bd_q={_b:.9e} relL2={_r:.4f}")

tests/test_1019_boundary_flux.py

@@ -39,17 +39,18 @@ def _heatflux_diagnostics(res=48):
     poisson.boundary_flux_field("Bottom", q)
     # field is symbolically usable
     bd_q = float(uw.maths.BdIntegral(mesh=mesh, fn=q.sym[0], boundary="Bottom").evaluate())
+    direct_integral = poisson.boundary_flux_integral("Bottom")
 
     xc = np.asarray(xs)[:, 0] if len(xs) else np.zeros(0)
     q_an = np.pi * np.sin(np.pi * xc) / np.sinh(np.pi)       # analytic outward flux
-    return np.asarray(flux), q_an, bd_q
+    return np.asarray(flux), q_an, bd_q, direct_integral
 
 
 @pytest.mark.skipif(uw.mpi.size > 1, reason="serial diagnostic: rank-local flux norms are 0/0 on non-owning ranks")
 def test_boundary_flux_scalar_heatflux_serial():
     """Surface heat flux reproduces the analytic flux to high accuracy, and its mean is
     the (analytic) Nusselt number — NOT removed."""
-    flux, q_an, bd_q = _heatflux_diagnostics(res=48)
+    flux, q_an, bd_q, direct_integral = _heatflux_diagnostics(res=48)
     corr = np.dot(flux, q_an) / (np.linalg.norm(flux) * np.linalg.norm(q_an))
     fa = flux if corr >= 0 else -flux
     relL2 = np.linalg.norm(fa - q_an) / np.linalg.norm(q_an)
@@ -59,6 +60,11 @@ def test_boundary_flux_scalar_heatflux_serial():
     assert np.isclose(abs(fa.mean()), 2.0 / np.sinh(np.pi), rtol=0.02), (
         f"mean flux {fa.mean():.4f} != Nusselt {2.0/np.sinh(np.pi):.4f}")
     assert abs(bd_q) > 0.0                                    # field populated + usable
+    # The reaction sum is the integral itself, whereas integrating the recovered
+    # pointwise field includes its finite-resolution projection error.
+    assert np.isclose(
+        abs(direct_integral), 2.0 / np.sinh(np.pi), rtol=1.0e-7, atol=0.0
+    )
 
 
 def _uniform_flux_3d(degree, mass):
@@ -103,7 +109,7 @@ def test_boundary_flux_3d_p2_lumped_rejected():
 
 
 if __name__ == "__main__":
-    _f, _a, _b = _heatflux_diagnostics()
+    _f, _a, _b, _i = _heatflux_diagnostics()
     c = np.dot(_f, _a) / (np.linalg.norm(_f) * np.linalg.norm(_a))
     print(f"corr={abs(c):.4f} relL2={np.linalg.norm((_f if c>=0 else -_f)-_a)/np.linalg.norm(_a):.4f}")
 

tests/test_1020_boundary_flux_integral.py

@@ -0,0 +1,53 @@
+"""Direct scalar reaction integrals, with exact signed conduction references.
+
+The same small tests run in serial or under mpirun; no serial golden file is
+required. Only Top and Bottom are driven, avoiding mixed corner reactions.
+"""
+
+from types import SimpleNamespace
+
+import numpy as np
+import pytest
+import underworld3 as uw
+
+pytestmark = [pytest.mark.level_2, pytest.mark.tier_b]
+
+
+@pytest.mark.parametrize("dim", [2, 3])
+@pytest.mark.parametrize("degree", [1, 2])
+def test_signed_reaction_integral(dim, degree):
+    mesh = uw.meshing.UnstructuredSimplexBox(
+        minCoords=(0.0,) * dim, maxCoords=(1.0,) * dim,
+        cellSize=0.25, regular=True, qdegree=3,
+    )
+    temperature = uw.discretisation.MeshVariable("T", mesh, 1, degree=degree)
+    poisson = uw.systems.Poisson(mesh, u_Field=temperature)
+    poisson.constitutive_model = uw.constitutive_models.DiffusionModel
+    poisson.constitutive_model.Parameters.diffusivity = 2.5
+    poisson.add_dirichlet_bc(0.0, "Bottom")
+    poisson.add_dirichlet_bc(1.0, "Top")
+    poisson.tolerance = 1e-11
+    poisson.petsc_options["snes_type"] = "ksponly"
+    poisson.solve()
+
+    fields = tuple(mesh.vars)
+    bottom = poisson.boundary_flux_integral("Bottom")
+    top = poisson.boundary_flux_integral("Top")
+    np.testing.assert_allclose([bottom, top], [-2.5, 2.5], rtol=0, atol=1e-8)
+    assert tuple(mesh.vars) == fields
+    assert abs(bottom + top) < 1e-8
+    with pytest.raises(ValueError, match="Unknown boundary"):
+        poisson.boundary_flux_integral("NotABoundary")
+    uw.pprint(f"DIRECT_FLUX dim={dim} degree={degree} ranks={uw.mpi.size} "
+              f"bottom={bottom:.12g} top={top:.12g}")
+
+
+@pytest.mark.level_1
+def test_vector_rejected_before_reaction_assembly():
+    from underworld3.utilities.boundary_flux import boundary_flux_integral
+
+    # No assembly method: validation must reject the vector before attempting it.
+    section = SimpleNamespace(getFieldComponents=lambda field: 2)
+    solver = SimpleNamespace(dm=SimpleNamespace(getLocalSection=lambda: section))
+    with pytest.raises(ValueError, match="requires a scalar solver field"):
+        boundary_flux_integral(solver, "Top")