diff --git a/parcels/field.py b/parcels/field.py index 43807a4518..41e9a08b02 100644 --- a/parcels/field.py +++ b/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"] @@ -343,7 +341,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 @@ -1018,10 +1015,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 @@ -1120,7 +1114,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 @@ -1150,62 +1143,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. diff --git a/parcels/grid.py b/parcels/grid.py index 560597c740..874be301cb 100644 --- a/parcels/grid.py +++ b/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 diff --git a/parcels/particledata.py b/parcels/particledata.py index 3d384dbf07..df8f9d99e2 100644 --- a/parcels/particledata.py +++ b/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 diff --git a/parcels/particleset.py b/parcels/particleset.py index 5585b4c9f4..9adf60ffca 100644 --- a/parcels/particleset.py +++ b/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 ==== #