Duplicates can happen when evaluating function at points - safe version

General
1

Hello all,

I have found an interesting detail when evaluating function at certain points as shown in this tutorial by @dokken.

I think this is what is happening (or at least this is one trigger of this):
If a point lays on a facet of two cells the collision mechanism finds two links. When these two cells are own by different processes the point is evaluated twice.

This line cells.append(colliding_cells.links(i)[0]) returns only one cell (the first it found) if the two cells are found on the same process, but cannot handle the case where the two cells are own by different processes.

This scenario is not likely but it happens.
I was getting bunch of unexpected behavior before I have found out.

Consider this MWE:

from dolfinx.mesh import CellType, create_rectangle
from mpi4py import MPI
from dolfinx import fem
from dolfinx import geometry
import numpy as np

mesh = create_rectangle(MPI.COMM_WORLD, [[0, 0], [1, 1]], [64, 64], CellType.triangle)
V = fem.FunctionSpace(mesh, ("CG", 1))
u = fem.Function(V)

def probe_function(p_coords, f, make_safe=False):
    domain = f.function_space.mesh
    points = np.array(p_coords)
    bb_tree = geometry.bb_tree(domain, domain.topology.dim)
    cell_candidates = geometry.compute_collisions_points(bb_tree, points)
    colliding_cells = geometry.compute_colliding_cells(domain, cell_candidates, points)
    cells = []
    points_on_proc = []
    idx = []
    for i, point in enumerate(points):
        if len(colliding_cells.links(i)) > 0:
            points_on_proc.append(point)
            cells.append(colliding_cells.links(i)[0])
            idx.append(i)
    points_on_proc = np.array(points_on_proc, dtype=np.float64)
    idx = np.array(idx, dtype=int)
    fv = f.eval(points_on_proc, cells)
    fv = domain.comm.allgather(fv)
    idx = domain.comm.allgather(idx)
    fv = np.vstack(fv)
    idx = np.concatenate(idx)
    if make_safe:
        sorted_idx_arg = np.argsort(idx) # sort ascending (might contain duplicates)
        fv = fv[sorted_idx_arg]    # reorder values to original order of the p_coords
        idx = idx[sorted_idx_arg]  # indexes to original order (duplicates still in)
        _, u_idx = np.unique(idx, return_index=True) # create idexing for unique values
        fv = fv[u_idx] # remove duplicates from values
    return fv.flatten()

points = [[0.0, 0.0, 0.0], [0.5, 0.5, 0.0]]
vals_unsafe = probe_function(points, u)
vals_safe = probe_function(points, u, make_safe=True)
if MPI.COMM_WORLD.rank == 0:
    print(vals_unsafe)
    print(vals_safe)

with mpirun -n 1

[0. 0.]
[0. 0.]

with mpirun -n 2

[0. 0. 0.]  # three values for two points!!
[0. 0.]

Also the “safe” version respects the order in points.
Perhaps this should be reflected in the tutorial @dokken ?

1 Like
2

There exists functionality to determine a unique owner, but it is way more costly.
You can use:

which is documented at:

It takes in:
A unique set of points per process, thus in the application below, points should only be sent in on rank 0.
It returns:
What points the current process now owns and where they should be sent to (if all points come from rank 0, they naturally have to go back to rank 0).

Using allgather, as done above, is not scalable on large systems.

2 Likes
3

Awesome, thank you for the reference. I will try out.