diff --git a/parcels/field.py b/parcels/field.py index de2dc7051c..0d3c8fb428 100644 --- a/parcels/field.py +++ b/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 diff --git a/parcels/grid.py b/parcels/grid.py index f955360fb8..891af2ecc9 100644 --- a/parcels/grid.py +++ b/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: diff --git a/parcels/kernel/basekernel.py b/parcels/kernel/basekernel.py index f72d38578c..4cc7a5bcf1 100644 --- a/parcels/kernel/basekernel.py +++ b/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): """