diff --git a/firedrake/evaluate.h b/firedrake/evaluate.h index 738a573867..39005bce62 100644 --- a/firedrake/evaluate.h +++ b/firedrake/evaluate.h @@ -9,13 +9,13 @@ extern "C" { struct Function { /* Number of cells in the base mesh */ - int n_cols; + PetscInt n_cols; /* 1 if extruded, 0 if not */ int extruded; /* number of layers for extruded, otherwise 1 */ - int n_layers; + PetscInt n_layers; /* Coordinate values and node mapping */ PetscScalar *coords; @@ -36,25 +36,25 @@ struct Function { typedef PetscReal (*ref_cell_l1_dist)(void *data_, struct Function *f, - int cell, + PetscInt cell, double *x); typedef PetscReal (*ref_cell_l1_dist_xtr)(void *data_, struct Function *f, - int cell, - int layer, + PetscInt cell, + PetscInt layer, double *x); -extern int locate_cell(struct Function *f, +extern PetscInt 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, + PetscReal *found_ref_cell_dist_l1, size_t ncells_ignore, - int* cells_ignore); + PetscInt* cells_ignore); extern int evaluate(struct Function *f, double *x, diff --git a/firedrake/function.py b/firedrake/function.py index b447a83174..1a21412939 100644 --- a/firedrake/function.py +++ b/firedrake/function.py @@ -37,9 +37,9 @@ class _CFunction(ctypes.Structure): r"""C struct collecting data from a :class:`Function`""" - _fields_ = [("n_cols", c_int), + _fields_ = [("n_cols", as_ctypes(IntType)), ("extruded", c_int), - ("n_layers", c_int), + ("n_layers", as_ctypes(IntType)), ("coords", c_void_p), ("coords_map", POINTER(as_ctypes(IntType))), ("f", c_void_p), diff --git a/firedrake/interpolation.py b/firedrake/interpolation.py index bf89a0c114..6a6a4fdedc 100644 --- a/firedrake/interpolation.py +++ b/firedrake/interpolation.py @@ -1552,8 +1552,9 @@ def _create_permutation_mat(self, mat_type: Literal["aij", "baij"]) -> PETSc.Mat end = start + self.source_size[0] contiguous_indices = numpy.arange(start, end, dtype=IntType) perm = numpy.zeros(self.nleaves, dtype=IntType) # result stored in here - self.sf.bcastBegin(MPI.INT, contiguous_indices, perm, MPI.REPLACE) - self.sf.bcastEnd(MPI.INT, contiguous_indices, perm, MPI.REPLACE) + mpi_int = MPI._typedict[numpy.dtype(IntType).char] + self.sf.bcastBegin(mpi_int, contiguous_indices, perm, MPI.REPLACE) + self.sf.bcastEnd(mpi_int, contiguous_indices, perm, MPI.REPLACE) rows = numpy.arange(self.target_size[0] + 1, dtype=IntType) # Vector and Tensor valued functions are stored in a flattened array, so # we need to space out the column indices according to the block size diff --git a/firedrake/locate.c b/firedrake/locate.c index a243e2119f..53350cb533 100644 --- a/firedrake/locate.c +++ b/firedrake/locate.c @@ -4,19 +4,19 @@ #include #include -int locate_cell(struct Function *f, +PetscInt 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, + PetscReal *found_ref_cell_dist_l1, size_t ncells_ignore, - int* cells_ignore) + PetscInt* cells_ignore) { RTError err; - int cell = -1; + PetscInt 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 @@ -25,8 +25,8 @@ int locate_cell(struct Function *f, 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; + PetscReal ref_cell_dist_l1 = PETSC_MAX_REAL; + PetscReal current_ref_cell_dist_l1 = -0.5; /* 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 */ @@ -45,7 +45,7 @@ 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++) { + for (size_t j = 0; j < ncells_ignore; j++) { if (ids[i] == cells_ignore[j]) { cell_ignore_found = 1; break; @@ -76,11 +76,11 @@ int locate_cell(struct Function *f, } else { for (uint64_t i = 0; i < nids; i++) { - int nlayers = f->n_layers; - int c = ids[i] / nlayers; - int l = ids[i] % nlayers; + PetscInt nlayers = f->n_layers; + PetscInt c = ids[i] / nlayers; + PetscInt 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++) { + for (size_t j = 0; j < ncells_ignore; j++) { if (ids[i] == cells_ignore[j]) { cell_ignore_found = 1; break; @@ -112,9 +112,9 @@ int locate_cell(struct Function *f, free(ids); } else { if (f->extruded == 0) { - for (int c = 0; c < f->n_cols; c++) { + for (PetscInt 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++) { + for (size_t j = 0; j < ncells_ignore; j++) { if (c == cells_ignore[j]) { cell_ignore_found = 1; break; @@ -144,10 +144,10 @@ int locate_cell(struct Function *f, } } else { - for (int c = 0; c < f->n_cols; c++) { - for (int l = 0; l < f->n_layers; l++) { + for (PetscInt c = 0; c < f->n_cols; c++) { + for (PetscInt 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++) { + for (size_t j = 0; j < ncells_ignore; j++) { if (l == cells_ignore[j]) { cell_ignore_found = 1; break; diff --git a/firedrake/mesh.py b/firedrake/mesh.py index 25ad2db025..9eb26aecc5 100644 --- a/firedrake/mesh.py +++ b/firedrake/mesh.py @@ -37,7 +37,7 @@ import firedrake.extrusion_utils as eutils import firedrake.cython.spatialindex as spatialindex import firedrake.utils as utils -from firedrake.utils import as_cstr, IntType, RealType +from firedrake.utils import IntType, IntType_c, RealType, RealType_c, as_ctypes from firedrake.logging import info_red, logger from firedrake.parameters import parameters from firedrake.petsc import PETSc, DEFAULT_PARTITIONER @@ -2706,11 +2706,12 @@ def locate_cells_ref_coords_and_dists(self, xs, tolerance=None, cells_ignore=Non tolerance = self.tolerance else: self.tolerance = tolerance - xs = np.asarray(xs, dtype=utils.ScalarType) - xs = xs.real.copy() + # `xs` are the physical coordinates we query the rtree with. + # libspatialindex requires these to be of type double + xs = np.asarray(xs).real.astype(np.float64, order="C") if xs.shape[1] != self.geometric_dimension: raise ValueError("Point coordinate dimension does not match mesh geometric dimension") - Xs = np.empty_like(xs) + Xs = np.empty_like(xs, dtype=RealType) 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") @@ -2719,14 +2720,14 @@ def locate_cells_ref_coords_and_dists(self, xs, tolerance=None, cells_ignore=Non 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) + ref_cell_dists_l1 = np.empty(npoints, dtype=RealType) cells = np.empty(npoints, dtype=IntType) assert xs.size == npoints * self.geometric_dimension run_c = self._c_locator(tolerance=tolerance) - cells_data = cells.ctypes.data_as(ctypes.POINTER(ctypes.c_int)) - ref_cells_dists = ref_cell_dists_l1.ctypes.data_as(ctypes.POINTER(ctypes.c_double)) + cells_data = cells.ctypes.data_as(ctypes.POINTER(as_ctypes(IntType))) + ref_cells_dists = ref_cell_dists_l1.ctypes.data_as(ctypes.POINTER(as_ctypes(RealType))) xs_data = xs.ctypes.data_as(ctypes.POINTER(ctypes.c_double)) - Xs_data = Xs.ctypes.data_as(ctypes.POINTER(ctypes.c_double)) + Xs_data = Xs.ctypes.data_as(ctypes.POINTER(as_ctypes(RealType))) with PETSc.Log.Event("c_locator_run"): run_c(self.coordinates._ctypes, xs_data, Xs_data, ref_cells_dists, cells_data, npoints, cells_ignore.shape[1], cells_ignore) return cells, Xs, ref_cell_dists_l1 @@ -2741,13 +2742,12 @@ def _c_locator(self, tolerance=None): try: return cache[tolerance] except KeyError: - 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, size_t ncells_ignore, int* cells_ignore) + int locator(struct Function *f, double *x, {RealType_c} *X, {RealType_c} *ref_cell_dists_l1, {IntType_c} *cells, {IntType_c} npoints, size_t ncells_ignore, {IntType_c}* cells_ignore) {{ - {IntTypeC} j = 0; /* index into x and X */ - for({IntTypeC} i=0; icoords, f->coords_map); return cell_dist_l1; } -%(RealType)s to_reference_coords_xtr(void *result_, struct Function *f, int cell, int layer, double *x) +%(RealType)s to_reference_coords_xtr(void *result_, struct Function *f, %(IntType)s cell, %(IntType)s layer, double *x) { %(RealType)s cell_dist_l1 = 0.0; %(non_extr_comment_out)s%(IntType)s layers[2] = {0, layer+2}; // +2 because the layer loop goes to layers[1]-1, which is nlayers-1 diff --git a/firedrake/preconditioners/fdm.py b/firedrake/preconditioners/fdm.py index 56468e4a1e..89040eda7a 100644 --- a/firedrake/preconditioners/fdm.py +++ b/firedrake/preconditioners/fdm.py @@ -1023,9 +1023,9 @@ class SchurComplementBlockLU(SchurComplementKernel): """Schur complement kernel builder that assumes a block-diagonal interior block, and uses its LU factorization to compute S = A11 - (A10 U^-1) (L^-1 A01).""" condense_code = dedent(""" - PetscBLASInt bn, lierr, lwork; + PetscBLASInt bn, lierr, lwork, *ipiv; PetscBool done; - PetscInt m, bsize, irow, icol, nnz, iswap, *ipiv, *perm; + PetscInt m, bsize, irow, icol, nnz, iswap, *perm; const PetscInt *ai; PetscScalar *vals, *work, *L, *U; Mat X; @@ -1123,9 +1123,9 @@ class SchurComplementBlockInverse(SchurComplementKernel): """Schur complement kernel builder that assumes a block-diagonal interior block, and uses its inverse to compute S = A11 - A10 A00^-1 A01.""" condense_code = dedent(""" - PetscBLASInt bn, lierr, lwork; + PetscBLASInt bn, lierr, lwork, *ipiv; PetscBool done; - PetscInt m, irow, bsize, *ipiv; + PetscInt m, irow, bsize; const PetscInt *ai; PetscScalar *vals, *work, *ainv, swork; PetscCall(MatProductNumeric(A11)); diff --git a/pyop2/codegen/builder.py b/pyop2/codegen/builder.py index d05e5aeeae..97adee8ab7 100644 --- a/pyop2/codegen/builder.py +++ b/pyop2/codegen/builder.py @@ -51,14 +51,14 @@ def __init__(self, interior_horizontal, layer_bounds, num_layers, values = Argument(shape, dtype=dtype, pfx="map") if offset is not None: assert type(offset) == tuple - offset = numpy.array(offset, dtype=numpy.int32) + offset = numpy.array(offset, dtype=dtype) if len(set(offset)) == 1: offset = Literal(offset[0], casting=True) else: offset = NamedLiteral(offset, parent=values, suffix="offset") if offset_quotient is not None: assert type(offset_quotient) == tuple - offset_quotient = numpy.array(offset_quotient, dtype=numpy.int32) + offset_quotient = numpy.array(offset_quotient, dtype=dtype) offset_quotient = NamedLiteral(offset_quotient, parent=values, suffix="offset_quotient") self.values = values