From 4dbd30bd837ff013f7643b326edf012b7a3f552c Mon Sep 17 00:00:00 2001 From: lmoresi Date: Wed, 1 Jul 2026 22:01:19 +1000 Subject: [PATCH 1/3] feat(solvers): general boundary_flux primitive (Consistent Boundary Flux) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Add solver.boundary_flux(boundary, mass="lumped", remove_mean=False, normal=None) and solver.boundary_flux_field(...) to SolverBaseClass, so the Consistent Boundary Flux recovery is a general capability of any solver, not a Stokes-only / topography-specific one. The essential-BC reaction (the assembled interior FEM residual at a Dirichlet node) is the consistent nodal flux; de-smearing it with the boundary mass gives a pointwise surface flux: - scalar diffusion / advection-diffusion -> surface heat flux -k dT/dn (mean = Nusselt) - Stokes -> boundary traction sigma.n (sigma_nn) New pieces: - SolverBaseClass._assemble_volume_reaction(): field-structure-agnostic (self.fields for Stokes, else the single Unknowns.u), globally assembled volume FEM residual. - utilities/boundary_flux.py: equation-independent boundary-node gathering, the lumped/ consistent boundary-mass de-smear (remove_mean flag: keep the mean for a physical flux = Nusselt, remove it for a gauge-free field = dynamic topography), and the field write. Naming: boundary_flux is the primitive; boundary_normal_traction / dynamic_topography (on the rotated-free-slip branch) become thin readings of it (traction; mean-removed topography). remove_mean defaults False here — the mean flux is physical. Validated on a harmonic Poisson (analytic flux): serial corr 1.0000, relL2 2e-4, converging 0.0008->0.0002->0.0004 at res 24/48/96; mean flux == Nusselt (2/sinh pi); np=2 machine-identical. test_1019 (serial). KNOWN LIMITATION (inline TODO): a partition that CUTS a flux boundary (np=4 on a box) has a ~few-% localized error at the cut node — the parallel assembly of the volume FEM residual there is not yet exact. The solve is partition-independent (verified to 1e-12); serial + boundary-uncut partitions are machine-exact. Parallel test gated to np=2 until the cut-node assembly is fixed. Not yet PR'd for that reason. Underworld development team with AI support from Claude Code --- .../cython/petsc_generic_snes_solvers.pyx | 105 +++++++++ src/underworld3/utilities/boundary_flux.py | 223 ++++++++++++++++++ tests/test_1019_boundary_flux.py | 67 ++++++ 3 files changed, 395 insertions(+) create mode 100644 src/underworld3/utilities/boundary_flux.py create mode 100644 tests/test_1019_boundary_flux.py diff --git a/src/underworld3/cython/petsc_generic_snes_solvers.pyx b/src/underworld3/cython/petsc_generic_snes_solvers.pyx index 296e162e..3dfae8aa 100644 --- a/src/underworld3/cython/petsc_generic_snes_solvers.pyx +++ b/src/underworld3/cython/petsc_generic_snes_solvers.pyx @@ -1819,6 +1819,111 @@ class SolverBaseClass(uw_object): return + def _assemble_volume_reaction(self, time=None, verbose=False): + """Globally-assembled FEM VOLUME residual (no boundary terms) in the DM-local + layout, as a numpy array. + + At an essential-BC (Dirichlet) node this residual IS the consistent boundary + reaction — the integrated nodal flux :math:`\\int_\\Gamma (F\\cdot\\hat n)\\phi_i` + (heat flux for a scalar diffusion solve, traction for Stokes). Interior nodes are + ~0. General across scalar / vector / Stokes solvers: the current solution is + gathered from ``self.fields`` when present (Stokes) else the single + ``Unknowns.u`` field, and the per-rank cell residual is summed across ranks + (``localToGlobal`` ADD) so boundary nodes shared across a partition cut carry the + complete flux. + """ + cdef DM dm + cdef Vec xvec + cdef Vec fvec + cdef DM _time_dm_reaction + + self._build(verbose, False, None) + + if time is not None: + if hasattr(time, 'magnitude') or hasattr(time, '_pint_qty'): + t_nd = float(uw.non_dimensionalise(time)) + else: + t_nd = float(time) + _time_dm_reaction = self.dm + UW_DMSetTime(_time_dm_reaction.dm, t_nd) + + self.mesh.update_lvec() + self.dm.setAuxiliaryVec(self.mesh.lvec, None) + self._update_constants() + + gvec = self.dm.getGlobalVec() + xlocal = self.dm.getLocalVec() + flocal = self.dm.getLocalVec() + gvec.setArray(0.0) + xlocal.setArray(0.0) + flocal.setArray(0.0) + + try: + # gather the current solution into the global vector (field-structure agnostic) + if getattr(self, "fields", None): + for name, var in self.fields.items(): + sgvec = gvec.getSubVector(self._subdict[name][0]) + self._subdict[name][1].localToGlobal(var.vec, sgvec) + gvec.restoreSubVector(self._subdict[name][0], sgvec) + else: + _names, _iss, _subdms = self.dm.createFieldDecomposition() + sgvec = gvec.getSubVector(_iss[0]) + _subdms[0].localToGlobal(self.Unknowns.u.vec, sgvec) + gvec.restoreSubVector(_iss[0], sgvec) + + self.dm.globalToLocal(gvec, xlocal) + + dm = self.dm + xvec = xlocal + fvec = flocal + CHKERRQ(DMPlexSNESComputeResidualFEM(dm.dm, xvec.vec, fvec.vec, NULL)) + + # assemble globally (sum shared-node contributions) then scatter back so + # every rank's boundary nodes — including ghosts — hold the complete reaction + # TODO(BUG): a partition that CUTS a flux boundary (e.g. np=4 on a box, bottom + # split across ranks) leaves a ~few-% localized error at the cut node — the + # localToGlobal(ADD) of the DMPlexSNESComputeResidualFEM output is not exact + # there (overlap / owned-cell interaction). Serial + boundary-uncut partitions + # (np<=2 here) are machine-exact; the SOLVE is partition-independent. Fix the + # parallel volume-residual assembly at cut boundary nodes. + gvec.setArray(0.0) + self.dm.localToGlobal(flocal, gvec, addv=PETSc.InsertMode.ADD_VALUES) + self.dm.globalToLocal(gvec, flocal) + return np.array(flocal.array, copy=True) + finally: + self.dm.restoreLocalVec(flocal) + self.dm.restoreLocalVec(xlocal) + self.dm.restoreGlobalVec(gvec) + + def boundary_flux(self, boundary, mass="lumped", remove_mean=False, normal=None): + r"""Consistent boundary flux on ``boundary``, recovered from the essential-BC + reaction of the last solve (the Consistent Boundary Flux method). + + Returns ``(xs, flux)`` with one entry per boundary node on this rank: for a + **scalar** solver the outward normal flux :math:`F\cdot\hat n` (e.g. surface heat + flux :math:`-k\,\partial T/\partial n`, whose boundary mean is the Nusselt + number); for a **vector** solver the traction :math:`\sigma\cdot\hat n` (pass + ``normal`` to get the scalar normal component :math:`\hat n\cdot\sigma\cdot\hat n`). + + ``mass`` de-smears the nodal reaction with the ``"lumped"`` (diagonal, monotone — + no overshoot at a flux jump) or ``"consistent"`` boundary mass. ``remove_mean`` + subtracts the boundary mean — leave ``False`` for a physical flux (the mean is + the Nusselt number); ``True`` gives a gauge-free field (e.g. dynamic topography). + Parallel-safe and partition-independent.""" + from underworld3.utilities.boundary_flux import boundary_flux as _bf + return _bf(self, boundary, mass=mass, remove_mean=remove_mean, normal=normal) + + def boundary_flux_field(self, boundary, field, mass="lumped", + remove_mean=False, scale=1.0, normal=None): + r"""Write the consistent boundary flux (see :meth:`boundary_flux`) onto a scalar + MeshVariable ``field`` at the boundary nodes (interior untouched), multiplied by + ``scale``. This is the field hand-off for downstream machinery (surface heat + flux for coupling, or — with ``remove_mean=True`` and ``scale=-1/(\Delta\rho g)`` + — dynamic topography). Returns ``field``.""" + from underworld3.utilities.boundary_flux import boundary_flux_to_field as _bff + return _bff(self, boundary, field, mass=mass, remove_mean=remove_mean, + scale=scale, normal=normal) + ## Specific to dimensionality diff --git a/src/underworld3/utilities/boundary_flux.py b/src/underworld3/utilities/boundary_flux.py new file mode 100644 index 00000000..b022e063 --- /dev/null +++ b/src/underworld3/utilities/boundary_flux.py @@ -0,0 +1,223 @@ +"""Consistent Boundary Flux (CBF) recovery — general across solvers. + +The residual of the assembled interior (volume) FEM problem, read at an essential-BC +boundary node, is the *consistent* nodal flux there (Gresho et al.). De-smearing that +nodal reaction with the boundary mass gives a pointwise surface flux: + + * scalar diffusion / advection-diffusion → surface heat flux -k dT/dn (Nusselt), + * Stokes → boundary traction sigma.n (sigma_nn). + +This module holds the parts that do not depend on the equation: the boundary-node +gathering, the boundary-mass de-smear (lumped / consistent), and the field hand-off. +The equation-specific bit — extracting the nodal reaction — is the solver method +``_assemble_volume_reaction`` (globally assembled, so shared boundary nodes are complete). + +``mass="lumped"`` (default) uses the diagonal boundary mass: being an M-matrix it cannot +overshoot where the flux jumps (no Gibbs wiggle) and is a purely local division. +``remove_mean=False`` (default) keeps the physical mean flux (the Nusselt number); +set ``remove_mean=True`` for a gauge-free field (e.g. dynamic topography). +""" +import numpy as np +from mpi4py import MPI + + +def _key(c, dim): + return tuple(round(float(t), 9) for t in np.asarray(c).ravel()[:dim]) + + +def _point_coord(dm, dim, cvec, csec, v0, v1, q): + """Coordinate of a DMPlex point (vertex → its coord; higher point → mean of its + closure vertices).""" + if v0 <= q < v1: + return cvec[csec.getOffset(q) // dim] + clo = dm.getTransitiveClosure(q)[0] + verts = [int(c) for c in clo if v0 <= c < v1] + return np.mean([cvec[csec.getOffset(v) // dim] for v in verts], axis=0) + + +def _boundary_field_nodes(solver, boundary, field_id=0): + """DMPlex points carrying `field_id` DOFs on `boundary`, with their coordinates. + Parallel-safe: a rank owning no part of the boundary gets a NULL stratum IS + (guarded); ghost nodes are included (their reaction is completed by the global + assembly in ``_assemble_volume_reaction``).""" + dm = solver.dm + dim = solver.mesh.dim + lsec = dm.getLocalSection() + csec = dm.getCoordinateSection() + cvec = np.asarray(dm.getCoordinatesLocal().array).reshape(-1, dim) + v0, v1 = dm.getDepthStratum(0) + fS, fE = dm.getHeightStratum(1) + bval = [b.value for b in solver.mesh.boundaries if b.name == boundary][0] + sis = dm.getStratumIS(boundary, bval) + if sis is None or sis.handle == 0: + return [], lsec, csec, cvec, v0, v1 + facets = [int(z) for z in sis.getIndices()] + seen = set(); out = [] + for f in facets: + if not (fS <= f < fE): + continue + for q in (int(c) for c in dm.getTransitiveClosure(f)[0]): + if q in seen or lsec.getFieldDof(q, field_id) <= 0: + continue + seen.add(q) + out.append((q, _point_coord(dm, dim, cvec, csec, v0, v1, q))) + return out, lsec, csec, cvec, v0, v1 + + +def _node_normals(solver, boundary, normal, nodes, dm, dim, cvec, csec, v0, v1): + """Per-node outward unit normal (only needed to project a vector reaction). + ``normal`` is None (geometric facet normal), a sympy 1×dim Matrix (analytic, + lambdified), or a constant (dim,) vector.""" + interior_ref = cvec.mean(axis=0) + sym_fn = const = None + if normal is not None: + try: + import sympy + if isinstance(normal, sympy.Matrix): + sym_fn = sympy.lambdify(list(solver.mesh.X), + [normal[0, k] for k in range(dim)], "numpy") + except Exception: + sym_fn = None + if sym_fn is None: + const = np.asarray(normal, dtype=float).ravel() + nmap = {} + coord = {q: c for q, c in nodes} + if normal is None: + # accumulate area-weighted facet normals to the closure nodes + bval = [b.value for b in solver.mesh.boundaries if b.name == boundary][0] + sis = dm.getStratumIS(boundary, bval) + facets = [] if (sis is None or sis.handle == 0) else [int(z) for z in sis.getIndices()] + fS, fE = dm.getHeightStratum(1) + acc = {} + for f in facets: + if not (fS <= f < fE): + continue + _, cent, nrm = dm.computeCellGeometryFVM(f) + ne = np.asarray(nrm, float); ne = ne / (np.linalg.norm(ne) + 1e-30) + if np.dot(ne, np.asarray(cent) - interior_ref) < 0: + ne = -ne + for q in (int(c) for c in dm.getTransitiveClosure(f)[0]): + if q in coord: + acc[q] = acc.get(q, np.zeros(dim)) + ne + for q in coord: + nn = acc.get(q, np.zeros(dim)) + nmap[q] = nn / (np.linalg.norm(nn) + 1e-30) + else: + for q, c in nodes: + ne = np.asarray(sym_fn(*c), float).ravel() if sym_fn is not None else const.copy() + nmap[q] = ne / (np.linalg.norm(ne) + 1e-30) + return nmap + + +def _desmear(solver, boundary, xs, R, mass, remove_mean): + """De-smear per-node reaction loads R (aligned with xs) into a pointwise flux via the + boundary mass, assembled globally by a coordinate-keyed allgather so every rank forms + the identical system. Returns the flux at this rank's local nodes (xs order).""" + dm = solver.dm; dim = solver.mesh.dim; comm = dm.comm.tompi4py() + csec = dm.getCoordinateSection() + cvec = np.asarray(dm.getCoordinatesLocal().array).reshape(-1, dim) + v0, v1 = dm.getDepthStratum(0) + if dim != 2: + # no line-mass geometry yet in 3D → global-mean lumped fallback + tot = comm.allreduce(float(np.sum(R)), op=MPI.SUM) + cnt = comm.allreduce(int(len(R)), op=MPI.SUM) + m = tot / max(cnt, 1) + return np.asarray(R) - (m if remove_mean else 0.0) + + e0, e1 = dm.getDepthStratum(1) + def vcoord(q): return cvec[csec.getOffset(q) // dim] + nodeR = {_key(x, dim): float(r) for x, r in zip(xs, R)} + bval = [b.value for b in solver.mesh.boundaries if b.name == boundary][0] + sis = dm.getStratumIS(boundary, bval) + strat = [] if (sis is None or sis.handle == 0) else [int(z) for z in sis.getIndices()] + local_elems = [] + for e in [q for q in strat if e0 <= q < e1]: + a, b = (int(c) for c in dm.getCone(e)) + cmid = _point_coord(dm, dim, cvec, csec, v0, v1, e) + h = float(np.hypot(*(vcoord(b) - vcoord(a)))) + local_elems.append((_key(vcoord(a), dim), _key(cmid, dim), _key(vcoord(b), dim), h)) + + R_by = {} + for d in comm.allgather(nodeR): + R_by.update(d) + uniq = {} + for lst in comm.allgather(local_elems): + for (ka, km, kb, h) in lst: + uniq[(ka, km, kb)] = h + keys = sorted(R_by.keys()); gi = {k: i for i, k in enumerate(keys)} + n = len(keys); Rg = np.zeros(n) + for k, i in gi.items(): + Rg[i] = R_by[k] + if mass == "lumped": + mL = np.zeros(n) + for (ka, km, kb), h in uniq.items(): + mL[gi[ka]] += h / 6.0; mL[gi[km]] += 2.0 * h / 3.0; mL[gi[kb]] += h / 6.0 + sig = Rg / mL + else: + M = np.zeros((n, n)) + Me = np.array([[4., 2, -1], [2, 16, 2], [-1, 2, 4]]) + for (ka, km, kb), h in uniq.items(): + tri = [gi[ka], gi[km], gi[kb]]; Mh = (h / 30.0) * Me + for ii in range(3): + for jj in range(3): + M[tri[ii], tri[jj]] += Mh[ii, jj] + sig = np.linalg.solve(M, Rg) + if remove_mean: + sig = sig - sig.mean() + return np.array([sig[gi[_key(x, dim)]] for x in xs]) + + +def boundary_flux(solver, boundary, mass="lumped", remove_mean=False, normal=None): + """See ``SolverBaseClass.boundary_flux``. Returns ``(xs, flux)`` for this rank's + boundary nodes; scalar solver → normal flux, vector solver → traction (or its normal + component if ``normal`` is given).""" + dm = solver.dm; dim = solver.mesh.dim + ra = np.asarray(solver._assemble_volume_reaction()).ravel() + nodes, lsec, csec, cvec, v0, v1 = _boundary_field_nodes(solver, boundary, field_id=0) + ncomp = lsec.getFieldComponents(0) + xs = np.array([c for _q, c in nodes]) if nodes else np.zeros((0, dim)) + + if ncomp == 1: + R = np.array([ra[lsec.getFieldOffset(q, 0)] for q, _c in nodes]) if nodes else np.zeros(0) + flux = _desmear(solver, boundary, xs, R, mass, remove_mean) + return xs, flux + + # vector reaction (traction sigma.n at each node) + Rvec = np.array([ra[lsec.getFieldOffset(q, 0):lsec.getFieldOffset(q, 0) + ncomp] + for q, _c in nodes]) if nodes else np.zeros((0, ncomp)) + if normal is not None: + # scalar NORMAL component sigma_nn = n.(sigma.n) + nmap = _node_normals(solver, boundary, normal, nodes, dm, dim, cvec, csec, v0, v1) + Rn = np.array([float(np.dot(nmap[q], Rvec[i])) for i, (q, _c) in enumerate(nodes)]) \ + if nodes else np.zeros(0) + return xs, _desmear(solver, boundary, xs, Rn, mass, remove_mean) + # full traction vector: de-smear each component independently + cols = [_desmear(solver, boundary, xs, Rvec[:, k] if len(Rvec) else np.zeros(0), + mass, remove_mean) for k in range(ncomp)] + return xs, (np.column_stack(cols) if nodes else np.zeros((0, ncomp))) + + +def boundary_flux_to_field(solver, boundary, field, mass="lumped", + remove_mean=False, scale=1.0, normal=None): + """See ``SolverBaseClass.boundary_flux_field``. Writes ``scale * flux`` onto the + scalar MeshVariable ``field`` at the boundary nodes (interior untouched).""" + dim = solver.mesh.dim + xs, flux = boundary_flux(solver, boundary, mass=mass, remove_mean=remove_mean, normal=normal) + fmap = {_key(x, dim): scale * float(f) for x, f in zip(np.asarray(xs), np.asarray(flux).ravel())} + fc = np.asarray(field.coords) + # Write the field ONCE from a local copy: a per-node write to var.data fires the + # write-callback each time, and the boundary-node count differs per rank (a rank may + # own none of the boundary), so per-node writes would desync the callback and hang. + newdata = np.asarray(field.data).copy() + for i in range(fc.shape[0]): + v = fmap.get(_key(fc[i], dim)) + if v is not None: + newdata[i, 0] = v + field.data[...] = newdata + base = getattr(field, "_base_var", field) + if hasattr(base, "_sync_lvec_to_gvec"): + base._sync_lvec_to_gvec() + if hasattr(base, "_canonical_data"): + base._canonical_data = None + solver.mesh._stale_lvec = True + return field diff --git a/tests/test_1019_boundary_flux.py b/tests/test_1019_boundary_flux.py new file mode 100644 index 00000000..4a6bfe6f --- /dev/null +++ b/tests/test_1019_boundary_flux.py @@ -0,0 +1,67 @@ +"""Consistent Boundary Flux (solver.boundary_flux) — scalar surface heat flux. + +The essential-BC reaction of a diffusion solve, de-smeared by the boundary mass, is the +consistent surface flux -k dT/dn (Gresho et al.); its boundary mean is the Nusselt +number. Validated against a harmonic manufactured solution with a known analytic flux. + +Serial (and np=2, where the flux boundary is not cut). NOTE: a partition that CUTS the +flux boundary (e.g. np=4 on this box) currently has a localized error in the parallel +assembly of the volume FEM residual at the cut node — tracked as a follow-up; the solve +itself is partition-independent. So the parallel test is gated to np=2. +""" +import numpy as np +import sympy +import pytest +import underworld3 as uw + +pytestmark = [pytest.mark.level_1, pytest.mark.tier_a] + + +def _heatflux_diagnostics(res=48): + mesh = uw.meshing.StructuredQuadBox( + elementRes=(res, res), minCoords=(0, 0), maxCoords=(1, 1), qdegree=3) + x, y = mesh.X + T = uw.discretisation.MeshVariable("Tbf", mesh, 1, degree=2) + q = uw.discretisation.MeshVariable("qbf", mesh, 1, degree=1, continuous=True) + poisson = uw.systems.Poisson(mesh, u_Field=T) + poisson.constitutive_model = uw.constitutive_models.DiffusionModel + poisson.constitutive_model.Parameters.diffusivity = 1 + poisson.f = 0.0 + poisson.add_dirichlet_bc(0.0, "Bottom") + poisson.add_dirichlet_bc(0.0, "Left") + poisson.add_dirichlet_bc(0.0, "Right") + poisson.add_dirichlet_bc(sympy.Matrix([sympy.sin(sympy.pi * x)]), "Top") + poisson.tolerance = 1e-11 + poisson.petsc_options["snes_type"] = "ksponly" + poisson.solve() + + # general API: boundary flux + field hand-off + xs, flux = poisson.boundary_flux("Bottom") # lumped, mean kept + 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()) + + 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 + + +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) + 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) + assert abs(corr) > 0.999, f"heat flux corr {corr:.4f} too low" + assert relL2 < 0.01, f"heat flux relL2 vs analytic {relL2:.4f} too large" + # physical MEAN flux preserved (Nusselt): |mean| ≈ 2/sinh(pi) + 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 + + +if __name__ == "__main__": + _f, _a, _b = _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}") From 40e8699638e7fff979429d7e660fae1efc381554 Mon Sep 17 00:00:00 2001 From: lmoresi Date: Thu, 2 Jul 2026 07:49:05 +1000 Subject: [PATCH 2/3] fix(boundary-flux): correct parallel reaction assembly at partition-cut boundary nodes MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The boundary flux was wrong (~7% localized) when a partition CUT the flux boundary (np=4 on a box). Root cause: _assemble_volume_reaction hand-rolled the assembly with localToGlobal(ADD) on the DMPlexSNESComputeResidualFEM output, which does not reproduce the true reaction at a cut node. Confirmed against the rock-solid VOLUME integral identity total_flux = integral(grad T . grad w): that is byte-identical at np=1/2/4, while the ADD-assembled reaction diverged at np=4. Fix (no hand-rolled global assembly): the DM has overlap=0, so PETSc's volume residual gives each rank only its OWNED cells' partial contribution at a shared boundary node. Return the raw local residual and assemble the complete reaction by SUMMING each rank's partial across ranks by coordinate — the same rock-solid coordinate gather already used for the boundary mass. This reproduces the volume integral exactly at cut nodes. Validated (harmonic Poisson, analytic flux): surface heat flux corr 1.0000, relL2 2e-4, and BdIntegral of the flux field byte-identical at np=1/2/4. New parallel test test_1065 passes at -n 2 AND -n 4 (boundary cut); serial test_1019 still passes. Underworld development team with AI support from Claude Code --- .../cython/petsc_generic_snes_solvers.pyx | 21 +++-- src/underworld3/utilities/boundary_flux.py | 6 +- .../test_1065_boundary_flux_parallel.py | 79 +++++++++++++++++++ tests/test_1019_boundary_flux.py | 8 +- 4 files changed, 97 insertions(+), 17 deletions(-) create mode 100644 tests/parallel/test_1065_boundary_flux_parallel.py diff --git a/src/underworld3/cython/petsc_generic_snes_solvers.pyx b/src/underworld3/cython/petsc_generic_snes_solvers.pyx index 3dfae8aa..0b6a01d4 100644 --- a/src/underworld3/cython/petsc_generic_snes_solvers.pyx +++ b/src/underworld3/cython/petsc_generic_snes_solvers.pyx @@ -1836,6 +1836,10 @@ class SolverBaseClass(uw_object): cdef Vec xvec cdef Vec fvec cdef DM _time_dm_reaction + cdef PetscFormKey key + cdef IS ccell_is + cdef PetscReal residual_time = 0.0 + cdef PetscReal implicit_form_time = -1.7976931348623157e308 self._build(verbose, False, None) @@ -1878,17 +1882,12 @@ class SolverBaseClass(uw_object): fvec = flocal CHKERRQ(DMPlexSNESComputeResidualFEM(dm.dm, xvec.vec, fvec.vec, NULL)) - # assemble globally (sum shared-node contributions) then scatter back so - # every rank's boundary nodes — including ghosts — hold the complete reaction - # TODO(BUG): a partition that CUTS a flux boundary (e.g. np=4 on a box, bottom - # split across ranks) leaves a ~few-% localized error at the cut node — the - # localToGlobal(ADD) of the DMPlexSNESComputeResidualFEM output is not exact - # there (overlap / owned-cell interaction). Serial + boundary-uncut partitions - # (np<=2 here) are machine-exact; the SOLVE is partition-independent. Fix the - # parallel volume-residual assembly at cut boundary nodes. - gvec.setArray(0.0) - self.dm.localToGlobal(flocal, gvec, addv=PETSc.InsertMode.ADD_VALUES) - self.dm.globalToLocal(gvec, flocal) + # Return the RAW local residual: each rank has computed its OWNED cells' + # contribution to its local nodes (the DM has overlap=0, so a boundary node + # shared across a partition cut holds only this rank's PARTIAL contribution). + # The caller assembles the complete reaction by summing these partials across + # ranks by coordinate (boundary_flux._desmear), consistent with the boundary- + # mass gather — this reproduces the rock-solid volume integral at cut nodes. return np.array(flocal.array, copy=True) finally: self.dm.restoreLocalVec(flocal) diff --git a/src/underworld3/utilities/boundary_flux.py b/src/underworld3/utilities/boundary_flux.py index b022e063..122515c5 100644 --- a/src/underworld3/utilities/boundary_flux.py +++ b/src/underworld3/utilities/boundary_flux.py @@ -137,9 +137,13 @@ def vcoord(q): return cvec[csec.getOffset(q) // dim] h = float(np.hypot(*(vcoord(b) - vcoord(a)))) local_elems.append((_key(vcoord(a), dim), _key(cmid, dim), _key(vcoord(b), dim), h)) + # SUM the nodal reaction across ranks by coordinate: with overlap=0 a boundary node + # shared across a partition cut holds only each rank's partial cell contribution, so + # summing them assembles the complete reaction (matches the rock-solid volume integral). R_by = {} for d in comm.allgather(nodeR): - R_by.update(d) + for k, v in d.items(): + R_by[k] = R_by.get(k, 0.0) + v uniq = {} for lst in comm.allgather(local_elems): for (ka, km, kb, h) in lst: diff --git a/tests/parallel/test_1065_boundary_flux_parallel.py b/tests/parallel/test_1065_boundary_flux_parallel.py new file mode 100644 index 00000000..f3cb42a8 --- /dev/null +++ b/tests/parallel/test_1065_boundary_flux_parallel.py @@ -0,0 +1,79 @@ +"""Parallel regression: solver.boundary_flux is partition-independent, including when +the flux boundary is CUT by the partition (np=4 on this box splits the bottom). + +The nodal reaction is PETSc's volume FEM residual (rock-solid); a boundary node shared +across a partition cut holds only each rank's partial contribution (DM overlap=0), and the +complete reaction is assembled by summing the partials across ranks by coordinate. This +test checks that the surface heat flux — and its BdIntegral — reproduce the serial +reference at np=2 and np=4. + +Run: + mpirun -n 2 python -m pytest --with-mpi tests/parallel/test_1065_boundary_flux_parallel.py + mpirun -n 4 python -m pytest --with-mpi tests/parallel/test_1065_boundary_flux_parallel.py +""" +import numpy as np +import sympy +import pytest +import underworld3 as uw + +pytestmark = [pytest.mark.mpi(min_size=2), pytest.mark.timeout(180)] + +# SERIAL reference: BdIntegral of the flux field over Bottom. `python `. +GOLDEN_BDFLUX = -1.731543e-01 + + +def _flux_diagnostics(res=48): + mesh = uw.meshing.StructuredQuadBox( + elementRes=(res, res), minCoords=(0, 0), maxCoords=(1, 1), qdegree=3) + x, y = mesh.X + T = uw.discretisation.MeshVariable("Tp", mesh, 1, degree=2) + q = uw.discretisation.MeshVariable("qp", mesh, 1, degree=1, continuous=True) + poisson = uw.systems.Poisson(mesh, u_Field=T) + poisson.constitutive_model = uw.constitutive_models.DiffusionModel + poisson.constitutive_model.Parameters.diffusivity = 1 + poisson.f = 0.0 + poisson.add_dirichlet_bc(0.0, "Bottom") + poisson.add_dirichlet_bc(0.0, "Left") + poisson.add_dirichlet_bc(0.0, "Right") + poisson.add_dirichlet_bc(sympy.Matrix([sympy.sin(sympy.pi * x)]), "Top") + poisson.tolerance = 1e-11 + poisson.petsc_options["snes_type"] = "ksponly" + poisson.solve() + + 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()) + + # gather + dedup for a whole-boundary relL2 vs analytic (on rank 0, then bcast) + comm = uw.mpi.comm + gx = comm.gather(np.asarray(xs).reshape(-1, 2), root=0) + gf = comm.gather(np.asarray(flux).reshape(-1), root=0) + relL2 = None + if uw.mpi.rank == 0: + seen = {} + for xb, fb in zip(gx, gf): + for xc, fc in zip(xb, fb): + seen[(round(float(xc[0]), 9),)] = (xc[0], fc) + X = np.array([v[0] for v in seen.values()]) + F = np.array([v[1] for v in seen.values()]) + q_an = np.pi * np.sin(np.pi * X) / np.sinh(np.pi) + 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) + + +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) + 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}" + + +if __name__ == "__main__": + _b, _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 4a6bfe6f..644a8af7 100644 --- a/tests/test_1019_boundary_flux.py +++ b/tests/test_1019_boundary_flux.py @@ -3,11 +3,9 @@ The essential-BC reaction of a diffusion solve, de-smeared by the boundary mass, is the consistent surface flux -k dT/dn (Gresho et al.); its boundary mean is the Nusselt number. Validated against a harmonic manufactured solution with a known analytic flux. - -Serial (and np=2, where the flux boundary is not cut). NOTE: a partition that CUTS the -flux boundary (e.g. np=4 on this box) currently has a localized error in the parallel -assembly of the volume FEM residual at the cut node — tracked as a follow-up; the solve -itself is partition-independent. So the parallel test is gated to np=2. +The nodal reaction is PETSc's (rock-solid) volume FEM residual; its complete value at a +partition-cut boundary node is assembled by summing each rank's partial contribution by +coordinate, so the flux is partition-independent (see the parallel test). """ import numpy as np import sympy From 787c0bde9f1757f90838e93f9d4a82f99714fd92 Mon Sep 17 00:00:00 2001 From: lmoresi Date: Thu, 2 Jul 2026 10:03:08 +1000 Subject: [PATCH 3/3] =?UTF-8?q?fix(boundary-flux):=20address=20Copilot=20r?= =?UTF-8?q?eview=20=E2=80=94=20robust=20boundary=20label,=20fail-fast,=20d?= =?UTF-8?q?ocs?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit - Use the consolidated "UW_Boundaries" DM label (distinguished by value) for boundary strata instead of a per-boundary label named like the boundary. The named per-boundary labels do not survive mesh adaptation (only UW_Boundaries is rebuilt), so the old lookup could return an empty stratum and silently give empty flux on adapted meshes. New _boundary_stratum_is() helper; raises a clear error for an unknown boundary name (was an unhelpful IndexError). Fixes all three lookup sites. - boundary_flux_field: fail fast when the target is a scalar MeshVariable but the flux is a vector (a vector solver with normal=None), instead of silently pairing a flattened vector with the nodes. - Correct docstrings (module, _boundary_field_nodes, _assemble_volume_reaction) to describe the actual parallel path: the reaction is the RAW per-rank volume residual, assembled at cut nodes by coordinate-summing in _desmear (not a global localToGlobal). - Note the consistent-mass de-smear is a dense O(n^3) solve; lumped (default) is O(n). Serial test_1019 and parallel test_1065 (-n 4) still pass. Underworld development team with AI support from Claude Code --- .../cython/petsc_generic_snes_solvers.pyx | 14 ++++-- src/underworld3/utilities/boundary_flux.py | 48 +++++++++++++++---- 2 files changed, 47 insertions(+), 15 deletions(-) diff --git a/src/underworld3/cython/petsc_generic_snes_solvers.pyx b/src/underworld3/cython/petsc_generic_snes_solvers.pyx index 0b6a01d4..d5188c4b 100644 --- a/src/underworld3/cython/petsc_generic_snes_solvers.pyx +++ b/src/underworld3/cython/petsc_generic_snes_solvers.pyx @@ -1820,17 +1820,21 @@ class SolverBaseClass(uw_object): return def _assemble_volume_reaction(self, time=None, verbose=False): - """Globally-assembled FEM VOLUME residual (no boundary terms) in the DM-local - layout, as a numpy array. + """RAW per-rank FEM VOLUME residual (no boundary terms) in the DM-local layout, + as a numpy array. At an essential-BC (Dirichlet) node this residual IS the consistent boundary reaction — the integrated nodal flux :math:`\\int_\\Gamma (F\\cdot\\hat n)\\phi_i` (heat flux for a scalar diffusion solve, traction for Stokes). Interior nodes are ~0. General across scalar / vector / Stokes solvers: the current solution is gathered from ``self.fields`` when present (Stokes) else the single - ``Unknowns.u`` field, and the per-rank cell residual is summed across ranks - (``localToGlobal`` ADD) so boundary nodes shared across a partition cut carry the - complete flux. + ``Unknowns.u`` field. + + NOTE: the returned array is NOT globally assembled. The DM has overlap=0, so each + rank computes only its OWNED cells' contribution; a boundary node shared across a + partition cut therefore holds only this rank's PARTIAL reaction. The complete + reaction is assembled by the caller (``utilities.boundary_flux._desmear``) by + SUMMING each rank's partial by coordinate — not by a hand-rolled localToGlobal. """ cdef DM dm cdef Vec xvec diff --git a/src/underworld3/utilities/boundary_flux.py b/src/underworld3/utilities/boundary_flux.py index 122515c5..fde9e725 100644 --- a/src/underworld3/utilities/boundary_flux.py +++ b/src/underworld3/utilities/boundary_flux.py @@ -10,7 +10,10 @@ This module holds the parts that do not depend on the equation: the boundary-node gathering, the boundary-mass de-smear (lumped / consistent), and the field hand-off. The equation-specific bit — extracting the nodal reaction — is the solver method -``_assemble_volume_reaction`` (globally assembled, so shared boundary nodes are complete). +``_assemble_volume_reaction``, which returns each rank's RAW (per-rank) volume FEM +residual; the complete reaction at a boundary node shared across a partition cut is then +assembled here in ``_desmear`` by SUMMING each rank's partial contribution by coordinate +(the same rock-solid gather used for the boundary mass — no hand-rolled global assembly). ``mass="lumped"`` (default) uses the diagonal boundary mass: being an M-matrix it cannot overshoot where the flux jumps (no Gibbs wiggle) and is a purely local division. @@ -25,6 +28,22 @@ def _key(c, dim): return tuple(round(float(t), 9) for t in np.asarray(c).ravel()[:dim]) +def _boundary_stratum_is(dm, mesh, boundary): + """Facet stratum IS for ``boundary`` via the CONSOLIDATED ``UW_Boundaries`` label + (boundaries are distinguished by value on this one label — the per-boundary labels + named like the boundary do not survive mesh adaptation; only ``UW_Boundaries`` is + rebuilt). Returns None if this rank owns no part of the boundary. Raises a clear + error for an unknown boundary name.""" + match = [b.value for b in mesh.boundaries if b.name == boundary] + if not match: + raise ValueError( + f"Unknown boundary {boundary!r}; known: {[b.name for b in mesh.boundaries]}") + label = dm.getLabel("UW_Boundaries") + if label is None: + return None + return label.getStratumIS(match[0]) + + def _point_coord(dm, dim, cvec, csec, v0, v1, q): """Coordinate of a DMPlex point (vertex → its coord; higher point → mean of its closure vertices).""" @@ -38,8 +57,8 @@ def _point_coord(dm, dim, cvec, csec, v0, v1, q): def _boundary_field_nodes(solver, boundary, field_id=0): """DMPlex points carrying `field_id` DOFs on `boundary`, with their coordinates. Parallel-safe: a rank owning no part of the boundary gets a NULL stratum IS - (guarded); ghost nodes are included (their reaction is completed by the global - assembly in ``_assemble_volume_reaction``).""" + (guarded); ghost/shared nodes are included and their partial per-rank reactions are + summed by coordinate in ``_desmear`` to form the complete reaction.""" dm = solver.dm dim = solver.mesh.dim lsec = dm.getLocalSection() @@ -47,8 +66,7 @@ def _boundary_field_nodes(solver, boundary, field_id=0): cvec = np.asarray(dm.getCoordinatesLocal().array).reshape(-1, dim) v0, v1 = dm.getDepthStratum(0) fS, fE = dm.getHeightStratum(1) - bval = [b.value for b in solver.mesh.boundaries if b.name == boundary][0] - sis = dm.getStratumIS(boundary, bval) + sis = _boundary_stratum_is(dm, solver.mesh, boundary) if sis is None or sis.handle == 0: return [], lsec, csec, cvec, v0, v1 facets = [int(z) for z in sis.getIndices()] @@ -84,8 +102,7 @@ def _node_normals(solver, boundary, normal, nodes, dm, dim, cvec, csec, v0, v1): coord = {q: c for q, c in nodes} if normal is None: # accumulate area-weighted facet normals to the closure nodes - bval = [b.value for b in solver.mesh.boundaries if b.name == boundary][0] - sis = dm.getStratumIS(boundary, bval) + sis = _boundary_stratum_is(dm, solver.mesh, boundary) facets = [] if (sis is None or sis.handle == 0) else [int(z) for z in sis.getIndices()] fS, fE = dm.getHeightStratum(1) acc = {} @@ -127,8 +144,7 @@ def _desmear(solver, boundary, xs, R, mass, remove_mean): e0, e1 = dm.getDepthStratum(1) def vcoord(q): return cvec[csec.getOffset(q) // dim] nodeR = {_key(x, dim): float(r) for x, r in zip(xs, R)} - bval = [b.value for b in solver.mesh.boundaries if b.name == boundary][0] - sis = dm.getStratumIS(boundary, bval) + sis = _boundary_stratum_is(dm, solver.mesh, boundary) strat = [] if (sis is None or sis.handle == 0) else [int(z) for z in sis.getIndices()] local_elems = [] for e in [q for q in strat if e0 <= q < e1]: @@ -158,6 +174,9 @@ def vcoord(q): return cvec[csec.getOffset(q) // dim] mL[gi[ka]] += h / 6.0; mL[gi[km]] += 2.0 * h / 3.0; mL[gi[kb]] += h / 6.0 sig = Rg / mL else: + # consistent P2 line mass — a dense (n×n) solve in the number of boundary nodes + # (O(n^3)); fine for a 1D boundary (n ~ resolution) but prefer the default lumped + # (O(n), monotone) for very large boundaries. M = np.zeros((n, n)) Me = np.array([[4., 2, -1], [2, 16, 2], [-1, 2, 4]]) for (ka, km, kb), h in uniq.items(): @@ -207,7 +226,16 @@ def boundary_flux_to_field(solver, boundary, field, mass="lumped", scalar MeshVariable ``field`` at the boundary nodes (interior untouched).""" dim = solver.mesh.dim xs, flux = boundary_flux(solver, boundary, mass=mass, remove_mean=remove_mean, normal=normal) - fmap = {_key(x, dim): scale * float(f) for x, f in zip(np.asarray(xs), np.asarray(flux).ravel())} + flux = np.asarray(flux) + # a SCALAR field can only hold a scalar flux — a vector solver returns a per-node + # traction VECTOR unless a `normal` is given to project it. Fail fast rather than + # silently pairing the flattened vector with the nodes. + if flux.ndim > 1 and flux.shape[1] != 1: + raise ValueError( + "boundary_flux_field target is a scalar MeshVariable but boundary_flux " + "returned a vector (traction). Pass normal= to project onto the normal " + "component, or use boundary_flux() directly for the full vector.") + fmap = {_key(x, dim): scale * float(f) for x, f in zip(np.asarray(xs), flux.ravel())} fc = np.asarray(field.coords) # Write the field ONCE from a local copy: a per-node write to var.data fires the # write-callback each time, and the boundary-node count differs per rank (a rank may