Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Fix point location in parallel #3804

Closed
wants to merge 7 commits into from
Closed

Fix point location in parallel #3804

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
7 commits
Select commit Hold shift + click to select a range
33dbc3d
having a go
connorjward Oct 18, 2024
4dafd10
fixup
connorjward Oct 18, 2024
0f4d524
cleanup
connorjward Oct 18, 2024
b7c0919
Things appear to work
connorjward Oct 18, 2024
2157267
add test
connorjward Oct 18, 2024
806b9be
cleanup test
connorjward Oct 18, 2024
b0a87cf
commenting
connorjward Oct 18, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 2 additions & 0 deletions firedrake/evaluate.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -48,6 +48,7 @@ typedef PetscReal (*ref_cell_l1_dist_xtr)(void *data_,
extern int locate_cell(struct Function *f,
double *x,
int dim,
int num_owned_cells,
ref_cell_l1_dist try_candidate,
ref_cell_l1_dist_xtr try_candidate_xtr,
void *temp_ref_coords,
Expand All @@ -56,6 +57,7 @@ extern int locate_cell(struct Function *f,

extern int evaluate(struct Function *f,
double *x,
int num_owned_cells,
PetscScalar *result);

#ifdef __cplusplus
Expand Down
11 changes: 10 additions & 1 deletion firedrake/function.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -576,7 +576,7 @@ def _c_evaluate(self, tolerance=None):
return cache[tolerance]
except KeyError:
result = make_c_evaluate(self, tolerance=tolerance)
result.argtypes = [POINTER(_CFunction), POINTER(c_double), c_void_p]
result.argtypes = [POINTER(_CFunction), POINTER(c_double), c_int, c_void_p]
result.restype = c_int
return cache.setdefault(tolerance, result)

Expand Down Expand Up @@ -650,8 +650,17 @@ def at(self, arg, *args, **kwargs):

def single_eval(x, buf):
r"""Helper function to evaluate at a single point."""
mesh = self.ufl_domain()
if mesh.extruded:
if mesh.variable_layers:
raise NotImplementedError("Current codegen for extruded meshes "
"assumes constant layers")
num_owned_cells = mesh.cell_set.size * mesh.layers
else:
num_owned_cells = mesh.cell_set.size
err = self._c_evaluate(tolerance=tolerance)(self._ctypes,
x.ctypes.data_as(POINTER(c_double)),
num_owned_cells,
buf.ctypes.data_as(c_void_p))
if err == -1:
raise PointNotInDomainError(self.function_space().mesh(), x.reshape(-1))
Expand Down
224 changes: 123 additions & 101 deletions firedrake/locate.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -4,140 +4,162 @@
#include <float.h>
#include <evaluate.h>

/**
* Check whether a point is either within, or sufficiently close to the current cell.
*
* To avoid the case where different ranks each claim the other to own the point,
* points are preferentially claimed by owned cells.
*/
bool check_cell(int64_t current_cell,
double current_ref_cell_dist_l1,
int num_owned_cells,
int64_t *found_cell,
bool *found_cell_is_owned,
void *found_ref_coords,
void *temp_ref_coords,
double *found_ref_cell_dist_l1)
{
if (current_cell < num_owned_cells) {
if (current_ref_cell_dist_l1 <= 0) {
// found the cell
*found_cell = current_cell;
memcpy(found_ref_coords, temp_ref_coords, sizeof(struct ReferenceCoords));
*found_ref_cell_dist_l1 = current_ref_cell_dist_l1;
return true;
}
else if (current_ref_cell_dist_l1 < tolerance &&
current_ref_cell_dist_l1 < *found_ref_cell_dist_l1)
{
/* Close to cell within tolerance so could be this cell */
*found_cell_is_owned = true;

*found_cell = current_cell;
memcpy(found_ref_coords, temp_ref_coords, sizeof(struct ReferenceCoords));
*found_ref_cell_dist_l1 = current_ref_cell_dist_l1;
}
}
else {
// Only consider ghost cells if an owned one is yet to be found
if (!(*found_cell_is_owned) &&
current_ref_cell_dist_l1 < tolerance &&
current_ref_cell_dist_l1 < *found_ref_cell_dist_l1)
{
*found_cell = current_cell;
memcpy(found_ref_coords, temp_ref_coords, sizeof(struct ReferenceCoords));
*found_ref_cell_dist_l1 = current_ref_cell_dist_l1;
}
}
return false;
}

int locate_cell(struct Function *f,
double *x,
int dim,
ref_cell_l1_dist try_candidate,
ref_cell_l1_dist_xtr try_candidate_xtr,
void *temp_ref_coords,
void *found_ref_coords,
double *found_ref_cell_dist_l1)
double *x,
int dim,
int num_owned_cells,
ref_cell_l1_dist try_candidate,
ref_cell_l1_dist_xtr try_candidate_xtr,
void *temp_ref_coords,
void *found_ref_coords,
double *found_ref_cell_dist_l1)
{
RTError err;
int cell = -1;
/* 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
ReferenceCoords struct are defined by python when the code which
surrounds this is declared in pointquery_utils.py. We cast when we use the
ref_coords_copy function and trust that the underlying memory which the
pointers refer to is updated as necessary. */
double ref_cell_dist_l1 = DBL_MAX;
double current_ref_cell_dist_l1 = -0.5;
struct ReferenceCoords but can't be declared as such in the function
signature because the dimensions of the reference coordinates in the
ReferenceCoords struct are defined by python when the code which
surrounds this is declared in pointquery_utils.py. We cast when we use the
ref_coords_copy function and trust that the underlying memory which the
pointers refer to is updated as necessary. */
/* NOTE: `tolerance`, which is used throughout this funciton, is a static
variable defined outside this function when putting together all the C
code that needs to be compiled - see pointquery_utils.py */

*found_ref_cell_dist_l1 = DBL_MAX;

RTError err;
int64_t found_cell = -1;
bool found_cell_is_owned = false;

if (f->sidx) {
/* We treat our list of candidate cells (ids) from libspatialindex's
Index_Intersects_id as our source of truth: the point must be in
one of the cells. */
int64_t *ids = NULL;
uint64_t nids = 0;
/* We treat our list of candidate cells (ids) from libspatialindex's
Index_Intersects_id as our source of truth: the point must be in
one of the cells. */
err = Index_Intersects_id(f->sidx, x, x, dim, &ids, &nids);
if (err != RT_None) {
fputs("ERROR: Index_Intersects_id failed in libspatialindex!", stderr);
return -1;
}
if (f->extruded == 0) {

if (!f->extruded) {
for (uint64_t i = 0; i < nids; i++) {
current_ref_cell_dist_l1 = (*try_candidate)(temp_ref_coords, f, ids[i], x);
if (current_ref_cell_dist_l1 <= 0.0) {
/* Found cell! */
cell = ids[i];
memcpy(found_ref_coords, temp_ref_coords, sizeof(struct ReferenceCoords));
found_ref_cell_dist_l1[0] = current_ref_cell_dist_l1;
break;
}
else if (current_ref_cell_dist_l1 < ref_cell_dist_l1) {
/* getting closer... */
ref_cell_dist_l1 = current_ref_cell_dist_l1;
if (ref_cell_dist_l1 < tolerance) {
/* Close to cell within tolerance so could be this cell */
cell = ids[i];
memcpy(found_ref_coords, temp_ref_coords, sizeof(struct ReferenceCoords));
found_ref_cell_dist_l1[0] = ref_cell_dist_l1;
}
}
int64_t current_cell = ids[i];
double current_ref_cell_dist_l1 = (*try_candidate)(temp_ref_coords, f, current_cell, x);
bool found = check_cell(current_cell,
current_ref_cell_dist_l1,
num_owned_cells,
&found_cell,
&found_cell_is_owned,
found_ref_coords,
temp_ref_coords,
found_ref_cell_dist_l1);
if (found) break;
}
}
else {
for (uint64_t i = 0; i < nids; i++) {
int64_t current_cell = ids[i];

int nlayers = f->n_layers;
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);
if (current_ref_cell_dist_l1 <= 0.0) {
/* Found cell! */
cell = ids[i];
memcpy(found_ref_coords, temp_ref_coords, sizeof(struct ReferenceCoords));
found_ref_cell_dist_l1[0] = current_ref_cell_dist_l1;
break;
}
else if (current_ref_cell_dist_l1 < ref_cell_dist_l1) {
/* getting closer... */
ref_cell_dist_l1 = current_ref_cell_dist_l1;
if (ref_cell_dist_l1 < tolerance) {
/* Close to cell within tolerance so could be this cell */
cell = ids[i];
memcpy(found_ref_coords, temp_ref_coords, sizeof(struct ReferenceCoords));
found_ref_cell_dist_l1[0] = ref_cell_dist_l1;
}
}
int c = current_cell / nlayers;
int l = current_cell % nlayers;
double current_ref_cell_dist_l1 = (*try_candidate_xtr)(temp_ref_coords, f, c, l, x);

bool found = check_cell(current_cell,
current_ref_cell_dist_l1,
num_owned_cells,
&found_cell,
&found_cell_is_owned,
found_ref_coords,
temp_ref_coords,
found_ref_cell_dist_l1);
if (found) break;
}
}
free(ids);
} else {
if (f->extruded == 0) {
if (!f->extruded) {
for (int c = 0; c < f->n_cols; c++) {
current_ref_cell_dist_l1 = (*try_candidate)(temp_ref_coords, f, c, x);
if (current_ref_cell_dist_l1 <= 0.0) {
/* Found cell! */
cell = c;
memcpy(found_ref_coords, temp_ref_coords, sizeof(struct ReferenceCoords));
found_ref_cell_dist_l1[0] = current_ref_cell_dist_l1;
break;
}
else if (current_ref_cell_dist_l1 < ref_cell_dist_l1) {
/* getting closer... */
ref_cell_dist_l1 = current_ref_cell_dist_l1;
if (ref_cell_dist_l1 < tolerance) {
/* Close to cell within tolerance so could be this cell */
cell = c;
memcpy(found_ref_coords, temp_ref_coords, sizeof(struct ReferenceCoords));
found_ref_cell_dist_l1[0] = ref_cell_dist_l1;
}
}
double current_ref_cell_dist_l1 = (*try_candidate)(temp_ref_coords, f, c, x);
bool found = check_cell(c,
current_ref_cell_dist_l1,
num_owned_cells,
&found_cell,
&found_cell_is_owned,
found_ref_coords,
temp_ref_coords,
found_ref_cell_dist_l1);
if (found) break;
}
}
else {
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);
if (current_ref_cell_dist_l1 <= 0.0) {
/* Found cell! */
cell = l;
memcpy(found_ref_coords, temp_ref_coords, sizeof(struct ReferenceCoords));
found_ref_cell_dist_l1[0] = current_ref_cell_dist_l1;
break;
}
else if (current_ref_cell_dist_l1 < ref_cell_dist_l1) {
/* getting closer... */
ref_cell_dist_l1 = current_ref_cell_dist_l1;
if (ref_cell_dist_l1 < tolerance) {
/* Close to cell within tolerance so could be this cell */
cell = l;
memcpy(found_ref_coords, temp_ref_coords, sizeof(struct ReferenceCoords));
found_ref_cell_dist_l1[0] = ref_cell_dist_l1;
}
}
}
if (cell != -1) {
cell = c * f->n_layers + cell;
break;
int64_t current_cell = c * f->n_layers + l;
double current_ref_cell_dist_l1 = (*try_candidate_xtr)(temp_ref_coords, f, c, l, x);
bool found = check_cell(current_cell,
current_ref_cell_dist_l1,
num_owned_cells,
&found_cell,
&found_cell_is_owned,
found_ref_coords,
temp_ref_coords,
found_ref_cell_dist_l1);
if (found) break;
}
if (found_cell != -1 && found_cell_is_owned) break;
}
}
}
return cell;
return found_cell;
}
19 changes: 15 additions & 4 deletions firedrake/mesh.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -2575,12 +2575,22 @@ def locate_cells_ref_coords_and_dists(self, xs, tolerance=None):
ref_cell_dists_l1 = np.empty(npoints, dtype=utils.RealType)
cells = np.empty(npoints, dtype=IntType)
assert xs.size == npoints * self.geometric_dimension()

