Speedup C-contiguous arrays

parcels/field.py

@@ -206,8 +206,8 @@ def __init__(self, name, data, lon=None, lat=None, depth=None, time=None, grid=N
         # (data_full_zdim = grid.zdim if no indices are used, for A- and C-grids and for some B-grids). It is used for the B-grid,
         # 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 = []
-        self.c_data_chunks = []
+        self.data_chunks = []   # the data buffer of the FileBuffer raw loaded data - shall be a list of C-contiguous arrays
+        self.c_data_chunks = []  # C-pointers to the data_chunks array
         self.nchunks = []
         self.chunk_set = False
         self.filebuffers = [None] * 2
@@ -1182,7 +1182,7 @@ def chunk_setup(self):
 
         self.data_chunks = [None] * npartitions
         self.c_data_chunks = [None] * npartitions
-        self.grid.load_chunk = np.zeros(npartitions, dtype=c_int)
+        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
         self.grid.chunk_info = [[len(self.nchunks)-1], list(self.nchunks[1:]), sum(list(list(ci) for ci in chunks[1:]), [])]
@@ -1198,7 +1198,7 @@ def chunk_data(self):
                 if g.load_chunk[block_id] == g.chunk_loading_requested \
                         or g.load_chunk[block_id] in g.chunk_loaded and self.data_chunks[block_id] is None:
                     block = self.get_block(block_id)
-                    self.data_chunks[block_id] = np.array(self.data.blocks[(slice(self.grid.tdim),) + block])
+                    self.data_chunks[block_id] = np.array(self.data.blocks[(slice(self.grid.tdim),) + block], order='C')
                 elif g.load_chunk[block_id] == g.chunk_not_loaded:
                     if isinstance(self.data_chunks, list):
                         self.data_chunks[block_id] = None
@@ -1212,7 +1212,7 @@ def chunk_data(self):
                 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)
+            self.data_chunks[0] = np.array(self.data, order='C')
 
     @property
     def ctypes_struct(self):
@@ -1235,8 +1235,8 @@ class CField(Structure):
             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] = self.data_chunks[i].copy()
+                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

parcels/grid.py

@@ -40,8 +40,14 @@ def __init__(self, lon, lat, time, time_origin, mesh):
         self.zi = None
         self.ti = -1
         self.lon = lon
+        if not self.lon.flags['C_CONTIGUOUS']:
+            self.lon = np.array(self.lon, order='C')
         self.lat = lat
+        if not self.lat.flags['C_CONTIGUOUS']:
+            self.lat = np.array(self.lat, order='C')
         self.time = np.zeros(1, dtype=np.float64) if time is None else time
+        if not self.time.flags['C_CONTIGUOUS']:
+            self.time = np.array(self.time, order='C')
         if not self.lon.dtype == np.float32:
             self.lon = self.lon.astype(np.float32)
         if not self.lat.dtype == np.float32:
@@ -317,6 +323,8 @@ def __init__(self, lon, lat, depth=None, time=None, time_origin=None, mesh='flat
 
         self.gtype = GridCode.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.zdim = self.depth.size
         self.z4d = -1  # only used in RectilinearSGrid
         if not self.depth.dtype == np.float32:
@@ -352,6 +360,8 @@ def __init__(self, lon, lat, depth, time=None, time_origin=None, mesh='flat'):
 
         self.gtype = GridCode.RectilinearSGrid
         self.depth = depth
+        if not self.depth.flags['C_CONTIGUOUS']:
+            self.depth = np.array(self.depth, order='C')
         self.zdim = self.depth.shape[-3]
         self.z4d = 1 if len(self.depth.shape) == 4 else 0
         if self.z4d:
@@ -447,6 +457,8 @@ def __init__(self, lon, lat, depth=None, time=None, time_origin=None, mesh='flat
 
         self.gtype = GridCode.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.zdim = self.depth.size
         self.z4d = -1  # only for SGrid
         if not self.depth.dtype == np.float32:
@@ -480,7 +492,9 @@ def __init__(self, lon, lat, depth, time=None, time_origin=None, mesh='flat'):
         assert (isinstance(depth, np.ndarray) and len(depth.shape) in [3, 4]), 'depth is not a 4D numpy array'
 
         self.gtype = GridCode.CurvilinearSGrid
-        self.depth = depth
+        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.zdim = self.depth.shape[-3]
         self.z4d = 1 if len(self.depth.shape) == 4 else 0
         if self.z4d:

parcels/kernel/basekernel.py

@@ -338,14 +338,14 @@ def load_fieldset_jit(self, pset):
                 g.load_chunk = np.where(g.load_chunk == g.chunk_loading_requested,
                                         g.chunk_loaded_touched, g.load_chunk)
                 if len(g.load_chunk) > g.chunk_not_loaded:  # not the case if a field in not called in the kernel
-                    if not g.load_chunk.flags.c_contiguous:
-                        g.load_chunk = g.load_chunk.copy()
+                    if not g.load_chunk.flags['C_CONTIGUOUS']:
+                        g.load_chunk = np.array(g.load_chunk, order='C')
                 if not g.depth.flags.c_contiguous:
-                    g.depth = g.depth.copy()
+                    g.depth = np.array(g.depth, order='C')
                 if not g.lon.flags.c_contiguous:
-                    g.lon = g.lon.copy()
+                    g.lon = np.array(g.lon, order='C')
                 if not g.lat.flags.c_contiguous:
-                    g.lat = g.lat.copy()
+                    g.lat = np.array(g.lat, order='C')
 
     def evaluate_particle(self, p, endtime, sign_dt, dt, analytical=False):
         """