diff --git a/src/underworld3/cython/petsc_generic_snes_solvers.pyx b/src/underworld3/cython/petsc_generic_snes_solvers.pyx index 8d07446b2..3b9c803fd 100644 --- a/src/underworld3/cython/petsc_generic_snes_solvers.pyx +++ b/src/underworld3/cython/petsc_generic_snes_solvers.pyx @@ -3070,8 +3070,10 @@ class SolverBaseClass(uw_object): ``normal`` to get the scalar normal component :math:`\hat n\cdot\sigma\cdot\hat n`). ``mass`` de-smears the nodal reaction with ``"lumped"`` or ``"consistent"`` - boundary mass. ``"auto"`` (default) selects lumped recovery for 2D traces and - 3D P1 triangles, and the required consistent solve for 3D P2 triangles. + boundary mass. ``"auto"`` (default) selects lumped recovery for 2D P1/P2 + traces and 3D P1 triangles, and the consistent solve for 3D P2 triangles and + 2D traces of degree >= 3 (where row-sum lumping is respectively invalid and + only O(h) pointwise). ``remove_mean`` subtracts the boundary mean — leave ``False`` for a physical flux (the mean is the Nusselt number); ``True`` gives a gauge-free field. diff --git a/src/underworld3/utilities/boundary_flux.py b/src/underworld3/utilities/boundary_flux.py index 6142ec692..0f3cfd10c 100644 --- a/src/underworld3/utilities/boundary_flux.py +++ b/src/underworld3/utilities/boundary_flux.py @@ -15,10 +15,14 @@ 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="auto"`` (default) uses a diagonal lumped mass where the trace basis admits -positive row sums (the 2D P2 line trace and the 3D P1 triangle trace), and the consistent -mass otherwise. A 3D P2 triangle has exactly zero row sum at every vertex, so its -pointwise recovery requires the consistent surface-mass solve. +``mass="auto"`` (default) uses a diagonal lumped mass where lumping is pointwise +sound (the 2D P1/P2 line traces and the 3D P1 triangle trace), and the consistent +mass otherwise. A 3D P2 triangle has exactly zero row sum at every vertex; a 2D trace +of degree ≥ 3 places its edge-interior nodes asymmetrically (Gauss-Jacobi), where +row-sum lumping is only O(h) pointwise — both take the consistent solve. A degree ≥ 3 +trace also carries several interpolation nodes per edge point, and each keeps its own +coordinate — keying both by the edge's single coordinate silently collapsed them +(issue #459). ``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). """ @@ -79,11 +83,88 @@ def _point_coord(dm, dim, cvec, csec, v0, v1, q): return np.mean([cvec[csec.getOffset(v) // dim] for v in verts], axis=0) +def _trace_interior_coords(solver, degree): + """Coordinates of the EDGE-INTERIOR interpolation nodes of a continuous Lagrange + field of ``degree``, keyed by DMPlex edge point: ``{edge: (degree-1, cdim) array}`` + in section slot order. A degree-3 trace carries TWO nodes per edge; the coordinate + section stores only one coordinate per point, so these must be built by + interpolating the mesh coordinate field into a matching-degree space (issue #459). + The space is created exactly as UW3 creates every field FE (``createDefault``, + ``node_endpoints=False`` — see ``Mesh._get_coords_for_basis``), so the per-point + node ordering matches the solver field's section by construction. + COLLECTIVE on the DM's communicator — every rank must call this, boundary or not.""" + from petsc4py import PETSc + + mesh = solver.mesh + cdim = mesh.cdim + dmold = solver.dm.getCoordinateDM() + dmold.createDS() + dmnew = dmold.clone() + prefix = "cbf_trace_coord_" + options = PETSc.Options() + options[prefix + "petscspace_degree"] = degree + options[prefix + "petscdualspace_lagrange_continuity"] = True + options[prefix + "petscdualspace_lagrange_node_endpoints"] = False + fe = PETSc.FE().createDefault( + mesh.dim, cdim, mesh.isSimplex, mesh.qdegree, prefix, PETSc.COMM_SELF) + dmnew.setField(0, fe) + dmnew.createDS() + mat_interp, vec_scale = dmold.createInterpolation(dmnew) + coords_new_g = dmnew.getGlobalVec() + coords_new_l = dmnew.getLocalVec() + mat_interp.mult(solver.dm.getCoordinates(), coords_new_g) + dmnew.globalToLocal(coords_new_g, coords_new_l) + arr = np.asarray(coords_new_l.array).reshape(-1, cdim).copy() + sec = dmnew.getLocalSection() + e0, e1 = dmnew.getDepthStratum(1) + out = {} + for e in range(e0, e1): + ndof = sec.getDof(e) + if ndof > 0: + row = sec.getOffset(e) // cdim + out[e] = arr[row: row + ndof // cdim] + dmnew.restoreGlobalVec(coords_new_g) + dmnew.restoreLocalVec(coords_new_l) + mat_interp.destroy() + if vec_scale is not None: + vec_scale.destroy() + fe.destroy() + dmnew.destroy() + return out + + +def _line_mass_1d(ts): + """Consistent 1-D line-element mass per unit length for a Lagrange basis with + nodes at parameters ``ts`` in [0, 1]: ``M_ij = ∫ L_i L_j dt`` (Gauss–Legendre, + exact for the polynomial integrand). The lumped row sums are ``∫ L_i`` by + partition of unity; their positivity is checked where the lumped path uses them.""" + t = np.asarray(ts, dtype=float) + if np.min(np.diff(np.sort(t))) < 1e-12: + raise RuntimeError( + "Line-trace interpolation nodes are not distinct — the per-node " + "coordinate build is inconsistent with the field layout (issue #459).") + degree = len(t) - 1 + xq, wq = np.polynomial.legendre.leggauss(degree + 1) + xq = 0.5 * (xq + 1.0) + wq = 0.5 * wq + L = np.ones((len(t), len(xq))) + for i, ti in enumerate(t): + for j, tj in enumerate(t): + if i != j: + L[i] *= (xq - tj) / (ti - tj) + return (L * wq) @ L.T + + 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/shared nodes are included and their partial per-rank reactions are - summed by coordinate in ``_desmear`` to form the complete reaction.""" + """Interpolation nodes carrying `field_id` DOFs on `boundary`, one entry per NODE + as ``(point, slot, coord)``. A DMPlex point can carry several nodes — a degree-3 + trace has two edge-interior nodes per edge point — and each keeps its OWN + coordinate: keying both by the point's single coordinate collapses them in the + de-smear and silently drops reactions (issue #459). + Parallel-safe and COLLECTIVE (the per-node coordinate build interpolates the mesh + coordinate field, and whether it is needed is agreed globally): a rank owning no + part of the boundary still participates; ghost/shared nodes are included and their + partial per-rank reactions are summed by coordinate in ``_desmear``.""" dm = solver.dm dim = solver.mesh.dim lsec = dm.getLocalSection() @@ -91,20 +172,41 @@ 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) + ncomp = lsec.getFieldComponents(field_id) sis = _boundary_stratum_is(dm, solver.mesh, boundary) - if not (sis and sis.getSize() > 0): - return [], lsec, csec, cvec, v0, v1 - facets = [int(z) for z in sis.getIndices()] - seen = set(); out = [] + facets = [] if not (sis and sis.getSize() > 0) else [ + int(z) for z in sis.getIndices() if fS <= int(z) < fE] + seen = set(); points = [] 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: + if q in seen: + continue + fdof = lsec.getFieldDof(q, field_id) + if fdof <= 0: continue seen.add(q) - out.append((q, _point_coord(dm, dim, cvec, csec, v0, v1, q))) - return out, lsec, csec, cvec, v0, v1 + if fdof % ncomp: + raise RuntimeError( + f"Field {field_id} carries {fdof} DOFs at point {q} with " + f"{ncomp} components — not a nodal (Lagrange) layout.") + points.append((q, fdof // ncomp)) + nnodes_max = dm.comm.tompi4py().allreduce( + max((m for _q, m in points), default=1), op=MPI.MAX) + interior = _trace_interior_coords(solver, nnodes_max + 1) if nnodes_max >= 2 else {} + out = [] + for q, m in points: + if m == 1: + out.append((q, 0, _point_coord(dm, dim, cvec, csec, v0, v1, q))) + else: + if q not in interior: + # multi-node points other than mesh edges (e.g. the interior nodes of + # a degree-4 face in 3D) have no per-node coordinate build yet + raise NotImplementedError( + f"Boundary point {q} carries {m} interpolation nodes but only " + "edge-interior nodes have per-node coordinates (issue #459).") + for slot, xc in enumerate(np.asarray(interior[q], dtype=float)[:m]): + out.append((q, slot, xc[:dim])) + return out, lsec, csec, cvec, v0, v1, interior def _node_normals(solver, boundary, normal, nodes, dm, dim, cvec, csec, v0, v1): @@ -124,9 +226,10 @@ def _node_normals(solver, boundary, normal, nodes, dm, dim, cvec, csec, v0, v1): if sym_fn is None: const = np.asarray(normal, dtype=float).ravel() nmap = {} - coord = {q: c for q, c in nodes} + pts = {q for q, _s, _c in nodes} if normal is None: - # accumulate area-weighted facet normals to the closure nodes + # accumulate area-weighted facet normals to the closure points; every node of + # a point (e.g. both P3 edge-interior nodes) shares its point's facet normal sis = _boundary_stratum_is(dm, solver.mesh, boundary) facets = [] if not (sis and sis.getSize() > 0) else [int(z) for z in sis.getIndices()] fS, fE = dm.getHeightStratum(1) @@ -139,19 +242,33 @@ def _node_normals(solver, boundary, normal, nodes, dm, dim, cvec, csec, v0, v1): 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: + if q in pts: acc[q] = acc.get(q, np.zeros(dim)) + ne - for q in coord: + for q, s, _c in nodes: nn = acc.get(q, np.zeros(dim)) - nmap[q] = nn / (np.linalg.norm(nn) + 1e-30) + nmap[(q, s)] = nn / (np.linalg.norm(nn) + 1e-30) else: - for q, c in nodes: + for q, s, 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) + nmap[(q, s)] = ne / (np.linalg.norm(ne) + 1e-30) return nmap -def _desmear(solver, boundary, xs, R, mass, remove_mean, partial_reaction=True): +def _node_reactions(xs, R, dim, boundary): + """Coordinate-keyed nodal reactions, refusing the #459 collapse: two reaction + nodes sharing one coordinate key would silently overwrite each other.""" + nodeR = {_key(x, dim): float(r) for x, r in zip(xs, R)} + if len(nodeR) != len(xs): + raise RuntimeError( + f"{len(xs)} boundary reaction nodes on {boundary!r} collapse onto " + f"{len(nodeR)} distinct coordinate keys — per-node coordinates are not " + "distinct (issue #459: a multi-node trace point needs true interpolation-" + "node coordinates, not the point's single coordinate).") + return nodeR + + +def _desmear(solver, boundary, xs, R, mass, remove_mean, partial_reaction=True, + edge_node_coords=None): """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). @@ -160,7 +277,11 @@ def _desmear(solver, boundary, xs, R, mass, remove_mean, partial_reaction=True): reconciled across ranks: ``True`` (default) SUMS each rank's contribution — correct when R is the RAW per-rank volume residual (``boundary_flux``, DM overlap=0); ``False`` OVERWRITES (all ranks already agree) — correct when R comes from an ASSEMBLED global - operator, e.g. the rotated free-slip reaction ``Q(A·u − b)`` (``rotated_bc``).""" + operator, e.g. the rotated free-slip reaction ``Q(A·u − b)`` (``rotated_bc``). + + ``edge_node_coords`` (2D, trace degree ≥ 3 only) is the ``_trace_interior_coords`` + map of per-edge interior node coordinates; ``boundary_flux`` passes the one it built + so the element keys match ``xs`` exactly. ``None`` builds it on demand (collective).""" dm = solver.dm; dim = solver.mesh.dim; comm = dm.comm.tompi4py() csec = dm.getCoordinateSection() cvec = np.asarray(dm.getCoordinatesLocal().array).reshape(-1, dim) @@ -177,7 +298,7 @@ def _desmear(solver, boundary, xs, R, mass, remove_mean, partial_reaction=True): def coord(q): return _point_coord(dm, dim, cvec, csec, v0, v1, q) - nodeR = {_key(x, dim): float(r) for x, r in zip(xs, R)} + nodeR = _node_reactions(xs, R, dim, boundary) sis = _boundary_stratum_is(dm, solver.mesh, boundary) facets = [] if not (sis and sis.getSize() > 0) else [ int(q) for q in sis.getIndices() if f0 <= int(q) < f1 @@ -363,20 +484,44 @@ def value_at(x): raise NotImplementedError( f"Boundary-flux recovery is not implemented for mesh dimension {dim}." ) - if mass in ("auto", "p1"): - mass = "lumped" # 2D: the lumped line mass is sound + if mass == "p1": + mass = "lumped" # P1 consumers read vertex values; vertex lumping is sound + lsec = dm.getLocalSection() + ncomp = lsec.getFieldComponents(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)} + nodeR = _node_reactions(xs, R, dim, boundary) sis = _boundary_stratum_is(dm, solver.mesh, boundary) strat = [] if not (sis and sis.getSize() > 0) else [int(z) for z in sis.getIndices()] + edges = [q for q in strat if e0 <= q < e1] + # Interior-node count per edge from the SECTION (structural — the trace order is + # never sniffed from coordinate keys): 0 → P1, 1 → P2, m → degree m+1. Whether the + # per-node coordinate build is needed is agreed globally (it is collective). + if edge_node_coords is None: + m_max = comm.allreduce( + max((lsec.getFieldDof(e, 0) // ncomp for e in edges), default=0), op=MPI.MAX) + edge_node_coords = _trace_interior_coords(solver, m_max + 1) if m_max >= 2 else {} local_elems = [] - for e in [q for q in strat if e0 <= q < e1]: + for e in edges: + m = lsec.getFieldDof(e, 0) // ncomp 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)) + xa, xb = vcoord(a), vcoord(b) + h = float(np.hypot(*(xb - xa))) + if m >= 2: + xin = np.asarray(edge_node_coords[e], dtype=float)[:m, :dim] + keys_e = (_key(xa, dim), *(_key(x, dim) for x in xin), _key(xb, dim)) + # node parameters along the edge chord, measured from the actual node + # coordinates so the element mass matches the basis in effect + ts_e = (0.0, *(float(np.dot(x - xa, xb - xa) / (h * h)) for x in xin), 1.0) + elif m == 1: + cmid = _point_coord(dm, dim, cvec, csec, v0, v1, e) + keys_e = (_key(xa, dim), _key(cmid, dim), _key(xb, dim)) + ts_e = (0.0, 0.5, 1.0) + else: + keys_e = (_key(xa, dim), _key(xb, dim)) + ts_e = (0.0, 1.0) + local_elems.append((keys_e, ts_e, h)) # Reconcile the nodal reaction across ranks by coordinate. partial_reaction=True: SUM # (raw per-rank residual, DM overlap=0, so a cut node holds only each rank's partial @@ -389,31 +534,47 @@ def vcoord(q): return cvec[csec.getOffset(q) // dim] R_by[k] = (R_by.get(k, 0.0) + v) if partial_reaction else v uniq = {} for lst in comm.allgather(local_elems): - for (ka, km, kb, h) in lst: - uniq[(ka, km, kb)] = h + for keys_e, ts_e, h in lst: + uniq[keys_e] = (ts_e, 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] - # Trace order from the DATA, not an assumption: a P2 trace has a reaction DOF at - # every edge point (the midpoint key is in R_by), a P1 trace has vertices only. - # Assembling P2 line masses against a P1 trace died with a bare KeyError on the - # missing midpoint (issue #413). A mix means the field layout is inconsistent - # with the trace — raise rather than guess. - mids_present = [km in R_by for (ka, km, kb) in uniq] - if mids_present and any(mids_present) != all(mids_present): + # A mixed trace (field DOFs on only part of the boundary's edges) means the field + # layout is inconsistent with the trace — raise rather than guess (issue #413). + orders = {len(keys_e) for keys_e in uniq} + if len(orders) > 1: raise NotImplementedError( - "2D boundary-flux recovery found edge-midpoint reactions on only part " - "of the boundary; mixed P1/P2 traces are not supported." + "2D boundary-flux recovery found different trace orders on parts of " + "the boundary; mixed traces are not supported." ) - trace_order = 2 if (mids_present and mids_present[0]) else 1 + if mass == "auto": + # The row-sum lumped de-smear is pointwise-exact for a linear flux only on + # the symmetric P1/P2 node layouts. A degree >= 3 trace places its interior + # nodes asymmetrically within each edge (PETSc's Gauss-Jacobi nodes), where + # lumping is only O(h) pointwise — the consistent line mass is exact there + # (up to the documented corner mixing, which it spreads over ~one element). + mass = "lumped" if max(orders, default=2) <= 3 else "consistent" + missing = {k for keys_e in uniq for k in keys_e if k not in gi} + if missing: + raise RuntimeError( + f"{len(missing)} trace nodes on {boundary!r} have no reaction entry — " + "the caller's node list does not cover the trace's interpolation nodes.") if mass == "lumped": mL = np.zeros(n) - for (ka, km, kb), h in uniq.items(): - if trace_order == 2: - mL[gi[ka]] += h / 6.0; mL[gi[km]] += 2.0 * h / 3.0; mL[gi[kb]] += h / 6.0 + for keys_e, (ts_e, h) in uniq.items(): + if len(ts_e) == 2: + w = (0.5, 0.5) + elif len(ts_e) == 3: + w = (1.0 / 6.0, 2.0 / 3.0, 1.0 / 6.0) else: - mL[gi[ka]] += h / 2.0; mL[gi[kb]] += h / 2.0 + w = _line_mass_1d(ts_e).sum(axis=1) # ∫ L_i by partition of unity + if np.any(w <= 0.0): + raise ValueError( + f"The degree-{len(ts_e) - 1} line trace has non-positive " + "lumped row sums; use mass='consistent'.") + for kk, wk in zip(keys_e, w): + mL[gi[kk]] += h * wk sig = Rg / mL else: # consistent line mass — a dense (n×n) solve in the number of boundary nodes @@ -422,14 +583,15 @@ def vcoord(q): return cvec[csec.getOffset(q) // dim] M = np.zeros((n, n)) Me2 = np.array([[4., 2, -1], [2, 16, 2], [-1, 2, 4]]) Me1 = np.array([[2., 1], [1, 2]]) - for (ka, km, kb), h in uniq.items(): - if trace_order == 2: - nodes = [gi[ka], gi[km], gi[kb]]; Mh = (h / 30.0) * Me2 + for keys_e, (ts_e, h) in uniq.items(): + if len(ts_e) == 3: + Mh = (h / 30.0) * Me2 + elif len(ts_e) == 2: + Mh = (h / 6.0) * Me1 else: - nodes = [gi[ka], gi[kb]]; Mh = (h / 6.0) * Me1 - for ii in range(len(nodes)): - for jj in range(len(nodes)): - M[nodes[ii], nodes[jj]] += Mh[ii, jj] + Mh = h * _line_mass_1d(ts_e) + idx = [gi[k] for k in keys_e] + M[np.ix_(idx, idx)] += Mh sig = np.linalg.solve(M, Rg) if remove_mean: sig = sig - sig.mean() @@ -442,27 +604,36 @@ def boundary_flux(solver, boundary, mass="auto", remove_mean=False, normal=None) 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) + nodes, lsec, csec, cvec, v0, v1, edge_nodes = _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)) + xs = np.array([c for _q, _s, 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) + # one reaction per interpolation node: slot s indexes within the point's + # field offset (two P3 edge-interior nodes → slots 0 and 1) + R = np.array([ra[lsec.getFieldOffset(q, 0) + s] for q, s, _c in nodes]) \ + if nodes else np.zeros(0) + flux = _desmear(solver, boundary, xs, R, mass, remove_mean, + edge_node_coords=edge_nodes) 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)) + # vector reaction (traction sigma.n at each node); per-point DOFs are node-major + # with components contiguous per node + Rvec = np.array([ra[lsec.getFieldOffset(q, 0) + s * ncomp: + lsec.getFieldOffset(q, 0) + (s + 1) * ncomp] + for q, s, _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) + Rn = np.array([float(np.dot(nmap[(q, s)], Rvec[i])) + for i, (q, s, _c) in enumerate(nodes)]) if nodes else np.zeros(0) + return xs, _desmear(solver, boundary, xs, Rn, mass, remove_mean, + edge_node_coords=edge_nodes) # 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)] + mass, remove_mean, edge_node_coords=edge_nodes) + for k in range(ncomp)] return xs, (np.column_stack(cols) if nodes else np.zeros((0, ncomp))) diff --git a/tests/parallel/test_1065_boundary_flux_parallel.py b/tests/parallel/test_1065_boundary_flux_parallel.py index ce21e1cea..7f20a9636 100644 --- a/tests/parallel/test_1065_boundary_flux_parallel.py +++ b/tests/parallel/test_1065_boundary_flux_parallel.py @@ -109,6 +109,29 @@ def test_boundary_flux_3d_pointwise_uniform_partition_independent(degree, mass): ) +def test_boundary_flux_degree3_partition_independent(): + """#459 at np >= 2: a degree-3 trace keeps one coordinate per edge-interior node, + and the per-node coordinate build is COLLECTIVE — ranks owning none of the flux + boundary must still participate. T = 1 - y is exact, so every trace node on every + partition reads the exact unit flux.""" + mesh = uw.meshing.StructuredQuadBox( + elementRes=(8, 8), minCoords=(0.0, 0.0), maxCoords=(1.0, 1.0), qdegree=4) + T = uw.discretisation.MeshVariable("T459p", mesh, 1, degree=3) + poisson = uw.systems.Poisson(mesh, u_Field=T) + poisson.constitutive_model = uw.constitutive_models.DiffusionModel + poisson.constitutive_model.Parameters.diffusivity = 1.0 + poisson.f = 0.0 + poisson.add_dirichlet_bc(1.0, "Bottom") + poisson.add_dirichlet_bc(0.0, "Top") + poisson.solve() + for wall, sign in (("Top", -1.0), ("Bottom", +1.0)): + _xs, flux = poisson.boundary_flux(wall) + local_error = float(np.max(np.abs(np.asarray(flux) - sign))) if len(flux) else 0.0 + max_error = uw.mpi.comm.allreduce(local_error, op=MPI.MAX) + assert max_error < 1e-3, ( + f"{wall}: degree-3 flux error {max_error:.3e} at np={uw.mpi.size}") + + if __name__ == "__main__": _b, _r = _flux_diagnostics() if uw.mpi.rank == 0: diff --git a/tests/test_1019_boundary_flux.py b/tests/test_1019_boundary_flux.py index 9e0ec9867..86e67faa1 100644 --- a/tests/test_1019_boundary_flux.py +++ b/tests/test_1019_boundary_flux.py @@ -194,6 +194,88 @@ def test_boundary_flux_p1_trace_2d(mass): ) +def _unit_flux_2d(degree, res=8): + """Unit-box conduction, T = 1 - y: exact at every degree, unit flux on Top/Bottom.""" + mesh = uw.meshing.StructuredQuadBox( + elementRes=(res, res), minCoords=(0.0, 0.0), maxCoords=(1.0, 1.0), qdegree=4) + T = uw.discretisation.MeshVariable(f"T459_{degree}", mesh, 1, degree=degree) + poisson = uw.systems.Poisson(mesh, u_Field=T) + poisson.constitutive_model = uw.constitutive_models.DiffusionModel + poisson.constitutive_model.Parameters.diffusivity = 1.0 + poisson.f = 0.0 + poisson.add_dirichlet_bc(1.0, "Bottom") + poisson.add_dirichlet_bc(0.0, "Top") + poisson.solve() + return poisson + + +@pytest.mark.parametrize("degree", (1, 2, 3)) +def test_boundary_flux_degree_sweep_2d(degree): + """#459: the trace degree must not change the answer. T = 1 - y is exact at every + degree, so each wall node must read the exact unit flux. A degree-3 trace carries + TWO edge-interior reactions per edge; both were keyed by the edge's single + coordinate, collapsed onto one dictionary key, and the recovery silently returned + 0.57-0.74 of the unit flux from 17 of the 25 trace nodes.""" + poisson = _unit_flux_2d(degree) + for wall, sign in (("Top", -1.0), ("Bottom", +1.0)): + xs, flux = poisson.boundary_flux(wall) + assert np.allclose(np.asarray(flux), sign, atol=1e-3), ( + f"{wall} (degree={degree}): flux range " + f"[{np.min(flux)}, {np.max(flux)}], expected {sign}") + if uw.mpi.size == 1: + # every trace node must report — the collapse also DROPPED nodes + assert len(np.asarray(xs)) == degree * 8 + 1 + + +def test_boundary_flux_p3_interior_node_placement(): + """#459: each degree-3 edge-interior reaction pairs with its OWN node coordinate. + T = xy gives dT/dn = x on Top, so the recovered flux is linear in x; swapping the + two interior nodes of an edge (or mis-placing them at the midpoint) displaces the + flux by O(edge/3) ~ 2e-2 at this resolution. Corner columns are excluded: corner + reactions mix both driven walls and the consistent mass spreads that mixture over + about one element (documented corner semantics).""" + mesh = uw.meshing.StructuredQuadBox( + elementRes=(16, 16), minCoords=(0.0, 0.0), maxCoords=(1.0, 1.0), qdegree=4) + T = uw.discretisation.MeshVariable("T459xy", mesh, 1, degree=3) + x, y = mesh.X + poisson = uw.systems.Poisson(mesh, u_Field=T) + poisson.constitutive_model = uw.constitutive_models.DiffusionModel + poisson.constitutive_model.Parameters.diffusivity = 1.0 + poisson.f = 0.0 + for wall in ("Top", "Bottom", "Left", "Right"): + poisson.add_dirichlet_bc(sympy.Matrix([x * y]), wall) + poisson.solve() + + xs, flux = poisson.boundary_flux("Top") # auto -> consistent at degree 3 + xs = np.asarray(xs) + flux = np.asarray(flux) + mid = (xs[:, 0] > 0.35) & (xs[:, 0] < 0.65) + if uw.mpi.size == 1: + assert np.count_nonzero(mid) > 0 + assert np.allclose(flux[mid], xs[mid, 0], atol=5e-3), ( + f"mid-wall flux error {np.max(np.abs(flux[mid] - xs[mid, 0])):.3e} — " + "degree-3 interior nodes mis-placed or mis-ordered") + + +def test_boundary_flux_p3_collapse_guard(monkeypatch): + """#459 negative control: force the pre-fix behaviour (both edge-interior nodes + keyed by the edge's one coordinate) and the de-smear must REFUSE — the silent + overwrite is exactly what returned wrong flux before the fix.""" + from underworld3.utilities import boundary_flux as bf + + poisson = _unit_flux_2d(3, res=4) + true_coords = bf._trace_interior_coords + + def collapsed(solver, degree): + full = true_coords(solver, degree) + return {e: np.repeat(a.mean(axis=0, keepdims=True), len(a), axis=0) + for e, a in full.items()} + + monkeypatch.setattr(bf, "_trace_interior_coords", collapsed) + with pytest.raises(RuntimeError, match="collapse"): + poisson.boundary_flux("Top") + + def test_volume_residual_fields_insert_essential_values(): """#411: compute_volume_residual_fields (Stokes-only diagnostic) missed the #407 insert — its residual must now match _assemble_volume_reaction (the