Merge pull request #1875 from OceanParcels/remove-jit-compat

Remove cstruct objects and remove c-contiguous requirement

Author: VeckoTheGecko

parcels/field.py

@@ -2,7 +2,7 @@
 import math
 import warnings
 from collections.abc import Iterable
-from ctypes import POINTER, Structure, c_float, c_int, pointer
+from ctypes import c_int
 from pathlib import Path
 from typing import TYPE_CHECKING, Literal
 
@@ -50,11 +50,9 @@
     DeferredNetcdfFileBuffer,
     NetcdfFileBuffer,
 )
-from .grid import CGrid, Grid, GridType, _calc_cell_areas
+from .grid import Grid, GridType, _calc_cell_areas
 
 if TYPE_CHECKING:
-    from ctypes import _Pointer as PointerType
-
     from parcels.fieldset import FieldSet
 
 __all__ = ["Field", "NestedField", "VectorField"]
@@ -341,7 +339,6 @@ def __init__(
         # since some datasets do not provide the deeper level of data (which is ignored by the interpolation).
         self.data_full_zdim = kwargs.pop("data_full_zdim", None)
         self._data_chunks = []  # type: ignore # the data buffer of the FileBuffer raw loaded data - shall be a list of C-contiguous arrays
-        self._c_data_chunks: list[PointerType | None] = []  # C-pointers to the data_chunks array
         self.nchunks: tuple[int, ...] = ()
         self._chunk_set: bool = False
         self.filebuffers = [None] * 2
@@ -1016,10 +1013,7 @@ def _time_index(self, time):
                 periods = int(
                     math.floor((time - self.grid.time_full[0]) / (self.grid.time_full[-1] - self.grid.time_full[0]))
                 )
-                if isinstance(self.grid.periods, c_int):
-                    self.grid.periods.value = periods
-                else:
-                    self.grid.periods = periods
+                self.grid.periods = periods
                 time -= periods * (self.grid.time_full[-1] - self.grid.time_full[0])
                 time_index = self.grid.time <= time
                 ti = time_index.argmin() - 1 if time_index.any() else 0
@@ -1106,7 +1100,6 @@ def _chunk_setup(self):
             return
 
         self._data_chunks = [None] * npartitions
-        self._c_data_chunks = [None] * npartitions
         self.grid._load_chunk = np.zeros(npartitions, dtype=c_int, order="C")
         # self.grid.chunk_info format: number of dimensions (without tdim); number of chunks per dimensions;
         #      chunksizes (the 0th dim sizes for all chunk of dim[0], then so on for next dims
@@ -1136,62 +1129,14 @@ def _chunk_data(self):
                         self._data_chunks[block_id] = None
                     else:
                         self._data_chunks[block_id, :] = None
-                    self._c_data_chunks[block_id] = None
         else:
             if isinstance(self._data_chunks, list):
                 self._data_chunks[0] = None
             else:
                 self._data_chunks[0, :] = None
-            self._c_data_chunks[0] = None
             self.grid._load_chunk[0] = g._chunk_loaded_touched
             self._data_chunks[0] = np.array(self.data, order="C")
 
-    @property
-    def ctypes_struct(self):
-        """Returns a ctypes struct object containing all relevant pointers and sizes for this field."""
-
-        # Ctypes struct corresponding to the type definition in parcels.h
-        class CField(Structure):
-            _fields_ = [
-                ("xdim", c_int),
-                ("ydim", c_int),
-                ("zdim", c_int),
-                ("tdim", c_int),
-                ("igrid", c_int),
-                ("allow_time_extrapolation", c_int),
-                ("time_periodic", c_int),
-                ("data_chunks", POINTER(POINTER(POINTER(c_float)))),
-                ("grid", POINTER(CGrid)),
-            ]
-
-        # Create and populate the c-struct object
-        allow_time_extrapolation = 1 if self.allow_time_extrapolation else 0
-        time_periodic = 1 if self.time_periodic else 0
-        for i in range(len(self.grid._load_chunk)):
-            if self.grid._load_chunk[i] == self.grid._chunk_loading_requested:
-                raise ValueError(
-                    "data_chunks should have been loaded by now if requested. grid._load_chunk[bid] cannot be 1"
-                )
-            if self.grid._load_chunk[i] in self.grid._chunk_loaded:
-                if not self._data_chunks[i].flags["C_CONTIGUOUS"]:
-                    self._data_chunks[i] = np.array(self._data_chunks[i], order="C")
-                self._c_data_chunks[i] = self._data_chunks[i].ctypes.data_as(POINTER(POINTER(c_float)))
-            else:
-                self._c_data_chunks[i] = None
-
-        cstruct = CField(
-            self.grid.xdim,
-            self.grid.ydim,
-            self.grid.zdim,
-            self.grid.tdim,
-            self.igrid,
-            allow_time_extrapolation,
-            time_periodic,
-            (POINTER(POINTER(c_float)) * len(self._c_data_chunks))(*self._c_data_chunks),
-            pointer(self.grid.ctypes_struct),
-        )
-        return cstruct
-
     def add_periodic_halo(self, zonal, meridional, halosize=5, data=None):
         """Add a 'halo' to all Fields in a FieldSet.
 

parcels/grid.py

@@ -1,6 +1,5 @@
 import functools
 import warnings
-from ctypes import POINTER, Structure, c_double, c_float, c_int, c_void_p, cast, pointer
 from enum import IntEnum
 
 import numpy as np
@@ -11,7 +10,6 @@
 from parcels.tools.warnings import FieldSetWarning
 
 __all__ = [
-    "CGrid",
     "CurvilinearSGrid",
     "CurvilinearZGrid",
     "Grid",
@@ -34,10 +32,6 @@ class GridType(IntEnum):
 GridCode = GridType
 
 
-class CGrid(Structure):
-    _fields_ = [("gtype", c_int), ("grid", c_void_p)]
-
-
 class Grid:
     """Grid class that defines a (spatial and temporal) grid on which Fields are defined."""
 
@@ -51,13 +45,10 @@ def __init__(
     ):
         self._ti = -1
         self._update_status: UpdateStatus | None = None
-        if not lon.flags["C_CONTIGUOUS"]:
-            lon = np.array(lon, order="C")
-        if not lat.flags["C_CONTIGUOUS"]:
-            lat = np.array(lat, order="C")
+        lon = np.array(lon)
+        lat = np.array(lat)
         time = np.zeros(1, dtype=np.float64) if time is None else time
-        if not time.flags["C_CONTIGUOUS"]:
-            time = np.array(time, order="C")
+        time = np.array(time)
         if not lon.dtype == np.float32:
             lon = lon.astype(np.float32)
         if not lat.dtype == np.float32:
@@ -76,7 +67,6 @@ def __init__(
         assert isinstance(self.time_origin, TimeConverter), "time_origin needs to be a TimeConverter object"
         assert_valid_mesh(mesh)
         self._mesh = mesh
-        self._cstruct = None
         self._cell_edge_sizes: dict[str, npt.NDArray] = {}
         self._zonal_periodic = False
         self._zonal_halo = 0
@@ -179,67 +169,6 @@ def create_grid(
             else:
                 return CurvilinearSGrid(lon, lat, depth, time, time_origin=time_origin, mesh=mesh, **kwargs)
 
-    @property
-    def ctypes_struct(self):
-        # This is unnecessary for the moment, but it could be useful when going will fully unstructured grids
-        self._cgrid = cast(pointer(self._child_ctypes_struct), c_void_p)
-        cstruct = CGrid(self._gtype, self._cgrid.value)
-        return cstruct
-
-    @property
-    def _child_ctypes_struct(self):
-        """Returns a ctypes struct object containing all relevant
-        pointers and sizes for this grid.
-        """
-
-        class CStructuredGrid(Structure):
-            # z4d is only to have same cstruct as RectilinearSGrid
-            _fields_ = [
-                ("xdim", c_int),
-                ("ydim", c_int),
-                ("zdim", c_int),
-                ("tdim", c_int),
-                ("z4d", c_int),
-                ("mesh_spherical", c_int),
-                ("zonal_periodic", c_int),
-                ("chunk_info", POINTER(c_int)),
-                ("load_chunk", POINTER(c_int)),
-                ("tfull_min", c_double),
-                ("tfull_max", c_double),
-                ("periods", POINTER(c_int)),
-                ("lonlat_minmax", POINTER(c_float)),
-                ("lon", POINTER(c_float)),
-                ("lat", POINTER(c_float)),
-                ("depth", POINTER(c_float)),
-                ("time", POINTER(c_double)),
-            ]
-
-        # Create and populate the c-struct object
-        if not self._cstruct:  # Not to point to the same grid various times if grid in various fields
-            if not isinstance(self.periods, c_int):
-                self.periods = c_int()
-                self.periods.value = 0
-            self._cstruct = CStructuredGrid(
-                self.xdim,
-                self.ydim,
-                self.zdim,
-                self.tdim,
-                self._z4d,
-                int(self.mesh == "spherical"),
-                int(self.zonal_periodic),
-                (c_int * len(self.chunk_info))(*self.chunk_info),
-                self._load_chunk.ctypes.data_as(POINTER(c_int)),
-                self.time_full[0],
-                self.time_full[-1],
-                pointer(self.periods),
-                self.lonlat_minmax.ctypes.data_as(POINTER(c_float)),
-                self.lon.ctypes.data_as(POINTER(c_float)),
-                self.lat.ctypes.data_as(POINTER(c_float)),
-                self.depth.ctypes.data_as(POINTER(c_float)),
-                self.time.ctypes.data_as(POINTER(c_double)),
-            )
-        return self._cstruct
-
     def _check_zonal_periodic(self):
         if self.zonal_periodic or self.mesh == "flat" or self.lon.size == 1:
             return
@@ -278,7 +207,7 @@ def _add_Sdepth_periodic_halo(self, zonal, meridional, halosize):
 
     def _computeTimeChunk(self, f, time, signdt):
         nextTime_loc = np.inf if signdt >= 0 else -np.inf
-        periods = self.periods.value if isinstance(self.periods, c_int) else self.periods
+        periods = self.periods
         prev_time_indices = self.time
         if self._update_status == "not_updated":
             if self._ti >= 0:
@@ -316,7 +245,7 @@ def _computeTimeChunk(self, f, time, signdt):
             if self._ti == -1:
                 self.time = self.time_full
                 self._ti, _ = f._time_index(time)
-                periods = self.periods.value if isinstance(self.periods, c_int) else self.periods
+                periods = self.periods
                 if (
                     signdt == -1
                     and self._ti == 0
@@ -483,8 +412,7 @@ def __init__(self, lon, lat, depth=None, time=None, time_origin=None, mesh: Mesh
 
         self._gtype = GridType.RectilinearZGrid
         self._depth = np.zeros(1, dtype=np.float32) if depth is None else depth
-        if not self.depth.flags["C_CONTIGUOUS"]:
-            self._depth = np.array(self.depth, order="C")
+        self._depth = np.array(self.depth)
         self._z4d = -1  # only used in RectilinearSGrid
         if not self.depth.dtype == np.float32:
             self._depth = self.depth.astype(np.float32)
@@ -539,8 +467,7 @@ def __init__(
 
         self._gtype = GridType.RectilinearSGrid
         self._depth = depth
-        if not self.depth.flags["C_CONTIGUOUS"]:
-            self._depth = np.array(self.depth, order="C")
+        self._depth = np.array(self.depth)
         self._z4d = 1 if len(self.depth.shape) == 4 else 0
         if self._z4d:
             # self.depth.shape[0] is 0 for S grids loaded from netcdf file
@@ -656,8 +583,6 @@ def __init__(
 
         self._gtype = GridType.CurvilinearZGrid
         self._depth = np.zeros(1, dtype=np.float32) if depth is None else depth
-        if not self.depth.flags["C_CONTIGUOUS"]:
-            self._depth = np.array(self.depth, order="C")
         self._z4d = -1  # only for SGrid
         if not self.depth.dtype == np.float32:
             self._depth = self.depth.astype(np.float32)
@@ -710,9 +635,7 @@ def __init__(
         assert isinstance(depth, np.ndarray) and len(depth.shape) in [3, 4], "depth is not a 4D numpy array"
 
         self._gtype = GridType.CurvilinearSGrid
-        self._depth = depth  # should be a C-contiguous array of floats
-        if not self.depth.flags["C_CONTIGUOUS"]:
-            self._depth = np.array(self.depth, order="C")
+        self._depth = depth
         self._z4d = 1 if len(self.depth.shape) == 4 else 0
         if self._z4d:
             # self.depth.shape[0] is 0 for S grids loaded from netcdf file

parcels/particledata.py

@@ -1,5 +1,4 @@
 import warnings
-from ctypes import POINTER, Structure
 from operator import attrgetter
 
 import numpy as np
@@ -42,8 +41,7 @@ def __init__(self, pclass, lon, lat, depth, time, lonlatdepth_dtype, pid_orig, n
         Parameters
         ----------
         ngrid :
-            number of grids in the fieldset of the overarching ParticleSet - required for initialising the
-            field references of the ctypes-link of particles that are allocated
+            number of grids in the fieldset of the overarching ParticleSet.
         """
         self._ncount = -1
         self._pu_indicators = None
@@ -328,21 +326,6 @@ def remove_multi_by_indices(self, indices):
 
         self._ncount -= len(indices)
 
-    def cstruct(self):
-        """Return the ctypes mapping of the particle data."""
-
-        class CParticles(Structure):
-            _fields_ = [(v.name, POINTER(np.ctypeslib.as_ctypes_type(v.dtype))) for v in self._ptype.variables]
-
-        def flatten_dense_data_array(vname):
-            data_flat = self._data[vname].view()
-            data_flat.shape = -1
-            return np.ctypeslib.as_ctypes(data_flat)
-
-        cdata = [flatten_dense_data_array(v.name) for v in self._ptype.variables]
-        cstruct = CParticles(*cdata)
-        return cstruct
-
     def _to_write_particles(self, time):
         """Return the Particles that need to be written at time: if particle.time is between time-dt/2 and time+dt (/2)"""
         pd = self._data

parcels/particleset.py

@@ -321,14 +321,6 @@ def lonlatdepth_dtype_from_field_interp_method(field):
                 return np.float64
         return np.float32
 
-    def cstruct(self):
-        cstruct = self.particledata.cstruct()
-        return cstruct
-
-    @property
-    def ctypes_struct(self):
-        return self.cstruct()
-
     @property
     def size(self):
         # ==== to change at some point - len and size are different things ==== #