Skip to content

Commit

Permalink
Fix voting algorithm (#3293)
Browse files Browse the repository at this point in the history 
* Voting algorithm: prioritise owned cells


---------

Co-authored-by: Connor <[email protected]>
  • Loading branch information
@ReubenHill @connorjward
ReubenHill and connorjward authored Oct 30, 2024
1 parent 76d44cf commit ec0329f
Show file tree
Hide file tree
Showing 7 changed files with 230 additions and 34 deletions.
4 changes: 3 additions & 1 deletion firedrake/evaluate.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -52,7 +52,9 @@ extern int locate_cell(struct Function *f,
ref_cell_l1_dist_xtr try_candidate_xtr,
void *temp_ref_coords,
void *found_ref_coords,
double *found_ref_cell_dist_l1);
double *found_ref_cell_dist_l1,
size_t ncells_ignore,
int* cells_ignore);

extern int evaluate(struct Function *f,
double *x,
Expand Down
45 changes: 44 additions & 1 deletion firedrake/locate.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -11,10 +11,13 @@ int locate_cell(struct Function *f,
ref_cell_l1_dist_xtr try_candidate_xtr,
void *temp_ref_coords,
void *found_ref_coords,
double *found_ref_cell_dist_l1)
double *found_ref_cell_dist_l1,
size_t ncells_ignore,
int* cells_ignore)
{
RTError err;
int cell = -1;
int cell_ignore_found = 0;
/* NOTE: temp_ref_coords and found_ref_coords are actually of type
struct ReferenceCoords but can't be declared as such in the function
signature because the dimensions of the reference coordinates in the
Expand Down Expand Up @@ -42,6 +45,16 @@ int locate_cell(struct Function *f,
if (f->extruded == 0) {
for (uint64_t i = 0; i < nids; i++) {
current_ref_cell_dist_l1 = (*try_candidate)(temp_ref_coords, f, ids[i], x);
for (uint64_t j = 0; j < ncells_ignore; j++) {
if (ids[i] == cells_ignore[j]) {
cell_ignore_found = 1;
break;
}
}
if (cell_ignore_found) {
cell_ignore_found = 0;
continue;
}
if (current_ref_cell_dist_l1 <= 0.0) {
/* Found cell! */
cell = ids[i];
Expand All @@ -67,6 +80,16 @@ int locate_cell(struct Function *f,
int c = ids[i] / nlayers;
int l = ids[i] % nlayers;
current_ref_cell_dist_l1 = (*try_candidate_xtr)(temp_ref_coords, f, c, l, x);
for (uint64_t j = 0; j < ncells_ignore; j++) {
if (ids[i] == cells_ignore[j]) {
cell_ignore_found = 1;
break;
}
}
if (cell_ignore_found) {
cell_ignore_found = 0;
continue;
}
if (current_ref_cell_dist_l1 <= 0.0) {
/* Found cell! */
cell = ids[i];
Expand All @@ -91,6 +114,16 @@ int locate_cell(struct Function *f,
if (f->extruded == 0) {
for (int c = 0; c < f->n_cols; c++) {
current_ref_cell_dist_l1 = (*try_candidate)(temp_ref_coords, f, c, x);
for (uint64_t j = 0; j < ncells_ignore; j++) {
if (c == cells_ignore[j]) {
cell_ignore_found = 1;
break;
}
}
if (cell_ignore_found) {
cell_ignore_found = 0;
continue;
}
if (current_ref_cell_dist_l1 <= 0.0) {
/* Found cell! */
cell = c;
Expand All @@ -114,6 +147,16 @@ int locate_cell(struct Function *f,
for (int c = 0; c < f->n_cols; c++) {
for (int l = 0; l < f->n_layers; l++) {
current_ref_cell_dist_l1 = (*try_candidate_xtr)(temp_ref_coords, f, c, l, x);
for (uint64_t j = 0; j < ncells_ignore; j++) {
if (l == cells_ignore[j]) {
cell_ignore_found = 1;
break;
}
}
if (cell_ignore_found) {
cell_ignore_found = 0;
continue;
}
if (current_ref_cell_dist_l1 <= 0.0) {
/* Found cell! */
cell = l;
Expand Down
131 changes: 104 additions & 27 deletions firedrake/mesh.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2494,20 +2494,21 @@ def spatial_index(self):
return spatialindex.from_regions(coords_min, coords_max)

@PETSc.Log.EventDecorator()
def locate_cell(self, x, tolerance=None):
def locate_cell(self, x, tolerance=None, cell_ignore=None):
"""Locate cell containing a given point.
:arg x: point coordinates
:kwarg tolerance: Tolerance for checking if a point is in a cell.
Default is this mesh's :attr:`tolerance` property. Changing
this from default will cause the spatial index to be rebuilt which
can take some time.
:kwarg cell_ignore: Cell number to ignore in the search.
:returns: cell number (int), or None (if the point is not
in the domain)
"""
return self.locate_cell_and_reference_coordinate(x, tolerance=tolerance)[0]
return self.locate_cell_and_reference_coordinate(x, tolerance=tolerance, cell_ignore=cell_ignore)[0]

def locate_reference_coordinate(self, x, tolerance=None):
def locate_reference_coordinate(self, x, tolerance=None, cell_ignore=None):
"""Get reference coordinates of a given point in its cell. Which
cell the point is in can be queried with the locate_cell method.
Expand All @@ -2516,12 +2517,13 @@ def locate_reference_coordinate(self, x, tolerance=None):
Default is this mesh's :attr:`tolerance` property. Changing
this from default will cause the spatial index to be rebuilt which
can take some time.
:kwarg cell_ignore: Cell number to ignore in the search.
:returns: reference coordinates within cell (numpy array) or
None (if the point is not in the domain)
"""
return self.locate_cell_and_reference_coordinate(x, tolerance=tolerance)[1]
return self.locate_cell_and_reference_coordinate(x, tolerance=tolerance, cell_ignore=cell_ignore)[1]

def locate_cell_and_reference_coordinate(self, x, tolerance=None):
def locate_cell_and_reference_coordinate(self, x, tolerance=None, cell_ignore=None):
"""Locate cell containing a given point and the reference
coordinates of the point within the cell.
Expand All @@ -2530,6 +2532,7 @@ def locate_cell_and_reference_coordinate(self, x, tolerance=None):
Default is this mesh's :attr:`tolerance` property. Changing
this from default will cause the spatial index to be rebuilt which
can take some time.
:kwarg cell_ignore: Cell number to ignore in the search.
:returns: tuple either
(cell number, reference coordinates) of type (int, numpy array),
or, when point is not in the domain, (None, None).
Expand All @@ -2538,12 +2541,12 @@ def locate_cell_and_reference_coordinate(self, x, tolerance=None):
if x.size != self.geometric_dimension():
raise ValueError("Point must have the same geometric dimension as the mesh")
x = x.reshape((1, self.geometric_dimension()))
cells, ref_coords, _ = self.locate_cells_ref_coords_and_dists(x, tolerance=tolerance)
cells, ref_coords, _ = self.locate_cells_ref_coords_and_dists(x, tolerance=tolerance, cells_ignore=[[cell_ignore]])
if cells[0] == -1:
return None, None
return cells[0], ref_coords[0]

def locate_cells_ref_coords_and_dists(self, xs, tolerance=None):
def locate_cells_ref_coords_and_dists(self, xs, tolerance=None, cells_ignore=None):
"""Locate cell containing a given point and the reference
coordinates of the point within the cell.
Expand All @@ -2552,6 +2555,11 @@ def locate_cells_ref_coords_and_dists(self, xs, tolerance=None):
Default is this mesh's :attr:`tolerance` property. Changing
this from default will cause the spatial index to be rebuilt which
can take some time.
:kwarg cells_ignore: Cell numbers to ignore in the search for each
point in xs. Shape should be (npoints, n_ignore_pts). Each column
corresponds to a single coordinate in xs. To not ignore any cells,
pass None. To ensure a full cell search for any given point, set
the corresponding entries to -1.
:returns: tuple either
(cell numbers array, reference coordinates array, ref_cell_dists_l1 array)
of type
Expand All @@ -2572,6 +2580,13 @@ def locate_cells_ref_coords_and_dists(self, xs, tolerance=None):
raise ValueError("Point coordinate dimension does not match mesh geometric dimension")
Xs = np.empty_like(xs)
npoints = len(xs)
if cells_ignore is None or cells_ignore[0][0] is None:
cells_ignore = np.full((npoints, 1), -1, dtype=IntType, order="C")
else:
cells_ignore = np.asarray(cells_ignore, dtype=IntType, order="C")
if cells_ignore.shape[0] != npoints:
raise ValueError("Number of cells to ignore does not match number of points")
assert cells_ignore.shape == (npoints, cells_ignore.shape[1])
ref_cell_dists_l1 = np.empty(npoints, dtype=utils.RealType)
cells = np.empty(npoints, dtype=IntType)
assert xs.size == npoints * self.geometric_dimension()
Expand All @@ -2580,7 +2595,9 @@ def locate_cells_ref_coords_and_dists(self, xs, tolerance=None):
Xs.ctypes.data_as(ctypes.POINTER(ctypes.c_double)),
ref_cell_dists_l1.ctypes.data_as(ctypes.POINTER(ctypes.c_double)),
cells.ctypes.data_as(ctypes.POINTER(ctypes.c_int)),
npoints)
npoints,
cells_ignore.shape[1],
cells_ignore)
return cells, Xs, ref_cell_dists_l1

def _c_locator(self, tolerance=None):
Expand All @@ -2596,7 +2613,7 @@ def _c_locator(self, tolerance=None):
IntTypeC = as_cstr(IntType)
src = pq_utils.src_locate_cell(self, tolerance=tolerance)
src += dedent(f"""
int locator(struct Function *f, double *x, double *X, double *ref_cell_dists_l1, {IntTypeC} *cells, {IntTypeC} npoints)
int locator(struct Function *f, double *x, double *X, double *ref_cell_dists_l1, {IntTypeC} *cells, {IntTypeC} npoints, size_t ncells_ignore, int* cells_ignore)
{{
{IntTypeC} j = 0; /* index into x and X */
for({IntTypeC} i=0; i<npoints; i++) {{
Expand All @@ -2609,7 +2626,9 @@ def _c_locator(self, tolerance=None):
/* to_reference_coords and to_reference_coords_xtr are defined in
pointquery_utils.py. If they contain python calls, this loop will
not run at c-loop speed. */
cells[i] = locate_cell(f, &x[j], {self.geometric_dimension()}, &to_reference_coords, &to_reference_coords_xtr, &temp_reference_coords, &found_reference_coords, &ref_cell_dists_l1[i]);
/* cells_ignore has shape (npoints, ncells_ignore) - find the ith row */
int *cells_ignore_i = cells_ignore + i*ncells_ignore;
cells[i] = locate_cell(f, &x[j], {self.geometric_dimension()}, &to_reference_coords, &to_reference_coords_xtr, &temp_reference_coords, &found_reference_coords, &ref_cell_dists_l1[i], ncells_ignore, cells_ignore_i);
for (int k = 0; k < {self.geometric_dimension()}; k++) {{
X[j] = found_reference_coords.X[k];
Expand Down Expand Up @@ -2643,7 +2662,9 @@ def _c_locator(self, tolerance=None):
ctypes.POINTER(ctypes.c_double),
ctypes.POINTER(ctypes.c_double),
ctypes.POINTER(ctypes.c_int),
ctypes.c_size_t]
ctypes.c_size_t,
ctypes.c_size_t,
np.ctypeslib.ndpointer(ctypes.c_int, flags="C_CONTIGUOUS")]
locator.restype = ctypes.c_int
return cache.setdefault(tolerance, locator)

Expand Down Expand Up @@ -4100,12 +4121,18 @@ def _parent_mesh_embedding(
if parent_mesh.geometric_dimension() > parent_mesh.topological_dimension():
# The reference coordinates contain an extra unnecessary dimension
# which we can safely delete
reference_coords = reference_coords[:, :parent_mesh.topological_dimension()]
reference_coords = reference_coords[:, : parent_mesh.topological_dimension()]

locally_visible[:] = parent_cell_nums != -1
ranks[locally_visible] = visible_ranks[parent_cell_nums[locally_visible]]
# see below for why np.inf is used here.
ref_cell_dists_l1[~locally_visible] = np.inf

# ensure that points which a rank thinks it owns are always chosen in a tie
# break by setting the rank to be negative. If multiple ranks think they
# own a point then the one with the highest rank will be chosen.
on_this_rank = ranks == parent_mesh.comm.rank
ranks[on_this_rank] = -parent_mesh.comm.rank
ref_cell_dists_l1_and_ranks = np.stack((ref_cell_dists_l1, ranks), axis=1)

# In parallel there will regularly be disagreements about which cell owns a
Expand All @@ -4124,23 +4151,73 @@ def _parent_mesh_embedding(
ref_cell_dists_l1_and_ranks, op=array_lexicographic_mpi_op
)

# switch ranks back to positive
owned_ref_cell_dists_l1_and_ranks[:, 1] = np.abs(
owned_ref_cell_dists_l1_and_ranks[:, 1]
)
ref_cell_dists_l1_and_ranks[:, 1] = np.abs(ref_cell_dists_l1_and_ranks[:, 1])
ranks = np.abs(ranks)

owned_ref_cell_dists_l1 = owned_ref_cell_dists_l1_and_ranks[:, 0]
owned_ranks = owned_ref_cell_dists_l1_and_ranks[:, 1]

# Any rows where owned_ref_cell_dists_l1_and_ranks and
# ref_cell_dists_l1_and_ranks differ in distance or rank correspond to
# points which are claimed by a cell that we cannot see. We should now
# update our information accordingly. This should only happen for points
# which we've already marked as being owned by a different rank.
extra_missing_points = ~np.all(
owned_ref_cell_dists_l1_and_ranks == ref_cell_dists_l1_and_ranks, axis=1
)
if any(owned_ranks[extra_missing_points] == parent_mesh.comm.rank):
raise RuntimeError(
"Some points have been claimed by a cell that we cannot see, "
"but which we think we own. This should not happen."
)
locally_visible[extra_missing_points] = False
parent_cell_nums[extra_missing_points] = -1
changed_ref_cell_dists_l1 = owned_ref_cell_dists_l1 != ref_cell_dists_l1
changed_ranks = owned_ranks != ranks

# If distance has changed the the point is not in local mesh partition
# since some other cell on another rank is closer.
locally_visible[changed_ref_cell_dists_l1] = False
parent_cell_nums[changed_ref_cell_dists_l1] = -1
# If the rank has changed but the distance hasn't then there was a tie
# break and we need to search for the point again, this time disallowing
# the previously identified cell: if we match the identified owned_rank AND
# the distance is the same then we have found the correct cell. If we
# cannot make a match to owned_rank and distance then we can't see the
# point.
changed_ranks_tied = changed_ranks & ~changed_ref_cell_dists_l1
if any(changed_ranks_tied):
cells_ignore_T = np.asarray([np.copy(parent_cell_nums)])
while any(changed_ranks_tied):
(
parent_cell_nums[changed_ranks_tied],
new_reference_coords,
ref_cell_dists_l1[changed_ranks_tied],
) = parent_mesh.locate_cells_ref_coords_and_dists(
coords_global[changed_ranks_tied],
tolerance,
cells_ignore=cells_ignore_T.T[changed_ranks_tied, :],
)
# delete extra dimension if necessary
if parent_mesh.geometric_dimension() > parent_mesh.topological_dimension():
new_reference_coords = new_reference_coords[:, : parent_mesh.topological_dimension()]
reference_coords[changed_ranks_tied, :] = new_reference_coords
# remove newly lost points
locally_visible[changed_ranks_tied] = (
parent_cell_nums[changed_ranks_tied] != -1
)
changed_ranks_tied &= locally_visible
# if new ref_cell_dists_l1 > owned_ref_cell_dists_l1 then we should
# disregard the point.
locally_visible[changed_ranks_tied] &= (
ref_cell_dists_l1[changed_ranks_tied]
<= owned_ref_cell_dists_l1[changed_ranks_tied]
)
changed_ranks_tied &= locally_visible
# update the identified rank
ranks[changed_ranks_tied] = visible_ranks[
parent_cell_nums[changed_ranks_tied]
]
# if the rank now matches then we have found the correct cell
locally_visible[changed_ranks_tied] &= (
owned_ranks[changed_ranks_tied] == ranks[changed_ranks_tied]
)
# remove these rank matches from changed_ranks_tied
changed_ranks_tied &= ~locally_visible
# add more cells to ignore
cells_ignore_T = np.vstack((
cells_ignore_T,
parent_cell_nums)
)

# Any ranks which are still np.inf are not in the mesh
missing_global_idxs = np.where(owned_ranks == np.inf)[0]
Expand Down
3 changes: 2 additions & 1 deletion firedrake/pointeval_utils.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -140,7 +140,8 @@ def predicate(index):
/* The type definitions and arguments used here are defined as statics in pointquery_utils.py */
double found_ref_cell_dist_l1 = DBL_MAX;
struct ReferenceCoords temp_reference_coords, found_reference_coords;
%(IntType)s cell = locate_cell(f, x, %(geometric_dimension)d, &to_reference_coords, &to_reference_coords_xtr, &temp_reference_coords, &found_reference_coords, &found_ref_cell_dist_l1);
int cells_ignore[1] = {-1};
%(IntType)s cell = locate_cell(f, x, %(geometric_dimension)d, &to_reference_coords, &to_reference_coords_xtr, &temp_reference_coords, &found_reference_coords, &found_ref_cell_dist_l1, 1, cells_ignore);
if (cell == -1) {
return -1;
}
Expand Down
35 changes: 35 additions & 0 deletions tests/regression/test_interpolate_cross_mesh.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -756,6 +756,41 @@ def test_exact_refinement_parallel():
test_exact_refinement()


def voting_algorithm_edgecases(nprocs):
# this triggers lots of cases where the VOM voting algorithm has to deal
# with points being claimed by multiple ranks: there are cases where each
# rank will claim another one owns a point, for example, and yet also all
# claim zero distance to the reference cell!
s = nprocs
nx = 2 * s
mx = 3 * nx
mh = [UnitCubeMesh(nx, nx, nx),
UnitCubeMesh(mx, mx, mx)]
family = "Lagrange"
degree = 1
Vc = FunctionSpace(mh[0], family, degree=degree)
Vf = FunctionSpace(mh[1], family, degree=degree)
uc = Function(Vc).interpolate(SpatialCoordinate(mh[0])[0])
uf = Function(Vf).interpolate(uc)
uf2 = Function(Vf).interpolate(SpatialCoordinate(mh[1])[0])
assert np.isclose(errornorm(uf, uf2), 0.0)


@pytest.mark.parallel(nprocs=2)
def test_voting_algorithm_edgecases_2_ranks():
voting_algorithm_edgecases(2)


@pytest.mark.parallel(nprocs=3)
def test_voting_algorithm_edgecases_3_ranks():
voting_algorithm_edgecases(3)


@pytest.mark.parallel(nprocs=4)
def test_voting_algorithm_edgecases_4_ranks():
voting_algorithm_edgecases(4)


@pytest.mark.parallel
@pytest.mark.parametrize('periodic', [False, True])
def test_interpolate_cross_mesh_interval(periodic):
Expand Down
Loading

0 comments on commit ec0329f

Please sign in to comment.