diff --git a/src/underworld3/function/functions_unit_system.py b/src/underworld3/function/functions_unit_system.py index ecfea55b..a92926bf 100644 --- a/src/underworld3/function/functions_unit_system.py +++ b/src/underworld3/function/functions_unit_system.py @@ -25,6 +25,8 @@ - Single source of truth for unit handling logic """ +import warnings + import numpy as np import underworld3 as uw @@ -293,17 +295,20 @@ def _evaluate_impl( evalf=evalf, rbf=rbf_flag, data_layout=data_layout, - check_extrapolated=check_extrapolated, + check_extrapolated=True, # always: free, and needed for the guard below force_l2=force_l2_flag, smoothing=smoothing, ) - # Step 4: Unpack extrapolation flag if needed - if check_extrapolated: - raw_values, extrapolated = raw_result_nondim - else: - raw_values = raw_result_nondim - extrapolated = None + # Step 4: Unpack the extrapolation flag. + # + # It is requested UNCONDITIONALLY. Measured cost of asking for it: none — + # 0.0073 s/call against 0.0075 s/call without, on 1969 points, i.e. inside the + # noise. The locator has to decide whether it located a point in order to fall + # back, so the mask is already known and returning it is free. + raw_values, extrapolated = raw_result_nondim + + _warn_if_points_are_not_owned(extrapolated) # Step 5: Re-dimensionalize and wrap with units # GATEWAY PRINCIPLE: evaluate() ALWAYS returns dimensional values when units are known @@ -763,6 +768,50 @@ def _apply_monotone_limit( @uw.timing.routine_timer_decorator +def _warn_if_points_are_not_owned(extrapolated): + """Warn when a parallel `evaluate` silently answered for points this rank does + not own. + + `evaluate` is RANK-LOCAL: it answers from this rank's portion of the mesh. In + serial that is the whole mesh and everything is exact. In parallel, a query + point owned by another rank is not located here, so the locator falls back to + extrapolation and returns a plausible number that is simply wrong — and two + ranks asked the identical question return different answers. + + Measured (#606), P2 field holding x^2+2y^2 on a unit box, querying every rank's + own DOF coordinates allgathered — so every point is a mesh node: + + np=1 all coords max error 8.9e-16 + np=2 own coords max error 6.7e-16 <- rank-local use is exact + np=2 all coords max error 1.48 <- on a field whose range is 3 + np=4 all coords max error 2.59, two thirds of points wrong + + The extrapolation mask is an EXACT detector of those points: at np=2, 69 of 166 + flagged and 69 wrong, with no wrong-but-unflagged and no flagged-but-fine; at + np=4, 120 and 120. Maximum error among unflagged points was 8.9e-16. So this + warning has neither false negatives nor false positives on the case that + motivated it. + + Serial is deliberately NOT warned about: there, an extrapolated point is one + genuinely outside the domain, which is a legitimate thing to ask for. + """ + if uw.mpi.size == 1 or extrapolated is None: + return + flagged = int(np.count_nonzero(extrapolated)) + if not flagged: + return + total = int(np.asarray(extrapolated).size) + warnings.warn( + f"evaluate() is rank-local: {flagged} of {total} query points are not " + f"located on this rank (rank {uw.mpi.rank} of {uw.mpi.size}), so their " + "values were EXTRAPOLATED and are wrong — different ranks will disagree " + "for the same query. Use global_evaluate() for points that may be owned " + "elsewhere, or restrict the query to this rank's own coordinates (#606).", + RuntimeWarning, + stacklevel=3, + ) + + def evaluate( expr, coords, @@ -787,6 +836,21 @@ def evaluate( automatically. With the default ``monotone=False`` the result is bit-identical to the historical ``evaluate``. + .. warning:: + **This is RANK-LOCAL.** It answers from this rank's portion of the + mesh. In serial that is the whole mesh. In parallel, a point owned by + another rank cannot be located here, so it is EXTRAPOLATED and the + value returned is wrong — and two ranks asked the same question return + different answers. Measured on a P2 field exact in P2 (#606): querying + every rank's own DOF coordinates allgathered gave a maximum error of + 1.48 at ``np=2`` and 2.59 at ``np=4``, on a field whose whole range is + 3, with two thirds of the points wrong at ``np=4`` — while the same + query restricted to each rank's own coordinates was exact to 1e-16. + + Use :func:`global_evaluate` for points that may be owned elsewhere. A + parallel call that extrapolates any point now emits a + ``RuntimeWarning`` naming the count. + Parameters ---------- expr : sympy expression or UWexpression diff --git a/tests/parallel/test_1075_evaluate_rank_locality.py b/tests/parallel/test_1075_evaluate_rank_locality.py new file mode 100644 index 00000000..4850fb58 --- /dev/null +++ b/tests/parallel/test_1075_evaluate_rank_locality.py @@ -0,0 +1,79 @@ +#!/usr/bin/env python3 +"""`evaluate` is rank-local, and must say so when that costs the caller (#606). + +`uw.function.evaluate` answers from the calling rank's portion of the mesh. In +serial that is the whole mesh. In parallel, a query point owned by another rank +cannot be located here, so the locator extrapolates and returns a plausible +number that is simply wrong — and two ranks asked the identical question return +different answers. + +Measured on a P2 field holding `x^2 + 2 y^2`, which P2 represents exactly, so any +discrepancy is location rather than approximation. Querying every rank's own DOF +coordinates allgathered — every point a mesh node — gave a maximum error of 1.48 +at np=2 and 2.59 at np=4 on a field whose range is [0, 3], while the same query +restricted to each rank's own coordinates was exact to 1e-16. + +The extrapolation mask is an exact detector of those points (69 flagged / 69 +wrong at np=2, 120 / 120 at np=4, with no wrong-but-unflagged and no +flagged-but-fine), which is what makes warning on it sound. +""" +import warnings + +import numpy as np +import pytest + +import underworld3 as uw + +pytestmark = [pytest.mark.level_1, pytest.mark.tier_a] + + +def _field(): + mesh = uw.meshing.UnstructuredSimplexBox( + minCoords=(0.0, 0.0), maxCoords=(1.0, 1.0), cellSize=0.2, qdegree=3) + field = uw.discretisation.MeshVariable("f1075", mesh, 1, degree=2) + coords = np.asarray(field.coords, dtype=float) + field.data[:, 0] = coords[:, 0] ** 2 + 2.0 * coords[:, 1] ** 2 + return field, coords + + +def _evaluate(field, points): + with warnings.catch_warnings(record=True) as caught: + warnings.simplefilter("always") + values = np.asarray( + uw.function.evaluate(field.sym[0], points), dtype=float).ravel() + warned = [w for w in caught if "rank-local" in str(w.message)] + truth = points[:, 0] ** 2 + 2.0 * points[:, 1] ** 2 + return float(np.abs(values - truth).max()), bool(warned) + + +def test_own_coordinates_are_exact_and_silent(): + """The rank-local contract honoured: no warning, and no error to warn about. + + This is the negative control. A warning that fired here would be noise on + every correct parallel use of `evaluate`. + """ + field, coords = _field() + error, warned = _evaluate(field, coords) + assert error < 1.0e-12, f"own coordinates should be exact, got {error:.2e}" + assert not warned, "warned about a query that was entirely rank-local" + + +def test_unowned_coordinates_warn(): + """A query spanning ranks must not answer silently. + + In serial there is nothing to warn about — every point is owned — so the + assertion is conditioned on rank count rather than skipped, which keeps the + exactness check alive at np=1. + """ + field, coords = _field() + every = np.vstack([g for g in uw.mpi.comm.allgather(coords) if len(g)]) + error, warned = _evaluate(field, every) + + if uw.mpi.size == 1: + assert error < 1.0e-12, f"serial must be exact, got {error:.2e}" + assert not warned, "warned in serial, where every point is owned" + return + + assert warned, ( + f"evaluate returned a max error of {error:.2e} for points this rank does " + "not own, and said nothing")