From 2275e2cbaeccd96309746d9a1cf096195dbff9ef Mon Sep 17 00:00:00 2001 From: Christian Kehl Date: Tue, 12 May 2020 13:34:17 +0200 Subject: [PATCH 1/4] avoid unnecessary ndarray copies by initialising them in C-contiguous order - to be tested --- parcels/field.py | 14 +++++++------- parcels/grid.py | 24 +++++++++++++++++++----- parcels/kernel/basekernel.py | 8 ++++---- 3 files changed, 30 insertions(+), 16 deletions(-) diff --git a/parcels/field.py b/parcels/field.py index ce1c27e635..896c4c5d57 100644 --- a/parcels/field.py +++ b/parcels/field.py @@ -198,8 +198,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 @@ -1142,7 +1142,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:]), [])] @@ -1158,7 +1158,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 @@ -1172,7 +1172,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): @@ -1195,8 +1195,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 861b59a06b..e0690083b8 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: @@ -62,7 +68,7 @@ def __init__(self, lon, lat, time, time_origin, mesh): self.defer_load = False self.lonlat_minmax = np.array([np.nanmin(lon), np.nanmax(lon), np.nanmin(lat), np.nanmax(lat)], dtype=np.float32) self.periods = 0 - self.load_chunk = [] + self.load_chunk = [] # should be a C-contiguous array of ints self.chunk_info = None self.chunksize = None self._add_last_periodic_data_timestep = False @@ -331,7 +337,9 @@ def __init__(self, lon, lat, depth=None, time=None, time_origin=None, mesh='flat assert(len(depth.shape) <= 1), 'depth is not a vector' self.gtype = GridCode.RectilinearZGrid - self.depth = np.zeros(1, dtype=np.float32) if depth is None else depth + self.depth = np.zeros(1, dtype=np.float32) if depth is None else 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.size self.z4d = -1 # only used in RectilinearSGrid if not self.depth.dtype == np.float32: @@ -366,7 +374,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 3D or 4D numpy array' self.gtype = GridCode.RectilinearSGrid - 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 = len(self.depth.shape) == 4 if self.z4d: @@ -461,7 +471,9 @@ def __init__(self, lon, lat, depth=None, time=None, time_origin=None, mesh='flat assert(len(depth.shape) == 1), 'depth is not a vector' self.gtype = GridCode.CurvilinearZGrid - self.depth = np.zeros(1, dtype=np.float32) if depth is None else depth + self.depth = np.zeros(1, dtype=np.float32) if depth is None else 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.size self.z4d = -1 # only for SGrid if not self.depth.dtype == np.float32: @@ -495,7 +507,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 = len(self.depth.shape) == 4 if self.z4d: diff --git a/parcels/kernel/basekernel.py b/parcels/kernel/basekernel.py index cb11a7df8c..f80bae6645 100644 --- a/parcels/kernel/basekernel.py +++ b/parcels/kernel/basekernel.py @@ -334,13 +334,13 @@ def load_fieldset_jit(self, pset): 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() + 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): """ From 54929dd022b265ed64d0ee9f571a1b52a5183da1 Mon Sep 17 00:00:00 2001 From: Erik van Sebille Date: Tue, 17 Aug 2021 13:35:06 +0200 Subject: [PATCH 2/4] Cleaning up some stray comments --- parcels/grid.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/parcels/grid.py b/parcels/grid.py index e0690083b8..8dad136fcc 100644 --- a/parcels/grid.py +++ b/parcels/grid.py @@ -68,7 +68,7 @@ def __init__(self, lon, lat, time, time_origin, mesh): self.defer_load = False self.lonlat_minmax = np.array([np.nanmin(lon), np.nanmax(lon), np.nanmin(lat), np.nanmax(lat)], dtype=np.float32) self.periods = 0 - self.load_chunk = [] # should be a C-contiguous array of ints + self.load_chunk = [] self.chunk_info = None self.chunksize = None self._add_last_periodic_data_timestep = False @@ -337,7 +337,7 @@ def __init__(self, lon, lat, depth=None, time=None, time_origin=None, mesh='flat assert(len(depth.shape) <= 1), 'depth is not a vector' self.gtype = GridCode.RectilinearZGrid - self.depth = np.zeros(1, dtype=np.float32) if depth is None else depth # should be a C-contiguous array of floats + 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 @@ -374,7 +374,7 @@ 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 3D or 4D numpy array' self.gtype = GridCode.RectilinearSGrid - self.depth = depth # should be a C-contiguous array of floats + self.depth = depth if not self.depth.flags['C_CONTIGUOUS']: self.depth = np.array(self.depth, order='C') self.zdim = self.depth.shape[-3] @@ -471,7 +471,7 @@ def __init__(self, lon, lat, depth=None, time=None, time_origin=None, mesh='flat assert(len(depth.shape) == 1), 'depth is not a vector' self.gtype = GridCode.CurvilinearZGrid - self.depth = np.zeros(1, dtype=np.float32) if depth is None else depth # should be a C-contiguous array of floats + 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 From ce9a855f864d925f70bdc512fa2575288e159632 Mon Sep 17 00:00:00 2001 From: Erik van Sebille Date: Tue, 4 Oct 2022 08:02:49 +0200 Subject: [PATCH 3/4] Changing 3D axes generation to fix windows plotting issue --- parcels/scripts/plottrajectoriesfile.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/parcels/scripts/plottrajectoriesfile.py b/parcels/scripts/plottrajectoriesfile.py index 0d181c4e4c..2f21e6539f 100644 --- a/parcels/scripts/plottrajectoriesfile.py +++ b/parcels/scripts/plottrajectoriesfile.py @@ -71,7 +71,7 @@ def plotTrajectoriesFile(filename, mode='2d', tracerfile=None, tracerfield='P', if mode == '3d': from mpl_toolkits.mplot3d import Axes3D # noqa plt.clf() # clear the figure - ax = fig.gca(projection='3d') + ax = plt.axes(projection='3d') for p in range(len(lon)): ax.plot(lon[p, :], lat[p, :], z[p, :], '.-') ax.set_xlabel('Longitude') From 9fc69a45e3cbfa0a51d37213b66a2399dc234ca9 Mon Sep 17 00:00:00 2001 From: Erik van Sebille Date: Sat, 25 Feb 2023 14:31:56 +0100 Subject: [PATCH 4/4] Update parcels/kernel/basekernel.py --- parcels/kernel/basekernel.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/parcels/kernel/basekernel.py b/parcels/kernel/basekernel.py index 266e43a072..4cc7a5bcf1 100644 --- a/parcels/kernel/basekernel.py +++ b/parcels/kernel/basekernel.py @@ -338,7 +338,7 @@ 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: + 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 = np.array(g.depth, order='C')