[finsberg]

Summary

Fixes evaluate_function with broadcast=True (the default) when run on more than one rank:

• NaN values come back as NaN. Previously they became -inf unless rank 0 found the point, so serial and parallel runs gave different results.
• Points outside the mesh come back as NaN instead of -inf.
• Two collective calls in total, instead of one allreduce per point and component.
• Functions whose value size differs from the dofmap block size work, e.g. Raviart–Thomas. Previously these failed with a shape mismatch.

The bug

Each rank set the points it did not find to -inf, then reduced every point with comm.allreduce(value, op=MPI.MAX). The lowercase allreduce reduces Pythonfloats with Python's max, whose result with a NaN depends on the argument order:

>>> max(float("nan"), float("-inf"))
nan
>>> max(float("-inf"), float("nan"))
-inf

A minimal example:

from mpi4py import MPI

import dolfinx
import numpy as np
import scifem

mesh = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, 8, 8)
V = dolfinx.fem.functionspace(mesh, ("Lagrange", 1))
u = dolfinx.fem.Function(V)
u.x.array[:] = np.nan

points = np.array([[0.1, 0.1], [0.9, 0.1], [0.1, 0.9], [0.9, 0.9]])
print(scifem.evaluate_function(u, points).ravel())

Run	Before	After
serial	[nan, nan, nan, nan]	[nan, nan, nan, nan]
mpirun -n 2	[-inf, -inf, nan, nan]	[nan, nan, nan, nan]
mpirun -n 3	[-inf, -inf, nan, -inf]	[nan, nan, nan, nan]

We hit this with activation-time fields, where NaN marks points that never activated.

The fix

• Each point is owned by the lowest rank that finds it, chosen with one integer Allreduce(MIN).
• Only the owner writes the point's value into an array that is zero elsewhere, and one Allreduce(SUM) combines them. Any value, NaN included, therefore passes through unchanged.
• Points that no rank finds are set to NaN.
• The output has prod(u.ufl_shape) columns, instead of the dofmap block size.

Behaviour changes

• Discontinuous functions. A point found on several ranks, e.g. on a process boundary, takes the owning rank's value instead of the maximum over ranks. Allranks still return the same values. The docstring is updated.
• Points outside the mesh are NaN instead of -inf.
• Complex functions return a complex array; previously the values were stored into a float64 array. Real functions still return float64.

Tests

New in tests/test_eval.py:

• NaN is preserved, for scalar and vector P1 functions;
• NaN next to finite values: points where every vertex of the containing cell is NaN are NaN, and points away from the NaN region have the exact value;
• a point outside the mesh is NaN;
• a DG0 function at the mesh vertices, which are shared between cells and ranks: every rank returns the same values;
• Raviart–Thomas (block size 1, value size 2).

Before the fix, the NaN tests fail on 4 ranks, the outside-the-mesh test fails on any number of ranks, and the Raviart–Thomas test fails with a shape mismatch.With the fix, all pass on 1, 2, 3 and 4 ranks.

Verification

Run locally with dolfinx 0.11.0.post0:

• python3 -m pytest: 878 passed, 1164 skipped, 8 xfailed.
• mpiexec -n 4 python3 -m pytest -m "not serial": 880 passed, 1162 skipped, 8 xfailed.
• pre-commit run --all-files: ruff and ruff-format pass. mypy reports 10 errors in src/scifem/solvers.py, which this PR does not touch; the same errors appear on main with this dolfinx version.

🤖 Generated with Claude Code