Petsc 64 bit int fixes

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,
-			   size_t ncells_ignore,
-			   int* cells_ignore);
+		       PetscReal *found_ref_cell_dist_l1,
+			   PetscInt ncells_ignore,
+			   PetscInt* cells_ignore);
 
 extern int evaluate(struct Function *f,
 		    double *x,

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),

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

firedrake/locate.c

@@ -4,19 +4,19 @@
 #include <float.h>
 #include <evaluate.h>
 
-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,
-        size_t ncells_ignore,
-        int* cells_ignore)
+        PetscReal *found_ref_cell_dist_l1,
+        PetscInt ncells_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 (PetscInt 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 (PetscInt 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 (PetscInt 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 (PetscInt j = 0; j < ncells_ignore; j++) {
                         if (l == cells_ignore[j]) {
                             cell_ignore_found = 1;
                             break;

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)
+                {IntType_c} locator(struct Function *f, double *x, {RealType_c} *X, {RealType_c} *ref_cell_dists_l1, {IntType_c} *cells, {IntType_c} npoints, {IntType_c} ncells_ignore, {IntType_c}* cells_ignore)
                 {{
-                    {IntTypeC} j = 0;  /* index into x and X */
-                    for({IntTypeC} i=0; i<npoints; i++) {{
+                    {IntType_c} j = 0;  /* index into x and X */
+                    for({IntType_c} i=0; i<npoints; i++) {{
                         /* i is the index into cells and ref_cell_dists_l1 */
 
                         /* The type definitions and arguments used here are defined as
@@ -2758,7 +2758,7 @@ def _c_locator(self, tolerance=None):
                         pointquery_utils.py. If they contain python calls, this loop will
                         not run at c-loop speed. */
                         /* cells_ignore has shape (npoints, ncells_ignore) - find the ith row */
-                        int *cells_ignore_i = cells_ignore + i*ncells_ignore;
+                        {IntType_c} *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++) {{
@@ -2791,13 +2791,13 @@ def _c_locator(self, tolerance=None):
             locator = getattr(dll, "locator")
             locator.argtypes = [ctypes.POINTER(function._CFunction),
                                 ctypes.POINTER(ctypes.c_double),
-                                ctypes.POINTER(ctypes.c_double),
-                                ctypes.POINTER(ctypes.c_double),
-                                ctypes.POINTER(ctypes.c_int),
-                                ctypes.c_size_t,
-                                ctypes.c_size_t,
-                                np.ctypeslib.ndpointer(ctypes.c_int, flags="C_CONTIGUOUS")]
-            locator.restype = ctypes.c_int
+                                ctypes.POINTER(as_ctypes(RealType)),
+                                ctypes.POINTER(as_ctypes(RealType)),
+                                ctypes.POINTER(as_ctypes(IntType)),
+                                as_ctypes(IntType),
+                                as_ctypes(IntType),
+                                np.ctypeslib.ndpointer(as_ctypes(IntType), flags="C_CONTIGUOUS")]
+            locator.restype = as_ctypes(IntType)
             return cache.setdefault(tolerance, locator)
 
     @cached_property  # TODO: Recalculate if mesh moves. Extend this for regular meshes.

firedrake/pointeval_utils.py

@@ -160,6 +160,7 @@ def predicate(index):
         "extruded_define": "1" if extruded else "0",
         "IntType": as_cstr(IntType),
         "scalar_type": utils.ScalarType_c,
+        "real_type": utils.RealType_c,
     }
     # if maps are the same, only need to pass one of them
     if coordinates.cell_node_map() == coefficient.cell_node_map():
@@ -177,9 +178,9 @@ def predicate(index):
 int evaluate(struct Function *f, double *x, %(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;
+    %(real_type)s found_ref_cell_dist_l1 = PETSC_MAX_REAL;
     struct ReferenceCoords temp_reference_coords, found_reference_coords;
-    int cells_ignore[1] = {-1};
+    %(IntType)s 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;

firedrake/pointquery_utils.py

@@ -283,14 +283,14 @@ def compile_coordinate_element(mesh: MeshGeometry, contains_eps: float, paramete
     void* const result_, double* const x, %(RealType)s* const cell_dist_l1, %(IntType)s const start, %(IntType)s const end%(extruded_arg)s,
     %(ScalarType)s const *__restrict__ coords, %(IntType)s const *__restrict__ coords_map);
 
-%(RealType)s to_reference_coords(void *result_, struct Function *f, int cell, double *x)
+%(RealType)s to_reference_coords(void *result_, struct Function *f, %(IntType)s cell, double *x)
 {
     %(RealType)s cell_dist_l1 = 0.0;
     %(extr_comment_out)swrap_to_reference_coords(result_, x, &cell_dist_l1, cell, cell+1, f->coords, 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

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));

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