if self.extruded:
if self.variable_layers:
raise NotImplementedError("Current codegen for extruded meshes "
"assumes constant layers")
num_owned_cells = self.cell_set.size * self.layers
else:
num_owned_cells = self.cell_set.size

self._c_locator(tolerance=tolerance)(self.coordinates._ctypes,
xs.ctypes.data_as(ctypes.POINTER(ctypes.c_double)),
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,
num_owned_cells)
return cells, Xs, ref_cell_dists_l1

def _c_locator(self, tolerance=None):
Expand All @@ -2596,7 +2606,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, int num_owned_cells)
{{
{IntTypeC} j = 0; /* index into x and X */
for({IntTypeC} i=0; i<npoints; i++) {{
Expand All @@ -2609,7 +2619,7 @@ 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[i] = locate_cell(f, &x[j], {self.geometric_dimension()}, num_owned_cells, &to_reference_coords, &to_reference_coords_xtr, &temp_reference_coords, &found_reference_coords, &ref_cell_dists_l1[i]);

for (int k = 0; k < {self.geometric_dimension()}; k++) {{
X[j] = found_reference_coords.X[k];
Expand Down Expand Up @@ -2643,7 +2653,8 @@ 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_int]
locator.restype = ctypes.c_int
return cache.setdefault(tolerance, locator)

Expand Down
4 changes: 2 additions & 2 deletions firedrake/pointeval_utils.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -135,12 +135,12 @@ def predicate(index):
%(scalar_type)s const *__restrict__ coords, %(scalar_type)s const *__restrict__ f, %(wrapper_map_args)s);


int evaluate(struct Function *f, double *x, %(scalar_type)s *result)
int evaluate(struct Function *f, double *x, int num_owned_cells, %(scalar_type)s *result)
{
/* 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);
%(IntType)s cell = locate_cell(f, x, %(geometric_dimension)d, num_owned_cells, &to_reference_coords, &to_reference_coords_xtr, &temp_reference_coords, &found_reference_coords, &found_ref_cell_dist_l1);
if (cell == -1) {
return -1;
}
Expand Down
Loading
Loading