From f6b780a5f1431a61be70c4506314f758a4f6861f Mon Sep 17 00:00:00 2001 From: Tyagi Date: Wed, 9 Sep 2026 11:32:54 +1000 Subject: [PATCH] Add MPI-complete scalar boundary-flux integrals for diagnostic reductions Expose boundary_flux_integral() as a collective raw scalar CBF reaction sum. Propagate boundary membership across the PETSc point SF before summing, so partial nodal reactions are retained on ranks without a labelled local boundary facet. Avoid pointwise recovery and temporary flux fields. Document sign, normalization, essential-boundary and corner constraints. Include exact signed P1/P2 triangle/tetrahedron tests and existing smooth-flux references. The clean upstream-based build passes 20 serial tests and nine tests on eight MPI ranks. No unrelated mantle solver changes are included. --- docs/api/solvers.md | 21 ++++++++ .../cython/petsc_generic_snes_solvers.pyx | 13 +++++ src/underworld3/utilities/boundary_flux.py | 51 ++++++++++++++++++ .../test_1065_boundary_flux_parallel.py | 18 +++++-- tests/test_1019_boundary_flux.py | 12 +++-- tests/test_1020_boundary_flux_integral.py | 53 +++++++++++++++++++ 6 files changed, 162 insertions(+), 6 deletions(-) create mode 100644 tests/test_1020_boundary_flux_integral.py diff --git a/docs/api/solvers.md b/docs/api/solvers.md index 55027e56e..ac7f51a53 100644 --- a/docs/api/solvers.md +++ b/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} diff --git a/src/underworld3/cython/petsc_generic_snes_solvers.pyx b/src/underworld3/cython/petsc_generic_snes_solvers.pyx index aa22e699e..238efa507 100644 --- a/src/underworld3/cython/petsc_generic_snes_solvers.pyx +++ b/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 diff --git a/src/underworld3/utilities/boundary_flux.py b/src/underworld3/utilities/boundary_flux.py index 326f7fcc3..6e5d707b2 100644 --- a/src/underworld3/utilities/boundary_flux.py +++ b/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``. diff --git a/tests/parallel/test_1065_boundary_flux_parallel.py b/tests/parallel/test_1065_boundary_flux_parallel.py index 7f20a9636..37038f1e4 100644 --- a/tests/parallel/test_1065_boundary_flux_parallel.py +++ b/tests/parallel/test_1065_boundary_flux_parallel.py @@ -21,6 +21,8 @@ # SERIAL reference: BdIntegral of the flux field over Bottom. `python `. 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}") diff --git a/tests/test_1019_boundary_flux.py b/tests/test_1019_boundary_flux.py index 805f86757..d9d4ce6a3 100644 --- a/tests/test_1019_boundary_flux.py +++ b/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}") diff --git a/tests/test_1020_boundary_flux_integral.py b/tests/test_1020_boundary_flux_integral.py new file mode 100644 index 000000000..a2a052b8a --- /dev/null +++ b/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")