From 6466cf1f017c7660ea1b37351c94e3f1bfed0f48 Mon Sep 17 00:00:00 2001 From: Christian Kehl Date: Tue, 7 Jul 2020 20:02:48 +0200 Subject: [PATCH 01/56] list_of_pset_array_trials - initial addenda for list-of-particles version --- parcels/particle.py | 4 + parcels/particleset.py | 304 +++++++++++++++++++++++++++++++++-------- 2 files changed, 250 insertions(+), 58 deletions(-) diff --git a/parcels/particle.py b/parcels/particle.py index ff093b13ba..1800b52b77 100644 --- a/parcels/particle.py +++ b/parcels/particle.py @@ -258,3 +258,7 @@ def __init__(self, *args, **kwargs): setattr(self, index, -1*np.ones((fieldset.gridset.size), dtype=np.int32)) setattr(self, index+'p', getattr(self, index).ctypes.data_as(c_void_p)) setattr(self, 'c'+index, getattr(self, index+'p').value) + + def update_cptr(self, new_ptr): + np.copyto(new_ptr, self._cptr) + self._cptr = new_ptr diff --git a/parcels/particleset.py b/parcels/particleset.py index 419666c97e..eb1f80212b 100644 --- a/parcels/particleset.py +++ b/parcels/particleset.py @@ -31,7 +31,6 @@ __all__ = ['ParticleSet'] - class ParticleSet(object): """Container class for storing particle and executing kernel over them. @@ -155,35 +154,138 @@ def convert_to_array(var): 'lon lat depth precision should be set to either np.float32 or np.float64' JITParticle.set_lonlatdepth_dtype(self.lonlatdepth_dtype) - self.particles = np.empty(lon.size, dtype=pclass) - self.ptype = pclass.getPType() + self._pid_mapping_bounds = {} # plist bracket_index -> (min_id, max_id, size_bracket) + self.nlist_limit = 4096 + # self.particles = np.empty(lon.size, dtype=pclass) + assert lon.shape[0] == lat.shape[0], ('Length of lon and lat do not match.') + self._plist = [] + start_index = 0 + end_index = min(self.nlist_limit , lon.shape[0]) + while end_index < lon.shape[0]: + self._plist.append(np.empty(end_index-start_index, dtype=pclass)) + bracket_index = len(self._plist)-1 + self._pid_mapping_bounds[bracket_index] = (np.iinfo(np.int32).max, np.iinfo(np.int32).min, end_index-start_index) + start_index = end_index + end_index = min(lon.shape[0], end_index+self.nlist_limit) + # self._particle_data = None + self._plist_c = None + self._pclass = pclass + self._ptype = self.pclass.getPType() self.kernel = None - if self.ptype.uses_jit: - # Allocate underlying data for C-allocated particles - self._particle_data = np.empty(lon.size, dtype=self.ptype.dtype) - - def cptr(i): - return self._particle_data[i] - else: - def cptr(i): - return None + # ==== ANNOTATION ==== # + # the particles themselves (self.particles, dtype: pclass) are already separate memory fields from + # their c-pointers (self._particle_data, dtype: PType.dtype [ie. padded memory block]). Thus, it makes little + # sense to require them being static over the whole ParticleSet lifetime. In other words: there should be a mechanism + # to update the particles, reallocate the c-pointer field, and then just value-copy the particle data into the c-memory. + # nparticles = self.allocate_cptrs() + if self._ptype.uses_jit: + self._plist_c = [] + for bracket_index in range(len(self._plist)): + self._plist_c.append(np.empty(self._pid_mapping_bounds[bracket_index][2], dtype=self.ptype.dtype)) if lon is not None and lat is not None: - # Initialise from arrays of lon/lat coordinates - assert self.particles.size == lon.size and self.particles.size == lat.size, ( - 'Size of ParticleSet does not match length of lon and lat.') + # == Initialise from arrays of lon/lat coordinates == # + # assert self.particles.size == lon.size and self.particles.size == lat.size, ('Size of ParticleSet does not match length of lon and lat.') for i in range(lon.size): - self.particles[i] = pclass(lon[i], lat[i], pid[i], fieldset=fieldset, depth=depth[i], cptr=cptr(i), time=time[i]) - # Set other Variables if provided + bracket_index = int(float(i)/self.nlist_limit) + slot_index = int(i % self.nlist_limit) + self._plist[bracket_index][slot_index] = pclass(lon[i], lat[i], pid[i], fieldset=fieldset, depth=depth[i], cptr=self.cptr(bracket_index, slot_index), time=time[i]) + bracket_info = self._pid_mapping_bounds[bracket_index] + self._pid_mapping_bounds[bracket_index] = (min(pid[i], bracket_info[0]), max(pid[i], bracket_info[1]), bracket_info[2]) + # self.particles[i] = pclass(lon[i], lat[i], pid[i], fieldset=fieldset, depth=depth[i], cptr=cptr(i), time=time[i]) + # == Set other Variables if provided == # for kwvar in kwargs: - if not hasattr(self.particles[i], kwvar): + if not hasattr(self._plist[bracket_index][slot_index], kwvar): raise RuntimeError('Particle class does not have Variable %s' % kwvar) - setattr(self.particles[i], kwvar, kwargs[kwvar][i]) + setattr(self._plist[bracket_index][slot_index], kwvar, kwargs[kwvar][i]) + # if not hasattr(self.particles[i], kwvar): + # raise RuntimeError('Particle class does not have Variable %s' % kwvar) + # setattr(self.particles[i], kwvar, kwargs[kwvar][i]) else: raise ValueError("Latitude and longitude required for generating ParticleSet") + @property + def particles(self): + return self._plist + + @property + def pid_mapping_bounds(self): + return self._pid_mapping_bounds + + @property + def particles_c(self): + return self._plist_c + + @property + def ptype(self): + return self._ptype + + @property + def pclass(self): + return self._pclass + + def cptr(self, bracket_index, slot_index): + # if self._particle_data is None: + if self._plist_c is None: + return None + if self.ptype.uses_jit: + return self._plist_c[bracket_index][slot_index] + else: + return None + + # def allocate_cptrs(self): + # nparticles = 0 + # for sublist in self.plist: + # nparticles += sublist.shape[0] + # # == Allocate underlying data for C-allocated particles == # + # if self.ptype.uses_jit: + # self._particle_data = np.empty(nparticles, dtype=self.ptype.dtype) + # return nparticles + + def get_cptr_index(self, pdata): + if isinstance(pdata, int): + # == parameter is the particle ID == # + located_particle = [(p, bracket_index, slot_index) for bracket_index, sublist in enumerate(self.plist) for slot_index, p in enumerate(sublist) if p.id == pdata] + if len(located_particle) < 1: + logger.warn_once("ParticleSet.get_cptr_index() - requested particle ID ({}) not existent.".format(pdata)) + return -1 + if len(located_particle) > 1: + logger.warn_once("ParticleSet.get_cptr_index() - requested particle ID ({}) ambiguous.".format(pdata)) + return -1 + located_particle = located_particle[0] + nparticle = 0 + for i in range(located_particle[1]): + nparticle += self.pid_mapping_bounds[i][2] + return nparticle+located_particle[2] + elif isinstance(pdata, self.pclass): + # == parameter is a particle itself == # + nparticle = 0 + bracket_index = 0 + bracket_info = self.pid_mapping_bounds[bracket_index] + while bracket_index < len(self.plist) and ((pdata.id < bracket_info[0]) or (pdata.id > bracket_info[1])): + nparticle += bracket_info[2] + bracket_index += 1 + bracket_info = self.pid_mapping_bounds[bracket_index] + if bracket_index >= len(self.plist): + logger.warn_once("ParticleSet.get_cptr_index() - requested particle ({}) not found.".format(pdata)) + return -1 + slot_index = 0 + while (slot_index < self.plist[bracket_index].shape[0]) and (self.plist[bracket_index][slot_index].id != pdata.id): + slot_index += 1 + nparticle += 1 + if slot_index >= self.plist[bracket_index].shape[0]: + logger.warn_once("ParticleSet.get_cptr_index() - requested particle ({}) not found.".format(pdata)) + return -1 + return nparticle + elif isinstance(pdata, tuple): + # == parameter is a list-array index tuple with n=2 == # + nparticle = 0 + for i in range(len(pdata[0])): + nparticle += self.pid_mapping_bounds[i][2] + return nparticle+pdata[1] + @classmethod def from_list(cls, fieldset, pclass, lon, lat, depth=None, time=None, repeatdt=None, lonlatdepth_dtype=None, **kwargs): """Initialise the ParticleSet from lists of lon and lat @@ -334,21 +436,26 @@ def lonlatdepth_dtype_from_field_interp_method(field): return np.float64 return np.float32 - @property - def size(self): - return self.particles.size + # @property + # def size(self): + # nparticles = 0 + # for bracket_index in self._pid_mapping_bounds.keys(): + # nparticles += self._pid_mapping_bounds[bracket_index][2] + # return nparticles + # # return self.particles.size def __repr__(self): - return "\n".join([str(p) for p in self]) + return "\n".join([str(p) for sublist in self._plist for p in sublist]) - def __len__(self): - return self.size + # def __len__(self): + # return len(self._plist) + # # return self.size - def __getitem__(self, key): - return self.particles[key] + # def __getitem__(self, key): + # return self.particles[key] - def __setitem__(self, key, value): - self.particles[key] = value + # def __setitem__(self, key, value): + # self.particles[key] = value def __iadd__(self, particles): self.add(particles) @@ -368,30 +475,106 @@ def _create_progressbar_(self, starttime, endtime): def add(self, particles): """Method to add particles to the ParticleSet""" if isinstance(particles, ParticleSet): - particles = particles.particles + offset = len(self._plist) + self._plist += particles.particles + # self.pid_mapping_bounds.update(particles.pid_mapping_bounds) + for i in particles.pid_mapping_bounds.keys(): + self._pid_mapping_bounds[offset+i] = particles.pid_mapping_bounds[i] else: raise NotImplementedError('Only ParticleSets can be added to a ParticleSet') - self.particles = np.append(self.particles, particles) - if self.ptype.uses_jit: - particles_data = [p._cptr for p in particles] - self._particle_data = np.append(self._particle_data, particles_data) - # Update C-pointer on particles - for p, pdata in zip(self.particles, self._particle_data): - p._cptr = pdata + + for sublist in particles.particles_c: + self._plist_c.append(sublist) + + # if isinstance(particles, ParticleSet): + # particles = particles.particles + # else: + # raise NotImplementedError('Only ParticleSets can be added to a ParticleSet') + # self.particles = np.append(self.particles, particles) + # if self.ptype.uses_jit: + # particles_data = [p._cptr for p in particles] + # self._particle_data = np.append(self._particle_data, particles_data) + # # Update C-pointer on particles + # for p, pdata in zip(self.particles, self._particle_data): + # p._cptr = pdata + + def _merge_brackets_(self): + lw_bound_nlist = int(self.nlist_limit / 2) + nmerges = 1 + while nmerges >= 0: + nmerges = 0 + trg_bracket = None + trg_bracket_index = -1 + src_bracket = None + src_bracket_index = -1 + for bracket_index in range(len(self._plist)): + if self._pid_mapping_bounds[bracket_index][2] < lw_bound_nlist: + if trg_bracket is None: + trg_bracket = self._plist[bracket_index] + trg_bracket_index = bracket_index + elif src_bracket is None: + src_bracket = self._plist[bracket_index] + src_bracket_index = bracket_index + if src_bracket is not None and trg_bracket is not None: + break + if src_bracket is not None and trg_bracket is not None: + self._plist[trg_bracket_index] = np.append(trg_bracket, src_bracket) + local_ids = np.array([p.id for p in self._plist[trg_bracket_index]]) + self.pid_mapping_bounds[trg_bracket_index] = (np.min(local_ids), np.max(local_ids), local_ids.shape[0]) + if self.ptype.uses_jit: + self._plist_c[trg_bracket_index] = np.append(self._plist_c[trg_bracket_index], self._plist_c[src_bracket_index]) + for p, pdata in zip(self._plist[trg_bracket_index], self._plist_c[trg_bracket_index]): + p._cptr = pdata + self._plist.remove(self._plist[src_bracket_index]) + self._plist_c.remove(self._plist_c[src_bracket_index]) + self._pid_mapping_bounds.pop(src_bracket_index) def remove(self, indices): """Method to remove particles from the ParticleSet, based on their `indices`""" - if isinstance(indices, collections.Iterable): - particles = [self.particles[i] for i in indices] - else: - particles = self.particles[indices] - self.particles = np.delete(self.particles, indices) - if self.ptype.uses_jit: - self._particle_data = np.delete(self._particle_data, indices) - # Update C-pointer on particles - for p, pdata in zip(self.particles, self._particle_data): - p._cptr = pdata - return particles + assert len(indices) == len(self._plist), ("ParticleSet.remove() - particle set length ({}) does not match index list length ({})".format( + len(self._plist), len(indices))) + for bracket_index in range(len(indices)): + local_indices = indices[bracket_index] + if isinstance(local_indices, collections.Iterable): + local_indices = np.array(local_indices) + self._remove_(bracket_index, local_indices) + self._merge_brackets_() + + # if isinstance(indices, collections.Iterable): + # particles = [self.particles[i] for i in indices] + # else: + # particles = self.particles[indices] + # self.particles = np.delete(self.particles, indices) + # if self.ptype.uses_jit: + # self._particle_data = np.delete(self._particle_data, indices) + # # == Update C-pointer on particles == # + # for p, pdata in zip(self.particles, self._particle_data): + # p._cptr = pdata + # return particles + + def remove_local_particles(self, bracket_index, local_indices): + return self._remove_(bracket_index, local_indices) + + def _remove_(self, bracket_index, local_indices): + # if isinstance(indices, collections.Iterable): + # particles = [self.particles[i] for i in indices] + # else: + # particles = self.particles[indices] + if isinstance(local_indices, collections.Iterable): + local_indices = np.array(local_indices) + return_p = self._plist[bracket_index][local_indices] + self._plist[bracket_index] = np.delete(self._plist[bracket_index], local_indices) + local_ids = np.array([p.id for p in self._plist[bracket_index]]) + self.pid_mapping_bounds[bracket_index] = (np.min(local_ids), np.max(local_ids), local_ids.shape[0]) + if self._ptype.uses_jit: + self._plist_c[bracket_index] = np.delete(self._plist_c[bracket_index], local_indices) + # # == Update C-pointer on particles == # + # for p, pdata in zip(self._plist[bracket_index], self._plist_c[bracket_index]): + # p._cptr = pdata + for slot_index in range(len(self._plist[bracket_index])): + self._plist[bracket_index][slot_index].update_cptr(self._plist_c[bracket_index][slot_index]) + # return particles + return return_p def execute(self, pyfunc=AdvectionRK4, endtime=None, runtime=None, dt=1., moviedt=None, recovery=None, output_file=None, movie_background_field=None, @@ -433,7 +616,7 @@ def execute(self, pyfunc=AdvectionRK4, endtime=None, runtime=None, dt=1., else: self.kernel = self.Kernel(pyfunc) # Prepare JIT kernel execution - if self.ptype.uses_jit: + if self._ptype.uses_jit: self.kernel.remove_lib() cppargs = ['-DDOUBLE_COORD_VARIABLES'] if self.lonlatdepth_dtype == np.float64 else None self.kernel.compile(compiler=GNUCompiler(cppargs=cppargs)) @@ -465,10 +648,11 @@ def execute(self, pyfunc=AdvectionRK4, endtime=None, runtime=None, dt=1., assert moviedt is None or moviedt >= 0, 'moviedt must be positive' # Set particle.time defaults based on sign of dt, if not set at ParticleSet construction - for p in self: - if np.isnan(p.time): - mintime, maxtime = self.fieldset.gridset.dimrange('time_full') - p.time = mintime if dt >= 0 else maxtime + for sublist in self._plist: + for p in sublist: + if np.isnan(p.time): + mintime, maxtime = self.fieldset.gridset.dimrange('time_full') + p.time = mintime if dt >= 0 else maxtime # Derive _starttime and endtime from arguments or fieldset defaults if runtime is not None and endtime is not None: @@ -492,8 +676,9 @@ def execute(self, pyfunc=AdvectionRK4, endtime=None, runtime=None, dt=1., execute_once = True # Initialise particle timestepping - for p in self: - p.dt = dt + for sublist in self._plist: + for p in sublist: + p.dt = dt # First write output_file, because particles could have been added if output_file: @@ -524,6 +709,7 @@ def execute(self, pyfunc=AdvectionRK4, endtime=None, runtime=None, dt=1., if verbose_progress: pbar = self._create_progressbar_(_starttime, endtime) while (time < endtime and dt > 0) or (time > endtime and dt < 0) or dt == 0: + pbar = None if verbose_progress is None and time_module.time() - walltime_start > 10: # Showing progressbar if runtime > 10 seconds if output_file: @@ -542,8 +728,9 @@ def execute(self, pyfunc=AdvectionRK4, endtime=None, runtime=None, dt=1., lat=self.repeatlat, depth=self.repeatdepth, pclass=self.repeatpclass, lonlatdepth_dtype=self.lonlatdepth_dtype, partitions=False, pid_orig=self.repeatpid, **self.repeatkwargs) - for p in pset_new: - p.dt = dt + for sublist in pset_new.particles: + for p in sublist: + p.dt = dt self.add(pset_new) next_prelease += self.repeatdt * np.sign(dt) if abs(time-next_output) < tol: @@ -607,8 +794,9 @@ def density(self, field=None, particle_val=None, relative=False, area_scale=Fals field = field if field else self.fieldset.U if isinstance(particle_val, str): - particle_val = [getattr(p, particle_val) for p in self.particles] + particle_val = [np.array([getattr(p, particle_val) for p in sublist]) for sublist in self._plist] else: +# =================================================================================================================== # particle_val = particle_val if particle_val else np.ones(len(self.particles)) density = np.zeros((field.grid.lat.size, field.grid.lon.size), dtype=np.float32) From 9f52d80906e1a3cd53a224b1b7364e40c38a161d Mon Sep 17 00:00:00 2001 From: Christian Kehl Date: Wed, 8 Jul 2020 16:10:52 +0200 Subject: [PATCH 02/56] list_of_pset_array_trials - wrapped up list-of-array adaptations in the ParticleSet class. --- parcels/particleset.py | 123 +++++++++++++++++++++++++++-------------- 1 file changed, 83 insertions(+), 40 deletions(-) diff --git a/parcels/particleset.py b/parcels/particleset.py index eb1f80212b..ff742444e2 100644 --- a/parcels/particleset.py +++ b/parcels/particleset.py @@ -244,47 +244,78 @@ def cptr(self, bracket_index, slot_index): # self._particle_data = np.empty(nparticles, dtype=self.ptype.dtype) # return nparticles - def get_cptr_index(self, pdata): + def get_list_array_index(self, pdata): + """Searches for the list-array indices for a given particle. + + :param pdata: :mod:`int` global index or :mod:`parcels.particle._Particle` (or its subclass) object to search for + :return :mod:`tuple` in format (, ), to index + """ if isinstance(pdata, int): - # == parameter is the particle ID == # - located_particle = [(p, bracket_index, slot_index) for bracket_index, sublist in enumerate(self.plist) for slot_index, p in enumerate(sublist) if p.id == pdata] - if len(located_particle) < 1: - logger.warn_once("ParticleSet.get_cptr_index() - requested particle ID ({}) not existent.".format(pdata)) - return -1 - if len(located_particle) > 1: - logger.warn_once("ParticleSet.get_cptr_index() - requested particle ID ({}) ambiguous.".format(pdata)) - return -1 - located_particle = located_particle[0] - nparticle = 0 - for i in range(located_particle[1]): - nparticle += self.pid_mapping_bounds[i][2] - return nparticle+located_particle[2] + # == parameter is the global particle index == # + bracket_index = 0 + start_index = 0 + end_index = self._pid_mapping_bounds[bracket_index][2] + while start_index bracket_info[1])): - nparticle += bracket_info[2] + while bracket_index < len(self._plist) and ((pdata.id < bracket_info[0]) or (pdata.id > bracket_info[1])): + # nparticle += bracket_info[2] bracket_index += 1 bracket_info = self.pid_mapping_bounds[bracket_index] - if bracket_index >= len(self.plist): + if bracket_index >= len(self._plist): logger.warn_once("ParticleSet.get_cptr_index() - requested particle ({}) not found.".format(pdata)) return -1 slot_index = 0 - while (slot_index < self.plist[bracket_index].shape[0]) and (self.plist[bracket_index][slot_index].id != pdata.id): + while (slot_index < self._plist[bracket_index].shape[0]) and (self._plist[bracket_index][slot_index].id != pdata.id): slot_index += 1 - nparticle += 1 - if slot_index >= self.plist[bracket_index].shape[0]: + # nparticle += 1 + if slot_index >= self._plist[bracket_index].shape[0]: logger.warn_once("ParticleSet.get_cptr_index() - requested particle ({}) not found.".format(pdata)) return -1 - return nparticle - elif isinstance(pdata, tuple): - # == parameter is a list-array index tuple with n=2 == # - nparticle = 0 - for i in range(len(pdata[0])): - nparticle += self.pid_mapping_bounds[i][2] - return nparticle+pdata[1] + # return nparticle + return (bracket_index, slot_index) + + def get_list_array_index_by_PID(self, pdata): + # == parameter is the particle ID == # + located_particle = [(p, bracket_index, slot_index) for bracket_index, sublist in enumerate(self._plist) for slot_index, p in enumerate(sublist) if p.id == pdata] + if len(located_particle) < 1: + logger.warn_once("ParticleSet.get_cptr_index() - requested particle ID ({}) not existent.".format(pdata)) + return -1 + if len(located_particle) > 1: + logger.warn_once("ParticleSet.get_cptr_index() - requested particle ID ({}) ambiguous.".format(pdata)) + return -1 + located_particle = located_particle[0] + # nparticle = 0 + # for i in range(located_particle[1]): + # nparticle += self.pid_mapping_bounds[i][2] + # return nparticle+located_particle[2] + return (located_particle[1], located_particle[2]) + + def covert_static_array_to_list_array(self, property_array): + if not isinstance(property_array, np.ndarray): + property_array = np.array(property_array) + nparticles = 0 + for bracket_index in self._pid_mapping_bounds.keys(): + nparticles += self._pid_mapping_bounds[bracket_index][2] + assert property_array.shape[0] == nparticles, ("Attaching the property array to the list of particle prohibited because num. of entries do not match.") + results = [] + bracket_index = 0 + start_index = 0 + end_index = self._pid_mapping_bounds[bracket_index][2] + while start_index Date: Wed, 8 Jul 2020 16:48:14 +0200 Subject: [PATCH 03/56] list_of_pset_array_trials - adapted the kernel to work with list-of-array particles. TODO: particle file. --- parcels/kernel.py | 201 +++++++++++++++++++++++++--------------------- 1 file changed, 108 insertions(+), 93 deletions(-) diff --git a/parcels/kernel.py b/parcels/kernel.py index 95a785e4f6..7635ddfd74 100644 --- a/parcels/kernel.py +++ b/parcels/kernel.py @@ -219,9 +219,13 @@ def load_lib(self): def execute_jit(self, pset, endtime, dt): """Invokes JIT engine to perform the core update loop""" - if len(pset.particles) > 0: - assert pset.fieldset.gridset.size == len(pset.particles[0].xi), \ + if len(pset.particles) > 0 and pset.particles[0].shape[0]: + assert pset.fieldset.gridset.size == len(pset.particles[0][0].xi), \ 'FieldSet has different amount of grids than Particle.xi. Have you added Fields after creating the ParticleSet?' + + # if len(pset.particles) > 0: + # assert pset.fieldset.gridset.size == len(pset.particles[0].xi), \ + # 'FieldSet has different amount of grids than Particle.xi. Have you added Fields after creating the ParticleSet?' for g in pset.fieldset.gridset.grids: g.cstruct = None # This force to point newly the grids from Python to C # Make a copy of the transposed array to enforce @@ -250,11 +254,18 @@ def execute_jit(self, pset, endtime, dt): fargs = [byref(f.ctypes_struct) for f in self.field_args.values()] fargs += [c_double(f) for f in self.const_args.values()] - particle_data = pset._particle_data.ctypes.data_as(c_void_p) - self._function(c_int(len(pset)), particle_data, - c_double(endtime), - c_double(dt), - *fargs) + + for subfield_c in pset.particles_c: + particle_data = subfield_c.ctypes.data_as(c_void_p) + self._function(c_int(subfield_c.shape[0]), particle_data, + c_double(endtime), c_double(dt), + *fargs) + + # particle_data = pset._particle_data.ctypes.data_as(c_void_p) + # self._function(c_int(len(pset)), particle_data, + # c_double(endtime), + # c_double(dt), + # *fargs) def execute_python(self, pset, endtime, dt): """Performs the core update loop via Python""" @@ -268,104 +279,106 @@ def execute_python(self, pset, endtime, dt): continue f.data = np.array(f.data) - for p in pset.particles: - ptype = p.getPType() - # Don't execute particles that aren't started yet - sign_end_part = np.sign(endtime - p.time) + # for p in pset.particles: + for subfield in pset.particles: + for p in subfield: + ptype = p.getPType() + # Don't execute particles that aren't started yet + sign_end_part = np.sign(endtime - p.time) - dt_pos = min(abs(p.dt), abs(endtime - p.time)) + dt_pos = min(abs(p.dt), abs(endtime - p.time)) - # ==== numerically stable; also making sure that continuously-recovered particles do end successfully, - # as they fulfil the condition here on entering at the final calculation here. ==== # - if ((sign_end_part != sign_dt) or np.isclose(dt_pos, 0)) and not np.isclose(dt, 0): - if abs(p.time) >= abs(endtime): - p.state = ErrorCode.Success - continue + # ==== numerically stable; also making sure that continuously-recovered particles do end successfully, + # as they fulfil the condition here on entering at the final calculation here. ==== # + if ((sign_end_part != sign_dt) or np.isclose(dt_pos, 0)) and not np.isclose(dt, 0): + if abs(p.time) >= abs(endtime): + p.state = ErrorCode.Success + continue - while p.state in [ErrorCode.Evaluate, ErrorCode.Repeat] or np.isclose(dt, 0): - - for var in ptype.variables: - p_var_back[var.name] = getattr(p, var.name) - try: - pdt_prekernels = sign_dt * dt_pos - p.dt = pdt_prekernels - state_prev = p.state - res = self.pyfunc(p, pset.fieldset, p.time) - if res is None: - res = ErrorCode.Success - - if res is ErrorCode.Success and p.state != state_prev: - res = p.state - - if res == ErrorCode.Success and not np.isclose(p.dt, pdt_prekernels): - res = ErrorCode.Repeat - - except FieldOutOfBoundError as fse_xy: - res = ErrorCode.ErrorOutOfBounds - p.exception = fse_xy - except FieldOutOfBoundSurfaceError as fse_z: - res = ErrorCode.ErrorThroughSurface - p.exception = fse_z - except TimeExtrapolationError as fse_t: - res = ErrorCode.ErrorTimeExtrapolation - p.exception = fse_t - except Exception as e: - res = ErrorCode.Error - p.exception = e - - # Handle particle time and time loop - if res in [ErrorCode.Success, ErrorCode.Delete]: - # Update time and repeat - p.time += p.dt - p.update_next_dt() - dt_pos = min(abs(p.dt), abs(endtime - p.time)) - - sign_end_part = np.sign(endtime - p.time) - if res != ErrorCode.Delete and not np.isclose(dt_pos, 0) and (sign_end_part == sign_dt): - res = ErrorCode.Evaluate - if sign_end_part != sign_dt: - dt_pos = 0 - - p.state = res - if np.isclose(dt, 0): - break - else: - p.state = res - # Try again without time update - for var in ptype.variables: - if var.name not in ['dt', 'state']: - setattr(p, var.name, p_var_back[var.name]) - dt_pos = min(abs(p.dt), abs(endtime - p.time)) + while p.state in [ErrorCode.Evaluate, ErrorCode.Repeat] or np.isclose(dt, 0): - sign_end_part = np.sign(endtime - p.time) - if sign_end_part != sign_dt: - dt_pos = 0 - break + for var in ptype.variables: + p_var_back[var.name] = getattr(p, var.name) + try: + pdt_prekernels = sign_dt * dt_pos + p.dt = pdt_prekernels + state_prev = p.state + res = self.pyfunc(p, pset.fieldset, p.time) + if res is None: + res = ErrorCode.Success + + if res is ErrorCode.Success and p.state != state_prev: + res = p.state + + if res == ErrorCode.Success and not np.isclose(p.dt, pdt_prekernels): + res = ErrorCode.Repeat + + except FieldOutOfBoundError as fse_xy: + res = ErrorCode.ErrorOutOfBounds + p.exception = fse_xy + except FieldOutOfBoundSurfaceError as fse_z: + res = ErrorCode.ErrorThroughSurface + p.exception = fse_z + except TimeExtrapolationError as fse_t: + res = ErrorCode.ErrorTimeExtrapolation + p.exception = fse_t + except Exception as e: + res = ErrorCode.Error + p.exception = e + + # Handle particle time and time loop + if res in [ErrorCode.Success, ErrorCode.Delete]: + # Update time and repeat + p.time += p.dt + p.update_next_dt() + dt_pos = min(abs(p.dt), abs(endtime - p.time)) + + sign_end_part = np.sign(endtime - p.time) + if res != ErrorCode.Delete and not np.isclose(dt_pos, 0) and (sign_end_part == sign_dt): + res = ErrorCode.Evaluate + if sign_end_part != sign_dt: + dt_pos = 0 + + p.state = res + if np.isclose(dt, 0): + break + else: + p.state = res + # Try again without time update + for var in ptype.variables: + if var.name not in ['dt', 'state']: + setattr(p, var.name, p_var_back[var.name]) + dt_pos = min(abs(p.dt), abs(endtime - p.time)) + + sign_end_part = np.sign(endtime - p.time) + if sign_end_part != sign_dt: + dt_pos = 0 + break def remove_deleted(self, pset, output_file, endtime): """Utility to remove all particles that signalled deletion""" - indices = [i for i, p in enumerate(pset.particles) - if p.state in [ErrorCode.Delete]] - if len(indices) > 0 and output_file is not None: - output_file.write(pset[indices], endtime, deleted_only=True) + indices = [] + for subfield in pset.particles: + local_indices = [i for i, p in enumerate(subfield) if p.state in [ErrorCode.Delete]] + indices.append(local_indices) + # ==== do / fix ParticleFile TODO ==== # + # indices = [i for i, p in enumerate(pset.particles) if p.state in [ErrorCode.Delete]] + # if len(indices) > 0 and output_file is not None: + # output_file.write(pset[indices], endtime, deleted_only=True) pset.remove(indices) def execute(self, pset, endtime, dt, recovery=None, output_file=None, execute_once=False): """Execute this Kernel over a ParticleSet for several timesteps""" - for p in pset.particles: - p.reset_state() + for subfield in pset.particles: + for p in subfield: + p.reset_state() + + # for p in pset.particles: + # p.reset_state() if abs(dt) < 1e-6 and not execute_once: logger.warning_once("'dt' is too small, causing numerical accuracy limit problems. Please chose a higher 'dt' and rather scale the 'time' axis of the field accordingly. (related issue #762)") - # def remove_deleted(pset): - # """Utility to remove all particles that signalled deletion""" - # indices = [i for i, p in enumerate(pset.particles) - # if p.state in [ErrorCode.Delete]] - # if len(indices) > 0 and output_file is not None: - # output_file.write(pset[indices], endtime, deleted_only=True) - # pset.remove(indices) - if recovery is None: recovery = {} elif ErrorCode.ErrorOutOfBounds in recovery and ErrorCode.ErrorThroughSurface not in recovery: @@ -386,8 +399,9 @@ def execute(self, pset, endtime, dt, recovery=None, output_file=None, execute_on # Remove all particles that signalled deletion self.remove_deleted(pset, output_file=output_file, endtime=endtime) - # Identify particles that threw errors - error_particles = [p for p in pset.particles if p.state not in [ErrorCode.Success, ErrorCode.Evaluate]] + # == Identify particles that threw errors == # + # error_particles = [p for p in pset.particles if p.state not in [ErrorCode.Success, ErrorCode.Evaluate]] + error_particles = [p for subfield in pset.particles for p in subfield if p.state not in [ErrorCode.Success, ErrorCode.Evaluate]] while len(error_particles) > 0: # Apply recovery kernel @@ -419,7 +433,8 @@ def execute(self, pset, endtime, dt, recovery=None, output_file=None, execute_on else: self.execute_python(pset, endtime, dt) - error_particles = [p for p in pset.particles if p.state not in [ErrorCode.Success, ErrorCode.Evaluate]] + # error_particles = [p for p in pset.particles if p.state not in [ErrorCode.Success, ErrorCode.Evaluate]] + error_particles = [p for subfield in pset.particles for p in subfield if p.state not in [ErrorCode.Success, ErrorCode.Evaluate]] def merge(self, kernel): funcname = self.funcname + kernel.funcname From 46f97d0cb970e381028bc6f799e94d15468dd59f Mon Sep 17 00:00:00 2001 From: Christian Kehl Date: Wed, 8 Jul 2020 17:55:10 +0200 Subject: [PATCH 04/56] list_of_pset_array_trials - adapted the particle file to work with the new pset (basically: change the iterator). structure to be wrapped up and tested now. --- parcels/particlefile.py | 12 ++++++++---- 1 file changed, 8 insertions(+), 4 deletions(-) diff --git a/parcels/particlefile.py b/parcels/particlefile.py index 9485cee97c..74bac98be9 100644 --- a/parcels/particlefile.py +++ b/parcels/particlefile.py @@ -229,7 +229,8 @@ def convert_pset_to_dict(self, pset, time, deleted_only=False): if self.lasttime_written != time and \ (self.write_ondelete is False or deleted_only is True): - if pset.size == 0: + # if pset.size == 0: + if len(pset.particles) == 0 or pset.particles.shape[0] == 0: logger.warning("ParticleSet is empty on writing as array at time %g" % time) else: if deleted_only: @@ -240,13 +241,15 @@ def convert_pset_to_dict(self, pset, time, deleted_only=False): # else: # pset_towrite = [p for p in pset if time + p.dt < p.time <= time - p.dt/2 and np.isfinite(p.id)] else: - pset_towrite = [p for p in pset if time - np.abs(p.dt/2) <= p.time < time + np.abs(p.dt) and np.isfinite(p.id)] + # pset_towrite = [p for p in pset if time - np.abs(p.dt/2) <= p.time < time + np.abs(p.dt) and np.isfinite(p.id)] + pset_towrite = [p for sublist in pset.particles for p in sublist if time - np.abs(p.dt / 2) <= p.time < time + np.abs(p.dt) and np.isfinite(p.id)] if len(pset_towrite) > 0: for var in self.var_names: data_dict[var] = np.array([getattr(p, var) for p in pset_towrite]) self.maxid_written = np.max([self.maxid_written, np.max(data_dict['id'])]) - pset_errs = [p for p in pset_towrite if p.state != ErrorCode.Delete and abs(time-p.time) > 1e-3] + # pset_errs = [p for p in pset_towrite if p.state != ErrorCode.Delete and abs(time-p.time) > 1e-3] + pset_errs = [p for sublist in pset_towrite.particles for p in sublist if p.state != ErrorCode.Delete and abs(time - p.time) > 1e-3] for p in pset_errs: logger.warning_once( 'time argument in pfile.write() is %g, but a particle has time % g.' % (time, p.time)) @@ -255,7 +258,8 @@ def convert_pset_to_dict(self, pset, time, deleted_only=False): self.time_written.append(time) if len(self.var_names_once) > 0: - first_write = [p for p in pset if (p.id not in self.written_once) and _is_particle_started_yet(p, time)] + # first_write = [p for p in pset if (p.id not in self.written_once) and _is_particle_started_yet(p, time)] + first_write = [p for sublist in pset.particles for p in sublist if (p.id not in self.written_once) and _is_particle_started_yet(p, time)] data_dict_once['id'] = np.array([p.id for p in first_write]) for var in self.var_names_once: data_dict_once[var] = np.array([getattr(p, var) for p in first_write]) From 9a41211ee486c88c4b4a2f6b8eef54082ffa65ec Mon Sep 17 00:00:00 2001 From: Christian Kehl Date: Wed, 8 Jul 2020 18:24:23 +0200 Subject: [PATCH 05/56] list_of_pset_array_trials - last adaptations before testing in order to write deleted particles as intended. --- parcels/kernel.py | 34 +++++++++++++++++++++++----------- parcels/particlefile.py | 12 ++++++++++-- parcels/particleset.py | 14 +++++++------- 3 files changed, 40 insertions(+), 20 deletions(-) diff --git a/parcels/kernel.py b/parcels/kernel.py index 7635ddfd74..db2ad903c7 100644 --- a/parcels/kernel.py +++ b/parcels/kernel.py @@ -356,16 +356,29 @@ def execute_python(self, pset, endtime, dt): break def remove_deleted(self, pset, output_file, endtime): - """Utility to remove all particles that signalled deletion""" - indices = [] - for subfield in pset.particles: - local_indices = [i for i, p in enumerate(subfield) if p.state in [ErrorCode.Delete]] - indices.append(local_indices) - # ==== do / fix ParticleFile TODO ==== # - # indices = [i for i, p in enumerate(pset.particles) if p.state in [ErrorCode.Delete]] - # if len(indices) > 0 and output_file is not None: - # output_file.write(pset[indices], endtime, deleted_only=True) - pset.remove(indices) + """Utility to remove all particles that signalled deletion""" + indices = [] + dparticles = [] + for subfield in pset.particles: + # del_items = [(i, p) for i, p in enumerate(subfield) if p.state in [ErrorCode.Delete]] # we could unzip that, but in Python <3.6, more than 255 elem would make it crash + local_indices = [] + local_particles = [] + for i, p in enumerate(subfield): + if p.state == ErrorCode.Delete: + local_indices.append(i) + local_particles.append(p) + indices.append(local_indices) + if len(local_particles) > 0: + local_particles = np.array(local_particles) + else: + local_particles = None + dparticles.append(local_particles) + # ==== do / fix ParticleFile TODO ==== # + output_file.write(dparticles, endtime, deleted_only=True) # => individual entries being None -> handled in ParticleFile + # indices = [i for i, p in enumerate(pset.particles) if p.state in [ErrorCode.Delete]] + # if len(indices) > 0 and output_file is not None: + # output_file.write(pset[indices], endtime, deleted_only=True) + pset.remove(indices) def execute(self, pset, endtime, dt, recovery=None, output_file=None, execute_once=False): """Execute this Kernel over a ParticleSet for several timesteps""" @@ -411,7 +424,6 @@ def execute(self, pset, endtime, dt, recovery=None, output_file=None, execute_on if p.state == ErrorCode.Repeat: p.reset_state() elif p.state == ErrorCode.Delete: - # p.delete() pass elif p.state in recovery_map: recovery_kernel = recovery_map[p.state] diff --git a/parcels/particlefile.py b/parcels/particlefile.py index 74bac98be9..1f31454477 100644 --- a/parcels/particlefile.py +++ b/parcels/particlefile.py @@ -242,7 +242,11 @@ def convert_pset_to_dict(self, pset, time, deleted_only=False): # pset_towrite = [p for p in pset if time + p.dt < p.time <= time - p.dt/2 and np.isfinite(p.id)] else: # pset_towrite = [p for p in pset if time - np.abs(p.dt/2) <= p.time < time + np.abs(p.dt) and np.isfinite(p.id)] - pset_towrite = [p for sublist in pset.particles for p in sublist if time - np.abs(p.dt / 2) <= p.time < time + np.abs(p.dt) and np.isfinite(p.id)] + pset_towrite = [p + for sublist in pset.particles + if sublist is not None + for p in sublist + if time - np.abs(p.dt / 2) <= p.time < time + np.abs(p.dt) and np.isfinite(p.id)] if len(pset_towrite) > 0: for var in self.var_names: data_dict[var] = np.array([getattr(p, var) for p in pset_towrite]) @@ -259,7 +263,11 @@ def convert_pset_to_dict(self, pset, time, deleted_only=False): if len(self.var_names_once) > 0: # first_write = [p for p in pset if (p.id not in self.written_once) and _is_particle_started_yet(p, time)] - first_write = [p for sublist in pset.particles for p in sublist if (p.id not in self.written_once) and _is_particle_started_yet(p, time)] + first_write = [p + for sublist in pset.particles + if sublist is not None + for p in sublist + if (p.id not in self.written_once) and _is_particle_started_yet(p, time)] data_dict_once['id'] = np.array([p.id for p in first_write]) for var in self.var_names_once: data_dict_once[var] = np.array([getattr(p, var) for p in first_write]) diff --git a/parcels/particleset.py b/parcels/particleset.py index ff742444e2..6270685492 100644 --- a/parcels/particleset.py +++ b/parcels/particleset.py @@ -467,13 +467,13 @@ def lonlatdepth_dtype_from_field_interp_method(field): return np.float64 return np.float32 - # @property - # def size(self): - # nparticles = 0 - # for bracket_index in self._pid_mapping_bounds.keys(): - # nparticles += self._pid_mapping_bounds[bracket_index][2] - # return nparticles - # # return self.particles.size + @property + def size(self): + nparticles = 0 + for bracket_index in self._pid_mapping_bounds.keys(): + nparticles += self._pid_mapping_bounds[bracket_index][2] + return nparticles + # return self.particles.size def __repr__(self): return "\n".join([str(p) for sublist in self._plist for p in sublist]) From c6535162bc2150017142c0a9fa358b20ce3b2d71 Mon Sep 17 00:00:00 2001 From: Christian Kehl Date: Wed, 8 Jul 2020 19:09:01 +0200 Subject: [PATCH 06/56] list_of_pset_array_trials - fixed minor issues due to testing. Now need to impelement a general iterator to make the np.allclose() things working. --- parcels/kernel.py | 7 +++++-- parcels/particleset.py | 14 ++++++++++---- 2 files changed, 15 insertions(+), 6 deletions(-) diff --git a/parcels/kernel.py b/parcels/kernel.py index db2ad903c7..4a068f91ae 100644 --- a/parcels/kernel.py +++ b/parcels/kernel.py @@ -367,14 +367,17 @@ def remove_deleted(self, pset, output_file, endtime): if p.state == ErrorCode.Delete: local_indices.append(i) local_particles.append(p) - indices.append(local_indices) if len(local_particles) > 0: local_particles = np.array(local_particles) + local_indices = np.array(local_indices) else: local_particles = None + local_indices = None + indices.append(local_indices) dparticles.append(local_particles) # ==== do / fix ParticleFile TODO ==== # - output_file.write(dparticles, endtime, deleted_only=True) # => individual entries being None -> handled in ParticleFile + if output_file is not None: + output_file.write(dparticles, endtime, deleted_only=True) # => individual entries being None -> handled in ParticleFile # indices = [i for i, p in enumerate(pset.particles) if p.state in [ErrorCode.Delete]] # if len(indices) > 0 and output_file is not None: # output_file.write(pset[indices], endtime, deleted_only=True) diff --git a/parcels/particleset.py b/parcels/particleset.py index 6270685492..d916140852 100644 --- a/parcels/particleset.py +++ b/parcels/particleset.py @@ -160,13 +160,14 @@ def convert_to_array(var): assert lon.shape[0] == lat.shape[0], ('Length of lon and lat do not match.') self._plist = [] start_index = 0 - end_index = min(self.nlist_limit , lon.shape[0]) + end_index = 0 #min(self.nlist_limit , lon.shape[0]) while end_index < lon.shape[0]: + end_index = min(lon.shape[0], start_index + self.nlist_limit) self._plist.append(np.empty(end_index-start_index, dtype=pclass)) bracket_index = len(self._plist)-1 self._pid_mapping_bounds[bracket_index] = (np.iinfo(np.int32).max, np.iinfo(np.int32).min, end_index-start_index) start_index = end_index - end_index = min(lon.shape[0], end_index+self.nlist_limit) + # end_index = min(lon.shape[0], end_index+self.nlist_limit) # self._particle_data = None self._plist_c = None self._pclass = pclass @@ -532,7 +533,7 @@ def add(self, particles): def _merge_brackets_(self): lw_bound_nlist = int(self.nlist_limit / 2) nmerges = 1 - while nmerges >= 0: + while nmerges > 0: nmerges = 0 trg_bracket = None trg_bracket_index = -1 @@ -559,6 +560,7 @@ def _merge_brackets_(self): self._plist.remove(self._plist[src_bracket_index]) self._plist_c.remove(self._plist_c[src_bracket_index]) self._pid_mapping_bounds.pop(src_bracket_index) + nmerges += 1 def remove(self, indices): """Method to remove particles from the ParticleSet, based on their `indices`""" @@ -566,6 +568,8 @@ def remove(self, indices): len(self._plist), len(indices))) for bracket_index in range(len(indices)): local_indices = indices[bracket_index] + if local_indices is None: + continue if isinstance(local_indices, collections.Iterable): local_indices = np.array(local_indices) self._remove_(bracket_index, local_indices) @@ -688,7 +692,8 @@ def execute(self, pyfunc=AdvectionRK4, endtime=None, runtime=None, dt=1., # Derive _starttime and endtime from arguments or fieldset defaults if runtime is not None and endtime is not None: raise RuntimeError('Only one of (endtime, runtime) can be specified') - _starttime = min([p.time for p in self]) if dt >= 0 else max([p.time for p in self]) + # _starttime = min([p.time for p in self]) if dt >= 0 else max([p.time for p in self]) + _starttime = min([p.time for sublist in self._plist for p in sublist]) if dt >= 0 else max([p.time for sublist in self._plist for p in sublist]) if self.repeatdt is not None and self.repeat_starttime is None: self.repeat_starttime = _starttime if runtime is not None: @@ -754,6 +759,7 @@ def execute(self, pyfunc=AdvectionRK4, endtime=None, runtime=None, dt=1., else: time = max(next_prelease, next_input, next_output, next_movie, next_callback, endtime) self.kernel.execute(self, endtime=time, dt=dt, recovery=recovery, output_file=output_file, execute_once=execute_once) + logger.info("Computed time = {}".format(time)) if abs(time-next_prelease) < tol: pset_new = ParticleSet(fieldset=self.fieldset, time=time, lon=self.repeatlon, lat=self.repeatlat, depth=self.repeatdepth, From 9997a98c62036fdef45bb5f5171fa9f2332b9b29 Mon Sep 17 00:00:00 2001 From: Christian Kehl Date: Wed, 8 Jul 2020 20:18:11 +0200 Subject: [PATCH 07/56] list_of_pset_array_trials - transformed PID mapping from dict into list --- parcels/particleset.py | 14 ++++++++------ 1 file changed, 8 insertions(+), 6 deletions(-) diff --git a/parcels/particleset.py b/parcels/particleset.py index d916140852..10febeedef 100644 --- a/parcels/particleset.py +++ b/parcels/particleset.py @@ -154,7 +154,9 @@ def convert_to_array(var): 'lon lat depth precision should be set to either np.float32 or np.float64' JITParticle.set_lonlatdepth_dtype(self.lonlatdepth_dtype) - self._pid_mapping_bounds = {} # plist bracket_index -> (min_id, max_id, size_bracket) + # == this should be a list not a dict! == # + self._pid_mapping_bounds = [] # plist bracket_index -> (min_id, max_id, size_bracket) + # ======================================= # self.nlist_limit = 4096 # self.particles = np.empty(lon.size, dtype=pclass) assert lon.shape[0] == lat.shape[0], ('Length of lon and lat do not match.') @@ -164,8 +166,8 @@ def convert_to_array(var): while end_index < lon.shape[0]: end_index = min(lon.shape[0], start_index + self.nlist_limit) self._plist.append(np.empty(end_index-start_index, dtype=pclass)) + self._pid_mapping_bounds.append((np.iinfo(np.int32).max, np.iinfo(np.int32).min, end_index-start_index)) bracket_index = len(self._plist)-1 - self._pid_mapping_bounds[bracket_index] = (np.iinfo(np.int32).max, np.iinfo(np.int32).min, end_index-start_index) start_index = end_index # end_index = min(lon.shape[0], end_index+self.nlist_limit) # self._particle_data = None @@ -304,7 +306,7 @@ def covert_static_array_to_list_array(self, property_array): if not isinstance(property_array, np.ndarray): property_array = np.array(property_array) nparticles = 0 - for bracket_index in self._pid_mapping_bounds.keys(): + for bracket_index in range(len(self._pid_mapping_bounds)): nparticles += self._pid_mapping_bounds[bracket_index][2] assert property_array.shape[0] == nparticles, ("Attaching the property array to the list of particle prohibited because num. of entries do not match.") results = [] @@ -471,7 +473,7 @@ def lonlatdepth_dtype_from_field_interp_method(field): @property def size(self): nparticles = 0 - for bracket_index in self._pid_mapping_bounds.keys(): + for bracket_index in range(len(self._pid_mapping_bounds)): nparticles += self._pid_mapping_bounds[bracket_index][2] return nparticles # return self.particles.size @@ -510,7 +512,7 @@ def add(self, particles): offset = len(self._plist) self._plist += particles.particles # self.pid_mapping_bounds.update(particles.pid_mapping_bounds) - for i in particles.pid_mapping_bounds.keys(): + for i in range(len(particles.pid_mapping_bounds)): self._pid_mapping_bounds[offset+i] = particles.pid_mapping_bounds[i] else: raise NotImplementedError('Only ParticleSets can be added to a ParticleSet') @@ -559,7 +561,7 @@ def _merge_brackets_(self): p._cptr = pdata self._plist.remove(self._plist[src_bracket_index]) self._plist_c.remove(self._plist_c[src_bracket_index]) - self._pid_mapping_bounds.pop(src_bracket_index) + self._pid_mapping_bounds.remove(src_bracket_index) nmerges += 1 def remove(self, indices): From 73ebfacc633c0bcb3a26ddaa613d2b28a06d548e Mon Sep 17 00:00:00 2001 From: Christian Kehl Date: Wed, 8 Jul 2020 20:26:39 +0200 Subject: [PATCH 08/56] list_of_pset_array_trials - atted __getitem__(self, key) func to ParticleSet for element iteration - slow, but working. --- parcels/particleset.py | 25 +++++++++++++++++-------- 1 file changed, 17 insertions(+), 8 deletions(-) diff --git a/parcels/particleset.py b/parcels/particleset.py index 10febeedef..a6d460703a 100644 --- a/parcels/particleset.py +++ b/parcels/particleset.py @@ -259,9 +259,9 @@ def get_list_array_index(self, pdata): start_index = 0 end_index = self._pid_mapping_bounds[bracket_index][2] while start_index Date: Thu, 9 Jul 2020 13:10:29 +0200 Subject: [PATCH 09/56] list_of_pset_array_trials - fixing address translation 1-array index -> LoA index --- parcels/particleset.py | 25 ++++++++++++++++++------- 1 file changed, 18 insertions(+), 7 deletions(-) diff --git a/parcels/particleset.py b/parcels/particleset.py index a6d460703a..6fafea011f 100644 --- a/parcels/particleset.py +++ b/parcels/particleset.py @@ -258,7 +258,7 @@ def get_list_array_index(self, pdata): bracket_index = 0 start_index = 0 end_index = self._pid_mapping_bounds[bracket_index][2] - while start_index= end_index: start_index = end_index end_index = start_index+self._pid_mapping_bounds[bracket_index][2] bracket_index += 1 @@ -313,7 +313,7 @@ def covert_static_array_to_list_array(self, property_array): bracket_index = 0 start_index = 0 end_index = self._pid_mapping_bounds[bracket_index][2] - while start_index= key: + while key >= end_index: start_index = end_index end_index = start_index+self._pid_mapping_bounds[bracket_index][2] bracket_index += 1 @@ -610,15 +612,24 @@ def _remove_(self, bracket_index, local_indices): local_indices = np.array(local_indices) return_p = self._plist[bracket_index][local_indices] self._plist[bracket_index] = np.delete(self._plist[bracket_index], local_indices) - local_ids = np.array([p.id for p in self._plist[bracket_index]]) - self.pid_mapping_bounds[bracket_index] = (np.min(local_ids), np.max(local_ids), local_ids.shape[0]) if self._ptype.uses_jit: self._plist_c[bracket_index] = np.delete(self._plist_c[bracket_index], local_indices) + if self._plist[bracket_index] is not None and self._plist[bracket_index].shape[0] > 0: + local_ids = np.array([p.id for p in self._plist[bracket_index]]) + self.pid_mapping_bounds[bracket_index] = (np.min(local_ids), np.max(local_ids), local_ids.shape[0]) + if self._ptype.uses_jit: # # == Update C-pointer on particles == # # for p, pdata in zip(self._plist[bracket_index], self._plist_c[bracket_index]): # p._cptr = pdata - for slot_index in range(len(self._plist[bracket_index])): - self._plist[bracket_index][slot_index].update_cptr(self._plist_c[bracket_index][slot_index]) + for slot_index in range(len(self._plist[bracket_index])): + self._plist[bracket_index][slot_index].update_cptr(self._plist_c[bracket_index][slot_index]) + else: + if len(self._plist) > 1: + self._plist.remove(bracket_index) + self._plist_c.remove(bracket_index) + self._pid_mapping_bounds.remove(bracket_index) + else: + self.pid_mapping_bounds[bracket_index] = (np.iinfo(np.int32).max, np.iinfo(np.int32).min, 0) # return particles return return_p From 518753efc705a93009f00db2d32788ed7710992f Mon Sep 17 00:00:00 2001 From: Christian Kehl Date: Tue, 14 Jul 2020 17:38:01 +0200 Subject: [PATCH 10/56] list_of_pset_array_trials - fixed errors emerging from tests --- parcels/kernel.py | 14 ++-- parcels/particle.py | 2 +- parcels/particlefile.py | 22 ++++-- parcels/particleset.py | 144 +++++++++++++++++++++++++++++++----- tests/test_particle_sets.py | 3 +- 5 files changed, 149 insertions(+), 36 deletions(-) diff --git a/parcels/kernel.py b/parcels/kernel.py index 4a068f91ae..e98e122540 100644 --- a/parcels/kernel.py +++ b/parcels/kernel.py @@ -219,8 +219,8 @@ def load_lib(self): def execute_jit(self, pset, endtime, dt): """Invokes JIT engine to perform the core update loop""" - if len(pset.particles) > 0 and pset.particles[0].shape[0]: - assert pset.fieldset.gridset.size == len(pset.particles[0][0].xi), \ + if len(pset.particles_a) > 0 and pset.particles_a[0].shape[0]: + assert pset.fieldset.gridset.size == len(pset.particles_a[0][0].xi), \ 'FieldSet has different amount of grids than Particle.xi. Have you added Fields after creating the ParticleSet?' # if len(pset.particles) > 0: @@ -280,7 +280,7 @@ def execute_python(self, pset, endtime, dt): f.data = np.array(f.data) # for p in pset.particles: - for subfield in pset.particles: + for subfield in pset.particles_a: for p in subfield: ptype = p.getPType() # Don't execute particles that aren't started yet @@ -359,7 +359,7 @@ def remove_deleted(self, pset, output_file, endtime): """Utility to remove all particles that signalled deletion""" indices = [] dparticles = [] - for subfield in pset.particles: + for subfield in pset.particles_a: # del_items = [(i, p) for i, p in enumerate(subfield) if p.state in [ErrorCode.Delete]] # we could unzip that, but in Python <3.6, more than 255 elem would make it crash local_indices = [] local_particles = [] @@ -385,7 +385,7 @@ def remove_deleted(self, pset, output_file, endtime): def execute(self, pset, endtime, dt, recovery=None, output_file=None, execute_once=False): """Execute this Kernel over a ParticleSet for several timesteps""" - for subfield in pset.particles: + for subfield in pset.particles_a: for p in subfield: p.reset_state() @@ -417,7 +417,7 @@ def execute(self, pset, endtime, dt, recovery=None, output_file=None, execute_on # == Identify particles that threw errors == # # error_particles = [p for p in pset.particles if p.state not in [ErrorCode.Success, ErrorCode.Evaluate]] - error_particles = [p for subfield in pset.particles for p in subfield if p.state not in [ErrorCode.Success, ErrorCode.Evaluate]] + error_particles = [p for subfield in pset.particles_a for p in subfield if p.state not in [ErrorCode.Success, ErrorCode.Evaluate]] while len(error_particles) > 0: # Apply recovery kernel @@ -449,7 +449,7 @@ def execute(self, pset, endtime, dt, recovery=None, output_file=None, execute_on self.execute_python(pset, endtime, dt) # error_particles = [p for p in pset.particles if p.state not in [ErrorCode.Success, ErrorCode.Evaluate]] - error_particles = [p for subfield in pset.particles for p in subfield if p.state not in [ErrorCode.Success, ErrorCode.Evaluate]] + error_particles = [p for subfield in pset.particles_a for p in subfield if p.state not in [ErrorCode.Success, ErrorCode.Evaluate]] def merge(self, kernel): funcname = self.funcname + kernel.funcname diff --git a/parcels/particle.py b/parcels/particle.py index 1800b52b77..0762b46749 100644 --- a/parcels/particle.py +++ b/parcels/particle.py @@ -260,5 +260,5 @@ def __init__(self, *args, **kwargs): setattr(self, 'c'+index, getattr(self, index+'p').value) def update_cptr(self, new_ptr): - np.copyto(new_ptr, self._cptr) + # np.copyto(new_ptr, self._cptr) self._cptr = new_ptr diff --git a/parcels/particlefile.py b/parcels/particlefile.py index 1f31454477..1b369a33bf 100644 --- a/parcels/particlefile.py +++ b/parcels/particlefile.py @@ -229,12 +229,19 @@ def convert_pset_to_dict(self, pset, time, deleted_only=False): if self.lasttime_written != time and \ (self.write_ondelete is False or deleted_only is True): + if not isinstance(pset, list): + pset = pset.particles_a + # if pset.size == 0: - if len(pset.particles) == 0 or pset.particles.shape[0] == 0: + if len(pset) == 0: # or pset.shape[0] == 0: logger.warning("ParticleSet is empty on writing as array at time %g" % time) else: if deleted_only: - pset_towrite = pset + pset_towrite = [p + for sublist in pset + if sublist is not None + for p in sublist] + # pset_towrite = pset # == commented due to git rebase to master 27 02 2020 == # # elif pset[0].dt > 0: # pset_towrite = [p for p in pset if time - p.dt/2 <= p.time < time + p.dt and np.isfinite(p.id)] @@ -243,7 +250,7 @@ def convert_pset_to_dict(self, pset, time, deleted_only=False): else: # pset_towrite = [p for p in pset if time - np.abs(p.dt/2) <= p.time < time + np.abs(p.dt) and np.isfinite(p.id)] pset_towrite = [p - for sublist in pset.particles + for sublist in pset if sublist is not None for p in sublist if time - np.abs(p.dt / 2) <= p.time < time + np.abs(p.dt) and np.isfinite(p.id)] @@ -252,8 +259,8 @@ def convert_pset_to_dict(self, pset, time, deleted_only=False): data_dict[var] = np.array([getattr(p, var) for p in pset_towrite]) self.maxid_written = np.max([self.maxid_written, np.max(data_dict['id'])]) - # pset_errs = [p for p in pset_towrite if p.state != ErrorCode.Delete and abs(time-p.time) > 1e-3] - pset_errs = [p for sublist in pset_towrite.particles for p in sublist if p.state != ErrorCode.Delete and abs(time - p.time) > 1e-3] + pset_errs = [p for p in pset_towrite if p.state != ErrorCode.Delete and abs(time-p.time) > 1e-3] + # pset_errs = [p for sublist in pset_towrite.particles_a for p in sublist if p.state != ErrorCode.Delete and abs(time - p.time) > 1e-3] for p in pset_errs: logger.warning_once( 'time argument in pfile.write() is %g, but a particle has time % g.' % (time, p.time)) @@ -264,7 +271,7 @@ def convert_pset_to_dict(self, pset, time, deleted_only=False): if len(self.var_names_once) > 0: # first_write = [p for p in pset if (p.id not in self.written_once) and _is_particle_started_yet(p, time)] first_write = [p - for sublist in pset.particles + for sublist in pset if sublist is not None for p in sublist if (p.id not in self.written_once) and _is_particle_started_yet(p, time)] @@ -291,7 +298,8 @@ def dump_dict_to_npy(self, data_dict, data_dict_once): self.file_list.append(tmpfilename) if len(data_dict_once) > 0: - tmpfilename = os.path.join(self.tempwritedir, str(len(self.file_list)) + '_once.npy') + # tmpfilename = os.path.join(self.tempwritedir, str(len(self.file_list)) + '_once.npy') + tmpfilename = os.path.join(self.tempwritedir, str(len(self.file_list_once)) + '_once.npy') with open(tmpfilename, 'wb') as f: np.save(f, data_dict_once) self.file_list_once.append(tmpfilename) diff --git a/parcels/particleset.py b/parcels/particleset.py index 6fafea011f..03d47a08d3 100644 --- a/parcels/particleset.py +++ b/parcels/particleset.py @@ -211,6 +211,17 @@ def convert_to_array(var): @property def particles(self): + if len(self._plist) == 0: + return None + result = self._plist[0] + if len(self._plist) > 1: + for bracket_index in range(0, len(self._plist)): + result = np.concatenate((result, self._plist[bracket_index]), axis=0) + return result + # return self._plist + + @property + def particles_a(self): return self._plist @property @@ -320,6 +331,22 @@ def covert_static_array_to_list_array(self, property_array): bracket_index += 1 return results + @staticmethod + def convert_list_of_arrays_to_static_list(pset): + return [p + for sublist in pset + if sublist is not None + for p in sublist] + + def convert_pset_to_static_list(self): + return [p + for sublist in self.particles_a + if sublist is not None + for p in sublist] + + def convert_pset_to_static_array(self): + return np.array(self.convert_pset_to_static_list(), dtype=self._pclass) + @classmethod def from_list(cls, fieldset, pclass, lon, lat, depth=None, time=None, repeatdt=None, lonlatdepth_dtype=None, **kwargs): """Initialise the ParticleSet from lists of lon and lat @@ -499,6 +526,9 @@ def __getitem__(self, key): return self._plist[bracket_index][slot_index] # return self.particles[key] + def __setitem__(self, key, value): + assert isinstance(key, tuple) and len(key) == 2, ("Error: trying to set/address data item without (bracket_index, slot_index) key.") + self._plist[key[0]][key[1]] = value # def __setitem__(self, key, value): # self.particles[key] = value @@ -520,17 +550,15 @@ def _create_progressbar_(self, starttime, endtime): def add(self, particles): """Method to add particles to the ParticleSet""" if isinstance(particles, ParticleSet): - offset = len(self._plist) - self._plist += particles.particles - # self.pid_mapping_bounds.update(particles.pid_mapping_bounds) - for i in range(len(particles.pid_mapping_bounds)): - self._pid_mapping_bounds[offset+i] = particles.pid_mapping_bounds[i] + self._plist += particles.particles_a + self._pid_mapping_bounds += particles.pid_mapping_bounds + + if self.ptype.uses_jit: + for sublist in particles.particles_c: + self._plist_c.append(sublist) else: raise NotImplementedError('Only ParticleSets can be added to a ParticleSet') - for sublist in particles.particles_c: - self._plist_c.append(sublist) - # if isinstance(particles, ParticleSet): # particles = particles.particles # else: @@ -570,17 +598,88 @@ def _merge_brackets_(self): self._plist_c[trg_bracket_index] = np.append(self._plist_c[trg_bracket_index], self._plist_c[src_bracket_index]) for p, pdata in zip(self._plist[trg_bracket_index], self._plist_c[trg_bracket_index]): p._cptr = pdata - self._plist.remove(self._plist[src_bracket_index]) - self._plist_c.remove(self._plist_c[src_bracket_index]) - self._pid_mapping_bounds.remove(src_bracket_index) + # self._plist.remove(self._plist[src_bracket_index]) + del self._plist[src_bracket_index] + if self.ptype.uses_jit: + # self._plist_c.remove(self._plist_c[src_bracket_index]) + del self._plist_c[src_bracket_index] + # self._pid_mapping_bounds.remove(self._pid_mapping_bounds[src_bracket_index]) + del self._pid_mapping_bounds[src_bracket_index] nmerges += 1 + def pop(self, indices=-1): + """ + Method to remove particles from the ParticleSet, based on their `indices`. + Returns numpy.ndarray of removed particles. + """ + lindices = indices + if isinstance(indices, np.ndarray) or not isinstance(indices, list): + ncounter = 0 + for bracket in self._pid_mapping_bounds: + ncounter += bracket[2] + if isinstance(indices, np.ndarray): + indices = indices.tolist() + if not isinstance(indices, list): + indices = [indices, ] + lindices = [] + for i in range(len(self._plist)): + lindices.append(None) + for gindex in indices: + while gindex < 0: + gindex = ncounter+gindex + (bracket_index, slot_index) = self.get_list_array_index(gindex) + if lindices[bracket_index] is None: + lindices[bracket_index] = [slot_index, ] + else: + lindices[bracket_index].append(slot_index) + assert len(lindices) == len(self._plist), ("ParticleSet.remove() - particle set length ({}) does not match index list length ({})".format( + len(self._plist), len(lindices))) + result = None + for bracket_index in range(len(lindices)): + local_indices = lindices[bracket_index] + if local_indices is None: + continue + if isinstance(local_indices, collections.Iterable): + local_indices = np.array(local_indices) + if result is None: + result = self._remove_(bracket_index, local_indices) + else: + result = np.concatenate((result, self._remove_(bracket_index, local_indices)), axis=0) + self._merge_brackets_() + if result is None: + return result + result = result.tolist() + if len(result) == 1: + return result[0] + return result + def remove(self, indices): """Method to remove particles from the ParticleSet, based on their `indices`""" - assert len(indices) == len(self._plist), ("ParticleSet.remove() - particle set length ({}) does not match index list length ({})".format( - len(self._plist), len(indices))) - for bracket_index in range(len(indices)): - local_indices = indices[bracket_index] + lindices = indices + if isinstance(indices, np.ndarray) or not isinstance(indices, list): + ncounter = 0 + for bracket in self._pid_mapping_bounds: + ncounter += bracket[2] + if isinstance(indices, np.ndarray): + indices = indices.tolist() + if not isinstance(indices, list): + indices = [indices, ] + lindices = [] + for i in range(len(self._plist)): + lindices.append(None) + for gindex in indices: + while gindex < 0: + gindex = ncounter+gindex + (bracket_index, slot_index) = self.get_list_array_index(gindex) + if lindices[bracket_index] is None: + lindices[bracket_index] = [slot_index, ] + else: + lindices[bracket_index].append(slot_index) + assert len(lindices) == len(self._plist), ("ParticleSet.remove() - particle set length ({}) does not match index list length ({})".format( + len(self._plist), len(lindices))) + + for bracket_index in range(len(lindices)): + local_indices = lindices[bracket_index] if local_indices is None: continue if isinstance(local_indices, collections.Iterable): @@ -625,9 +724,10 @@ def _remove_(self, bracket_index, local_indices): self._plist[bracket_index][slot_index].update_cptr(self._plist_c[bracket_index][slot_index]) else: if len(self._plist) > 1: - self._plist.remove(bracket_index) - self._plist_c.remove(bracket_index) - self._pid_mapping_bounds.remove(bracket_index) + self._plist.remove(self._plist[bracket_index]) + if self._ptype.uses_jit: + self._plist_c.remove(self._plist_c[bracket_index]) + self._pid_mapping_bounds.remove(self._pid_mapping_bounds[bracket_index]) else: self.pid_mapping_bounds[bracket_index] = (np.iinfo(np.int32).max, np.iinfo(np.int32).min, 0) # return particles @@ -787,7 +887,7 @@ def execute(self, pyfunc=AdvectionRK4, endtime=None, runtime=None, dt=1., lat=self.repeatlat, depth=self.repeatdepth, pclass=self.repeatpclass, lonlatdepth_dtype=self.lonlatdepth_dtype, partitions=False, pid_orig=self.repeatpid, **self.repeatkwargs) - for sublist in pset_new.particles: + for sublist in pset_new.particles_a: for p in sublist: p.dt = dt self.add(pset_new) @@ -882,7 +982,11 @@ def density(self, field=None, particle_val=None, relative=False, area_scale=Fals # density[yi, xi] += particle_val[pi] if relative: - density /= np.sum(particle_val) + psum = 0 + for bracket_index, subfield in enumerate(particle_val): + psum += np.sum(particle_val[bracket_index]) + density /= psum + # density /= np.sum(particle_val) if area_scale: density /= field.cell_areas() diff --git a/tests/test_particle_sets.py b/tests/test_particle_sets.py index bda1e55398..84fada24df 100644 --- a/tests/test_particle_sets.py +++ b/tests/test_particle_sets.py @@ -306,7 +306,8 @@ def test_pset_remove_index(fieldset, mode, npart=100): lat = np.linspace(1, 0, npart) pset = ParticleSet(fieldset, lon=lon, lat=lat, pclass=ptype[mode], lonlatdepth_dtype=np.float64) for ilon, ilat in zip(lon[::-1], lat[::-1]): - p = pset.remove(-1) + # p = pset.remove(-1) + p = pset.pop() assert(p.lon == ilon) assert(p.lat == ilat) assert(pset.size == 0) From 49f505577ea2d38311391c04fb80ed89aa0c66fb Mon Sep 17 00:00:00 2001 From: Christian Kehl Date: Tue, 14 Jul 2020 18:20:29 +0200 Subject: [PATCH 11/56] list_of_pset_array_trials - fixing linter errors --- parcels/particlefile.py | 2 +- parcels/particleset.py | 23 +++++----- parcels/particleset_benchmark.py | 77 ++++++++++++++++++++++++++++++++ tests/test_fieldset_sampling.py | 52 ++++++--------------- 4 files changed, 102 insertions(+), 52 deletions(-) diff --git a/parcels/particlefile.py b/parcels/particlefile.py index 1b369a33bf..ea4895b8a1 100644 --- a/parcels/particlefile.py +++ b/parcels/particlefile.py @@ -233,7 +233,7 @@ def convert_pset_to_dict(self, pset, time, deleted_only=False): pset = pset.particles_a # if pset.size == 0: - if len(pset) == 0: # or pset.shape[0] == 0: + if len(pset) == 0: logger.warning("ParticleSet is empty on writing as array at time %g" % time) else: if deleted_only: diff --git a/parcels/particleset.py b/parcels/particleset.py index 03d47a08d3..bba5488c6d 100644 --- a/parcels/particleset.py +++ b/parcels/particleset.py @@ -16,7 +16,6 @@ from parcels.kernels.advection import AdvectionRK4 from parcels.particle import JITParticle from parcels.particlefile import ParticleFile -from parcels.tools.error import ErrorCode from parcels.tools.loggers import logger try: from mpi4py import MPI @@ -31,6 +30,7 @@ __all__ = ['ParticleSet'] + class ParticleSet(object): """Container class for storing particle and executing kernel over them. @@ -155,14 +155,14 @@ def convert_to_array(var): JITParticle.set_lonlatdepth_dtype(self.lonlatdepth_dtype) # == this should be a list not a dict! == # - self._pid_mapping_bounds = [] # plist bracket_index -> (min_id, max_id, size_bracket) + self._pid_mapping_bounds = [] # plist bracket_index -> (min_id, max_id, size_bracket) # ======================================= # self.nlist_limit = 4096 # self.particles = np.empty(lon.size, dtype=pclass) assert lon.shape[0] == lat.shape[0], ('Length of lon and lat do not match.') self._plist = [] start_index = 0 - end_index = 0 #min(self.nlist_limit , lon.shape[0]) + end_index = 0 while end_index < lon.shape[0]: end_index = min(lon.shape[0], start_index + self.nlist_limit) self._plist.append(np.empty(end_index-start_index, dtype=pclass)) @@ -189,7 +189,7 @@ def convert_to_array(var): if lon is not None and lat is not None: # == Initialise from arrays of lon/lat coordinates == # - # assert self.particles.size == lon.size and self.particles.size == lat.size, ('Size of ParticleSet does not match length of lon and lat.') + # assert self.particles.size == lon.size and self.particles.size == lat.size, ('Size of ParticleSet does not match length of lon and lat.') for i in range(lon.size): bracket_index = int(float(i)/self.nlist_limit) @@ -197,7 +197,7 @@ def convert_to_array(var): self._plist[bracket_index][slot_index] = pclass(lon[i], lat[i], pid[i], fieldset=fieldset, depth=depth[i], cptr=self.cptr(bracket_index, slot_index), time=time[i]) bracket_info = self._pid_mapping_bounds[bracket_index] self._pid_mapping_bounds[bracket_index] = (min(pid[i], bracket_info[0]), max(pid[i], bracket_info[1]), bracket_info[2]) - # self.particles[i] = pclass(lon[i], lat[i], pid[i], fieldset=fieldset, depth=depth[i], cptr=cptr(i), time=time[i]) + # self.particles[i] = pclass(lon[i], lat[i], pid[i], fieldset=fieldset, depth=depth[i], cptr=cptr(i), time=time[i]) # == Set other Variables if provided == # for kwvar in kwargs: if not hasattr(self._plist[bracket_index][slot_index], kwvar): @@ -277,11 +277,9 @@ def get_list_array_index(self, pdata): return (bracket_index, slot_index) elif isinstance(pdata, self.pclass): # == parameter is a particle itself == # - # nparticle = 0 bracket_index = 0 bracket_info = self.pid_mapping_bounds[bracket_index] while bracket_index < len(self._plist) and ((pdata.id < bracket_info[0]) or (pdata.id > bracket_info[1])): - # nparticle += bracket_info[2] bracket_index += 1 bracket_info = self.pid_mapping_bounds[bracket_index] if bracket_index >= len(self._plist): @@ -290,11 +288,9 @@ def get_list_array_index(self, pdata): slot_index = 0 while (slot_index < self._plist[bracket_index].shape[0]) and (self._plist[bracket_index][slot_index].id != pdata.id): slot_index += 1 - # nparticle += 1 if slot_index >= self._plist[bracket_index].shape[0]: logger.warn_once("ParticleSet.get_cptr_index() - requested particle ({}) not found.".format(pdata)) return -1 - # return nparticle return (bracket_index, slot_index) def get_list_array_index_by_PID(self, pdata): @@ -717,9 +713,11 @@ def _remove_(self, bracket_index, local_indices): local_ids = np.array([p.id for p in self._plist[bracket_index]]) self.pid_mapping_bounds[bracket_index] = (np.min(local_ids), np.max(local_ids), local_ids.shape[0]) if self._ptype.uses_jit: - # # == Update C-pointer on particles == # - # for p, pdata in zip(self._plist[bracket_index], self._plist_c[bracket_index]): - # p._cptr = pdata + # == Update C-pointer on particles == # + # -- OLD -- # + # for p, pdata in zip(self._plist[bracket_index], self._plist_c[bracket_index]): + # p._cptr = pdata + # -- END -- # for slot_index in range(len(self._plist[bracket_index])): self._plist[bracket_index][slot_index].update_cptr(self._plist_c[bracket_index][slot_index]) else: @@ -970,7 +968,6 @@ def density(self, field=None, particle_val=None, relative=False, area_scale=Fals _, _, _, xi, yi, _ = field.search_indices(p.lon, p.lat, p.depth, 0, 0, search2D=True) density[yi, xi] += particle_val[bracket_index][slot_index] - # for pi, p in enumerate(self.particles): # try: # breaks if either p.xi, p.yi, p.zi, p.ti do not exist (in scipy) or field not in fieldset # if p.ti[field.igrid] < 0: # xi, yi, zi, ti, not initialised diff --git a/parcels/particleset_benchmark.py b/parcels/particleset_benchmark.py index 5bc903acbf..69ed754d0c 100644 --- a/parcels/particleset_benchmark.py +++ b/parcels/particleset_benchmark.py @@ -51,6 +51,82 @@ def measure_mem_usage(): USE_ASYNC_MEMLOG = False USE_RUSE_SYNC_MEMLOG = False # can be faulty +class ParticleSet_TimingLog(): + stime = 0 + etime = 0 + mtime = 0 + samples = [] + times_steps = [] + _iter = 0 + + def start_timing(self): + if MPI: + mpi_comm = MPI.COMM_WORLD + mpi_rank = mpi_comm.Get_rank() + if mpi_rank == 0: + # self.stime = MPI.Wtime() + # self.stime = time_module.perf_counter() + self.stime = time_module.process_time() + else: + self.stime = time_module.perf_counter() + + def stop_timing(self): + if MPI: + mpi_comm = MPI.COMM_WORLD + mpi_rank = mpi_comm.Get_rank() + if mpi_rank == 0: + # self.etime = MPI.Wtime() + # self.etime = time_module.perf_counter() + self.etime = time_module.process_time() + else: + self.etime = time_module.perf_counter() + + def accumulate_timing(self): + if MPI: + mpi_comm = MPI.COMM_WORLD + mpi_rank = mpi_comm.Get_rank() + if mpi_rank == 0: + self.mtime += (self.etime-self.stime) + else: + self.mtime = 0 + else: + self.mtime += (self.etime-self.stime) + + def advance_iteration(self): + if MPI: + mpi_comm = MPI.COMM_WORLD + mpi_rank = mpi_comm.Get_rank() + if mpi_rank == 0: + self.times_steps.append(self.mtime) + self.samples.append(self._iter) + self._iter += 1 + self.mtime = 0 + else: + self.times_steps.append(self.mtime) + self.samples.append(self._iter) + self._iter += 1 + self.mtime = 0 + + +class ParticleSet_ParamLogging(): + samples = [] + params = [] + _iter = 0 + + def advance_iteration(self, param): + if MPI: + mpi_comm = MPI.COMM_WORLD + mpi_rank = mpi_comm.Get_rank() + if mpi_rank == 0: + self.params.append(param) + self.samples.append(self._iter) + self._iter += 1 + else: + self.params.append(param) + self.samples.append(self._iter) + self._iter += 1 + + class ParticleSet_Benchmark(ParticleSet): def __init__(self, fieldset, pclass=JITParticle, lon=None, lat=None, depth=None, time=None, repeatdt=None, @@ -481,3 +557,4 @@ def plot_and_log(self, total_times = None, compute_times = None, io_times = None max_mem = max(max_mem_sync, max_mem_async) data_string += "{:10.4f}, {:10.4f}, {:10.4f}, {:10.4f}, {}".format(total_times.sum(), compute_times.sum(), io_times.sum(), plot_times.sum(), max_mem) f.write(data_string) + diff --git a/tests/test_fieldset_sampling.py b/tests/test_fieldset_sampling.py index c9abfc5d63..8fdb9ee9aa 100644 --- a/tests/test_fieldset_sampling.py +++ b/tests/test_fieldset_sampling.py @@ -451,13 +451,13 @@ def test_sample(particle, fieldset, time): if mode == 'jit': assert np.all(np.array([len(p.xi) == fieldset.gridset.size for p in pset])) assert np.all(np.array([[p.xi[i] >= 0 for i in range(0, len(p.xi))] for p in pset])) - #assert np.all(pset.xi[:, fieldset.B.igrid] < xdim * 4) - #assert np.all(pset.xi[:, 0] < xdim) - #assert pset.yi.shape[0] == len(pset.lon) - #assert pset.yi.shape[1] == fieldset.gridset.size - #assert np.all(pset.yi >= 0) - #assert np.all(pset.yi[:, fieldset.B.igrid] < ydim * 3) - #assert np.all(pset.yi[:, 0] < ydim) + # assert np.all(pset.xi[:, fieldset.B.igrid] < xdim * 4) + # assert np.all(pset.xi[:, 0] < xdim) + # assert pset.yi.shape[0] == len(pset.lon) + # assert pset.yi.shape[1] == fieldset.gridset.size + # assert np.all(pset.yi >= 0) + # assert np.all(pset.yi[:, fieldset.B.igrid] < ydim * 3) + # assert np.all(pset.yi[:, 0] < ydim) @pytest.mark.parametrize('mode', ['jit', 'scipy']) @@ -493,13 +493,13 @@ def test_sample(particle, fieldset, time): if mode == 'jit': assert np.all(np.array([len(p.xi) == fieldset.gridset.size for p in pset])) assert np.all(np.array([[p.xi[i] >= 0 for i in range(0, len(p.xi))] for p in pset])) - #assert np.all(pset.xi[:, fieldset.B.igrid] < xdim * 4) - #assert np.all(pset.xi[:, 0] < xdim) - #assert pset.yi.shape[0] == len(pset.lon) - #assert pset.yi.shape[1] == fieldset.gridset.size - #assert np.all(pset.yi >= 0) - #assert np.all(pset.yi[:, fieldset.B.igrid] < ydim * 3) - #assert np.all(pset.yi[:, 0] < ydim) + # assert np.all(pset.xi[:, fieldset.B.igrid] < xdim * 4) + # assert np.all(pset.xi[:, 0] < xdim) + # assert pset.yi.shape[0] == len(pset.lon) + # assert pset.yi.shape[1] == fieldset.gridset.size + # assert np.all(pset.yi >= 0) + # assert np.all(pset.yi[:, fieldset.B.igrid] < ydim * 3) + # assert np.all(pset.yi[:, 0] < ydim) @pytest.mark.parametrize('mode', ['jit', 'scipy']) @@ -549,30 +549,6 @@ def test_multiple_grid_addlater_error(): assert fail -@pytest.mark.parametrize('mode', ['jit', 'scipy']) -def test_sampling_multiple_grid_sizes(mode): - """Sampling test that tests for FieldSet with different grid sizes - - While this currently works fine in Scipy mode, it fails in JIT mode with - an out_of_bounds_error because there is only one (xi, yi, zi) for each particle - A solution would be to define xi, yi, zi for each field separately - """ - xdim = 10 - ydim = 20 - gf = 10 # factor by which the resolution of U is higher than of V - U = Field('U', np.zeros((ydim*gf, xdim*gf), dtype=np.float32), - lon=np.linspace(0., 1., xdim*gf, dtype=np.float32), - lat=np.linspace(0., 1., ydim*gf, dtype=np.float32)) - V = Field('V', np.zeros((ydim, xdim), dtype=np.float32), - lon=np.linspace(0., 1., xdim, dtype=np.float32), - lat=np.linspace(0., 1., ydim, dtype=np.float32)) - fieldset = FieldSet(U, V) - pset = ParticleSet(fieldset, pclass=pclass(mode), lon=[0.8], lat=[0.9]) - - pset.execute(AdvectionRK4, runtime=10, dt=1) - assert np.isclose(pset[0].lon, 0.8) - - @pytest.mark.parametrize('mode', ['jit', 'scipy']) @pytest.mark.parametrize('with_W', [True, False]) @pytest.mark.parametrize('mesh', ['flat', 'spherical']) From 60c40f311db252b671cbd8560078f23deeb04773 Mon Sep 17 00:00:00 2001 From: Christian Kehl Date: Thu, 16 Jul 2020 17:07:21 +0200 Subject: [PATCH 12/56] list_of_pset_array_trials - added benchmark refactoring. tested and working (stommel). --- parcels/particleset_benchmark.py | 103 ++++++++++++++++++++++++------- performance/benchmark_stommel.py | 6 +- 2 files changed, 81 insertions(+), 28 deletions(-) diff --git a/parcels/particleset_benchmark.py b/parcels/particleset_benchmark.py index 69ed754d0c..32e6cc49db 100644 --- a/parcels/particleset_benchmark.py +++ b/parcels/particleset_benchmark.py @@ -55,10 +55,35 @@ class ParticleSet_TimingLog(): stime = 0 etime = 0 mtime = 0 - samples = [] - times_steps = [] + _samples = [] + _times_steps = [] _iter = 0 + def __init__(self): + self.stime = 0 + self.etime = 0 + self.mtime = 0 + self._samples = [] + self._times_steps = [] + self._iter = 0 + + @property + def timing(self): + return self._times_steps + + @property + def samples(self): + return self._samples + + def __len__(self): + return len(self._samples) + + def get_values(self): + return self._times_steps + + def get_value(self, index): + return self._times_steps[index] + def start_timing(self): if MPI: mpi_comm = MPI.COMM_WORLD @@ -97,33 +122,55 @@ def advance_iteration(self): mpi_comm = MPI.COMM_WORLD mpi_rank = mpi_comm.Get_rank() if mpi_rank == 0: - self.times_steps.append(self.mtime) - self.samples.append(self._iter) + self._times_steps.append(self.mtime) + self._samples.append(self._iter) self._iter += 1 self.mtime = 0 else: - self.times_steps.append(self.mtime) - self.samples.append(self._iter) + self._times_steps.append(self.mtime) + self._samples.append(self._iter) self._iter += 1 self.mtime = 0 class ParticleSet_ParamLogging(): - samples = [] - params = [] + _samples = [] + _params = [] _iter = 0 + def __init__(self): + self._samples = [] + self._params = [] + self._iter = 0 + + @property + def samples(self): + return self._samples + + @property + def params(self): + return self._params + + def get_params(self): + return self._params + + def get_param(self, index): + return self._params[index] + + def __len__(self): + return len(self._samples) + def advance_iteration(self, param): if MPI: mpi_comm = MPI.COMM_WORLD mpi_rank = mpi_comm.Get_rank() if mpi_rank == 0: - self.params.append(param) - self.samples.append(self._iter) + self._params.append(param) + self._samples.append(self._iter) self._iter += 1 else: - self.params.append(param) - self.samples.append(self._iter) + self._params.append(param) + self._samples.append(self._iter) self._iter += 1 @@ -174,8 +221,11 @@ def execute(self, pyfunc=AdvectionRK4, endtime=None, runtime=None, dt=1., :param movie_background_field: field plotted as background in the movie if moviedt is set. 'vector' shows the velocity as a vector field. :param verbose_progress: Boolean for providing a progress bar for the kernel execution loop. - :param postIterationCallbacks: (Optional) Array of functions that are to be called after each iteration (post-process, non-Kernel) - :param callbackdt: (Optional, in conjecture with 'postIterationCallbacks) timestep inverval to (latestly) interrupt the running kernel and invoke post-iteration callbacks from 'postIterationCallbacks' + :param postIterationCallbacks: (Optional) Array of functions that are to be called after each + iteration (post-process, non-Kernel) + :param callbackdt: (Optional, in conjecture with 'postIterationCallbacks) timestep inverval + to (latestly) interrupt the running kernel and invoke post-iteration + callbacks from 'postIterationCallbacks'. """ # check if pyfunc has changed since last compile. If so, recompile @@ -218,10 +268,11 @@ def execute(self, pyfunc=AdvectionRK4, endtime=None, runtime=None, dt=1., assert moviedt is None or moviedt >= 0, 'moviedt must be positive' # Set particle.time defaults based on sign of dt, if not set at ParticleSet construction - for p in self: - if np.isnan(p.time): - mintime, maxtime = self.fieldset.gridset.dimrange('time_full') - p.time = mintime if dt >= 0 else maxtime + for sublist in self._plist: + for p in sublist: + if np.isnan(p.time): + mintime, maxtime = self.fieldset.gridset.dimrange('time_full') + p.time = mintime if dt >= 0 else maxtime # Derive _starttime and endtime from arguments or fieldset defaults if runtime is not None and endtime is not None: @@ -229,7 +280,8 @@ def execute(self, pyfunc=AdvectionRK4, endtime=None, runtime=None, dt=1., # ====================================== # # ==== EXPENSIVE LIST COMPREHENSION ==== # # ====================================== # - _starttime = min([p.time for p in self]) if dt >= 0 else max([p.time for p in self]) + # _starttime = min([p.time for p in self]) if dt >= 0 else max([p.time for p in self]) + _starttime = min([p.time for sublist in self._plist for p in sublist]) if dt >= 0 else max([p.time for sublist in self._plist for p in sublist]) if self.repeatdt is not None and self.repeat_starttime is None: self.repeat_starttime = _starttime if runtime is not None: @@ -248,8 +300,10 @@ def execute(self, pyfunc=AdvectionRK4, endtime=None, runtime=None, dt=1., execute_once = True # Initialise particle timestepping - for p in self: - p.dt = dt + # for p in self: + for sublist in self._plist: + for p in sublist: + p.dt = dt # First write output_file, because particles could have been added if output_file: @@ -319,8 +373,11 @@ def execute(self, pyfunc=AdvectionRK4, endtime=None, runtime=None, dt=1., lat=self.repeatlat, depth=self.repeatdepth, pclass=self.repeatpclass, lonlatdepth_dtype=self.lonlatdepth_dtype, partitions=False, pid_orig=self.repeatpid, **self.repeatkwargs) - for p in pset_new: - p.dt = dt + # for p in pset_new: + # p.dt = dt + for sublist in pset_new.particles_a: + for p in sublist: + p.dt = dt self.add(pset_new) self.mem_io_log.stop_timing() self.mem_io_log.accumulate_timing() diff --git a/performance/benchmark_stommel.py b/performance/benchmark_stommel.py index a336cacf2d..eb0e125b64 100644 --- a/performance/benchmark_stommel.py +++ b/performance/benchmark_stommel.py @@ -6,8 +6,8 @@ from parcels import AdvectionEE, AdvectionRK45, AdvectionRK4 from parcels import FieldSet, ScipyParticle, JITParticle, Variable, AdvectionRK4, ErrorCode from parcels.particleset_benchmark import ParticleSet_Benchmark as ParticleSet -# from parcels.kernel_benchmark import Kernel_Benchmark as Kernel from parcels.field import VectorField, NestedField, SummedField +# from parcels.kernel_benchmark import Kernel_Benchmark as Kernel # from parcels import plotTrajectoriesFile_loadedField from datetime import timedelta as delta import math @@ -426,11 +426,8 @@ def Age(particle, fieldset, time): mpi_comm = MPI.COMM_WORLD mpi_rank = mpi_comm.Get_rank() if mpi_rank==0: - # endtime = ostime.time() - # endtime = MPI.Wtime() endtime = ostime.process_time() else: - #endtime = ostime.time() endtime = ostime.process_time() if args.write_out: @@ -494,7 +491,6 @@ def Age(particle, fieldset, time): if MPI: mpi_comm = MPI.COMM_WORLD - # mpi_comm.Barrier() Nparticles = mpi_comm.reduce(np.array(pset.nparticle_log.get_params()), op=MPI.SUM, root=0) Nmem = mpi_comm.reduce(np.array(pset.mem_log.get_params()), op=MPI.SUM, root=0) if mpi_comm.Get_rank() == 0: From 3dafd246e2dd5fca7c4eae567ffe6afa66deae7b Mon Sep 17 00:00:00 2001 From: Christian Kehl Date: Mon, 20 Jul 2020 23:06:58 +0200 Subject: [PATCH 13/56] list_of_pset_array_trials - rebase to benchmarking branch rev 54e4c18 --- parcels/kernel_benchmark.py | 159 +++++++++--------- parcels/particleset_benchmark.py | 123 -------------- performance/benchmark_CMEMS.py | 1 - .../benchmark_deep_migration_NPacific.py | 1 - performance/benchmark_galapagos_backwards.py | 1 - performance/benchmark_palaeo_Y2K.py | 6 +- performance/benchmark_perlin.py | 8 - performance/benchmark_stommel.py | 3 + tests/test_fieldset_sampling.py | 4 +- 9 files changed, 91 insertions(+), 215 deletions(-) diff --git a/parcels/kernel_benchmark.py b/parcels/kernel_benchmark.py index 07c430ab71..b64b129c14 100644 --- a/parcels/kernel_benchmark.py +++ b/parcels/kernel_benchmark.py @@ -84,8 +84,8 @@ def __del__(self): def execute_jit(self, pset, endtime, dt): """Invokes JIT engine to perform the core update loop""" self._io_timings.start_timing() - if len(pset.particles) > 0: - assert pset.fieldset.gridset.size == len(pset.particles[0].xi), \ + if len(pset.particles_a) > 0 and pset.particles_a[0].shape[0]: + assert pset.fieldset.gridset.size == len(pset.particles_a[0][0].xi), \ 'FieldSet has different amount of grids than Particle.xi. Have you added Fields after creating the ParticleSet?' for g in pset.fieldset.gridset.grids: g.cstruct = None # This force to point newly the grids from Python to C @@ -118,15 +118,16 @@ def execute_jit(self, pset, endtime, dt): fargs = [byref(f.ctypes_struct) for f in self.field_args.values()] fargs += [c_double(f) for f in self.const_args.values()] - particle_data = pset._particle_data.ctypes.data_as(c_void_p) self._mem_io_timings.stop_timing() self._mem_io_timings.accumulate_timing() self._compute_timings.start_timing() - self._function(c_int(len(pset)), particle_data, - c_double(endtime), - c_double(dt), - *fargs) + for subfield_c in pset.particles_c: + particle_data = subfield_c.ctypes.data_as(c_void_p) + self._function(c_int(subfield_c.shape[0]), particle_data, + c_double(endtime), c_double(dt), + *fargs) + self._compute_timings.stop_timing() self._compute_timings.accumulate_timing() @@ -155,79 +156,81 @@ def execute_python(self, pset, endtime, dt): self._mem_io_timings.accumulate_timing() self._compute_timings.start_timing() - for p in pset.particles: - ptype = p.getPType() - # Don't execute particles that aren't started yet - sign_end_part = np.sign(endtime - p.time) - - dt_pos = min(abs(p.dt), abs(endtime - p.time)) - - # ==== numerically stable; also making sure that continuously-recovered particles do end successfully, - # as they fulfil the condition here on entering at the final calculation here. ==== # - if ((sign_end_part != sign_dt) or np.isclose(dt_pos, 0)) and not np.isclose(dt, 0): - if abs(p.time) >= abs(endtime): - p.state = ErrorCode.Success - continue + # for p in pset.particles: + for subfield in pset.particles_a: + for p in subfield: + ptype = p.getPType() + # Don't execute particles that aren't started yet + sign_end_part = np.sign(endtime - p.time) + + dt_pos = min(abs(p.dt), abs(endtime - p.time)) + + # ==== numerically stable; also making sure that continuously-recovered particles do end successfully, + # as they fulfil the condition here on entering at the final calculation here. ==== # + if ((sign_end_part != sign_dt) or np.isclose(dt_pos, 0)) and not np.isclose(dt, 0): + if abs(p.time) >= abs(endtime): + p.state = ErrorCode.Success + continue + + while p.state in [ErrorCode.Evaluate, ErrorCode.Repeat] or np.isclose(dt, 0): - while p.state in [ErrorCode.Evaluate, ErrorCode.Repeat] or np.isclose(dt, 0): - - for var in ptype.variables: - p_var_back[var.name] = getattr(p, var.name) - try: - pdt_prekernels = sign_dt * dt_pos - p.dt = pdt_prekernels - state_prev = p.state - res = self.pyfunc(p, pset.fieldset, p.time) - if res is None: - res = ErrorCode.Success - - if res is ErrorCode.Success and p.state != state_prev: - res = p.state - - if res == ErrorCode.Success and not np.isclose(p.dt, pdt_prekernels): - res = ErrorCode.Repeat - - except FieldOutOfBoundError as fse_xy: - res = ErrorCode.ErrorOutOfBounds - p.exception = fse_xy - except FieldOutOfBoundSurfaceError as fse_z: - res = ErrorCode.ErrorThroughSurface - p.exception = fse_z - except TimeExtrapolationError as fse_t: - res = ErrorCode.ErrorTimeExtrapolation - p.exception = fse_t - except Exception as e: - res = ErrorCode.Error - p.exception = e - - # Handle particle time and time loop - if res in [ErrorCode.Success, ErrorCode.Delete]: - # Update time and repeat - p.time += p.dt - p.update_next_dt() - dt_pos = min(abs(p.dt), abs(endtime - p.time)) - - sign_end_part = np.sign(endtime - p.time) - if res != ErrorCode.Delete and not np.isclose(dt_pos, 0) and (sign_end_part == sign_dt): - res = ErrorCode.Evaluate - if sign_end_part != sign_dt: - dt_pos = 0 - - p.state = res - if np.isclose(dt, 0): - break - else: - p.state = res - # Try again without time update for var in ptype.variables: - if var.name not in ['dt', 'state']: - setattr(p, var.name, p_var_back[var.name]) - dt_pos = min(abs(p.dt), abs(endtime - p.time)) - - sign_end_part = np.sign(endtime - p.time) - if sign_end_part != sign_dt: - dt_pos = 0 - break + p_var_back[var.name] = getattr(p, var.name) + try: + pdt_prekernels = sign_dt * dt_pos + p.dt = pdt_prekernels + state_prev = p.state + res = self.pyfunc(p, pset.fieldset, p.time) + if res is None: + res = ErrorCode.Success + + if res is ErrorCode.Success and p.state != state_prev: + res = p.state + + if res == ErrorCode.Success and not np.isclose(p.dt, pdt_prekernels): + res = ErrorCode.Repeat + + except FieldOutOfBoundError as fse_xy: + res = ErrorCode.ErrorOutOfBounds + p.exception = fse_xy + except FieldOutOfBoundSurfaceError as fse_z: + res = ErrorCode.ErrorThroughSurface + p.exception = fse_z + except TimeExtrapolationError as fse_t: + res = ErrorCode.ErrorTimeExtrapolation + p.exception = fse_t + except Exception as e: + res = ErrorCode.Error + p.exception = e + + # Handle particle time and time loop + if res in [ErrorCode.Success, ErrorCode.Delete]: + # Update time and repeat + p.time += p.dt + p.update_next_dt() + dt_pos = min(abs(p.dt), abs(endtime - p.time)) + + sign_end_part = np.sign(endtime - p.time) + if res != ErrorCode.Delete and not np.isclose(dt_pos, 0) and (sign_end_part == sign_dt): + res = ErrorCode.Evaluate + if sign_end_part != sign_dt: + dt_pos = 0 + + p.state = res + if np.isclose(dt, 0): + break + else: + p.state = res + # Try again without time update + for var in ptype.variables: + if var.name not in ['dt', 'state']: + setattr(p, var.name, p_var_back[var.name]) + dt_pos = min(abs(p.dt), abs(endtime - p.time)) + + sign_end_part = np.sign(endtime - p.time) + if sign_end_part != sign_dt: + dt_pos = 0 + break self._compute_timings.stop_timing() self._compute_timings.accumulate_timing() diff --git a/parcels/particleset_benchmark.py b/parcels/particleset_benchmark.py index 32e6cc49db..6e5c53c504 100644 --- a/parcels/particleset_benchmark.py +++ b/parcels/particleset_benchmark.py @@ -51,129 +51,6 @@ def measure_mem_usage(): USE_ASYNC_MEMLOG = False USE_RUSE_SYNC_MEMLOG = False # can be faulty -class ParticleSet_TimingLog(): - stime = 0 - etime = 0 - mtime = 0 - _samples = [] - _times_steps = [] - _iter = 0 - - def __init__(self): - self.stime = 0 - self.etime = 0 - self.mtime = 0 - self._samples = [] - self._times_steps = [] - self._iter = 0 - - @property - def timing(self): - return self._times_steps - - @property - def samples(self): - return self._samples - - def __len__(self): - return len(self._samples) - - def get_values(self): - return self._times_steps - - def get_value(self, index): - return self._times_steps[index] - - def start_timing(self): - if MPI: - mpi_comm = MPI.COMM_WORLD - mpi_rank = mpi_comm.Get_rank() - if mpi_rank == 0: - # self.stime = MPI.Wtime() - # self.stime = time_module.perf_counter() - self.stime = time_module.process_time() - else: - self.stime = time_module.perf_counter() - - def stop_timing(self): - if MPI: - mpi_comm = MPI.COMM_WORLD - mpi_rank = mpi_comm.Get_rank() - if mpi_rank == 0: - # self.etime = MPI.Wtime() - # self.etime = time_module.perf_counter() - self.etime = time_module.process_time() - else: - self.etime = time_module.perf_counter() - - def accumulate_timing(self): - if MPI: - mpi_comm = MPI.COMM_WORLD - mpi_rank = mpi_comm.Get_rank() - if mpi_rank == 0: - self.mtime += (self.etime-self.stime) - else: - self.mtime = 0 - else: - self.mtime += (self.etime-self.stime) - - def advance_iteration(self): - if MPI: - mpi_comm = MPI.COMM_WORLD - mpi_rank = mpi_comm.Get_rank() - if mpi_rank == 0: - self._times_steps.append(self.mtime) - self._samples.append(self._iter) - self._iter += 1 - self.mtime = 0 - else: - self._times_steps.append(self.mtime) - self._samples.append(self._iter) - self._iter += 1 - self.mtime = 0 - - -class ParticleSet_ParamLogging(): - _samples = [] - _params = [] - _iter = 0 - - def __init__(self): - self._samples = [] - self._params = [] - self._iter = 0 - - @property - def samples(self): - return self._samples - - @property - def params(self): - return self._params - - def get_params(self): - return self._params - - def get_param(self, index): - return self._params[index] - - def __len__(self): - return len(self._samples) - - def advance_iteration(self, param): - if MPI: - mpi_comm = MPI.COMM_WORLD - mpi_rank = mpi_comm.Get_rank() - if mpi_rank == 0: - self._params.append(param) - self._samples.append(self._iter) - self._iter += 1 - else: - self._params.append(param) - self._samples.append(self._iter) - self._iter += 1 - - class ParticleSet_Benchmark(ParticleSet): def __init__(self, fieldset, pclass=JITParticle, lon=None, lat=None, depth=None, time=None, repeatdt=None, diff --git a/performance/benchmark_CMEMS.py b/performance/benchmark_CMEMS.py index 6959eed652..af6af774b0 100644 --- a/performance/benchmark_CMEMS.py +++ b/performance/benchmark_CMEMS.py @@ -366,7 +366,6 @@ def perIterGC(): if MPI: mpi_comm = MPI.COMM_WORLD - # mpi_comm.Barrier() Nparticles = mpi_comm.reduce(np.array(pset.nparticle_log.get_params()), op=MPI.SUM, root=0) Nmem = mpi_comm.reduce(np.array(pset.mem_log.get_params()), op=MPI.SUM, root=0) if mpi_comm.Get_rank() == 0: diff --git a/performance/benchmark_deep_migration_NPacific.py b/performance/benchmark_deep_migration_NPacific.py index 17b5fa2973..2edff5bb07 100644 --- a/performance/benchmark_deep_migration_NPacific.py +++ b/performance/benchmark_deep_migration_NPacific.py @@ -442,7 +442,6 @@ def Profiles(particle, fieldset, time): if MPI: mpi_comm = MPI.COMM_WORLD - # mpi_comm.Barrier() Nparticles = mpi_comm.reduce(np.array(pset.nparticle_log.get_params()), op=MPI.SUM, root=0) Nmem = mpi_comm.reduce(np.array(pset.mem_log.get_params()), op=MPI.SUM, root=0) if mpi_comm.Get_rank() == 0: diff --git a/performance/benchmark_galapagos_backwards.py b/performance/benchmark_galapagos_backwards.py index 317292cbaf..be53489051 100644 --- a/performance/benchmark_galapagos_backwards.py +++ b/performance/benchmark_galapagos_backwards.py @@ -208,7 +208,6 @@ def periodicBC(particle, fieldSet, time): if MPI: mpi_comm = MPI.COMM_WORLD - # mpi_comm.Barrier() Nparticles = mpi_comm.reduce(np.array(pset.nparticle_log.get_params()), op=MPI.SUM, root=0) Nmem = mpi_comm.reduce(np.array(pset.mem_log.get_params()), op=MPI.SUM, root=0) if mpi_comm.Get_rank() == 0: diff --git a/performance/benchmark_palaeo_Y2K.py b/performance/benchmark_palaeo_Y2K.py index 4150891ee2..981b4b9bd8 100644 --- a/performance/benchmark_palaeo_Y2K.py +++ b/performance/benchmark_palaeo_Y2K.py @@ -399,8 +399,10 @@ class DinoParticle(JITParticle): # pset.execute(kernels, runtime=delta(days=365*9), dt=delta(minutes=-20), output_file=pfile, verbose_progress=False, # recovery={ErrorCode.ErrorOutOfBounds: DeleteParticle}, postIterationCallbacks=postProcessFuncs) # postIterationCallbacks=postProcessFuncs, callbackdt=delta(hours=12) - pset.execute(kernels, runtime=delta(days=time_in_days), dt=delta(hours=-12), output_file=pfile, verbose_progress=False, recovery={ErrorCode.ErrorOutOfBounds: DeleteParticle}, postIterationCallbacks=postProcessFuncs, callbackdt=np.infty) - + pset.execute(kernels, runtime=delta(days=time_in_days), dt=delta(hours=-12), output_file=pfile, + verbose_progress=False, recovery={ErrorCode.ErrorOutOfBounds: DeleteParticle}, + postIterationCallbacks=postProcessFuncs, callbackdt=np.infty) + if MPI: mpi_comm = MPI.COMM_WORLD mpi_rank = mpi_comm.Get_rank() diff --git a/performance/benchmark_perlin.py b/performance/benchmark_perlin.py index 67622b4909..e0f62c2062 100644 --- a/performance/benchmark_perlin.py +++ b/performance/benchmark_perlin.py @@ -91,18 +91,13 @@ def perlin_fieldset_from_numpy(periodic_wrap=False, write_out=False): # Coordinates of the test fieldset (on A-grid in deg) lon = np.linspace(-a*0.5, a*0.5, img_shape[0], dtype=np.float32) - # sys.stdout.write("lon field: {}\n".format(lon.size)) lat = np.linspace(-b*0.5, b*0.5, img_shape[1], dtype=np.float32) - # sys.stdout.write("lat field: {}\n".format(lat.size)) totime = tsteps*tscale*24.0*60.0*60.0 time = np.linspace(0., totime, tsteps, dtype=np.float64) - # sys.stdout.write("time field: {}\n".format(time.size)) # Define arrays U (zonal), V (meridional) U = perlin2d.generate_fractal_noise_temporal2d(img_shape, tsteps, (perlinres[1], perlinres[2]), noctaves, perlin_persistence, max_shift=((-1, 2), (-1, 2))) U = np.transpose(U, (0,2,1)) - # U = np.swapaxes(U, 1, 2) - # print("U-statistics - min: {:10.7f}; max: {:10.7f}; avg. {:10.7f}; std_dev: {:10.7f}".format(U.min(initial=0), U.max(initial=0), U.mean(), U.std())) V = perlin2d.generate_fractal_noise_temporal2d(img_shape, tsteps, (perlinres[1], perlinres[2]), noctaves, perlin_persistence, max_shift=((-1, 2), (-1, 2))) V = np.transpose(V, (0,2,1)) # V = np.swapaxes(V, 1, 2) @@ -450,12 +445,9 @@ def Age(particle, fieldset, time): if MPI: mpi_comm = MPI.COMM_WORLD - # mpi_comm.Barrier() Nparticles = mpi_comm.reduce(np.array(pset.nparticle_log.get_params()), op=MPI.SUM, root=0) Nmem = mpi_comm.reduce(np.array(pset.mem_log.get_params()), op=MPI.SUM, root=0) if mpi_comm.Get_rank() == 0: pset.plot_and_log(memory_used=Nmem, nparticles=Nparticles, target_N=target_N, imageFilePath=imageFileName, odir=odir, xlim_range=[0, 730], ylim_range=[0, 120]) else: pset.plot_and_log(target_N=target_N, imageFilePath=imageFileName, odir=odir, xlim_range=[0, 730], ylim_range=[0, 150]) - - diff --git a/performance/benchmark_stommel.py b/performance/benchmark_stommel.py index eb0e125b64..6bf63192c4 100644 --- a/performance/benchmark_stommel.py +++ b/performance/benchmark_stommel.py @@ -426,8 +426,11 @@ def Age(particle, fieldset, time): mpi_comm = MPI.COMM_WORLD mpi_rank = mpi_comm.Get_rank() if mpi_rank==0: + # endtime = ostime.time() + # endtime = MPI.Wtime() endtime = ostime.process_time() else: + #endtime = ostime.time() endtime = ostime.process_time() if args.write_out: diff --git a/tests/test_fieldset_sampling.py b/tests/test_fieldset_sampling.py index 8fdb9ee9aa..54fb346fa5 100644 --- a/tests/test_fieldset_sampling.py +++ b/tests/test_fieldset_sampling.py @@ -521,7 +521,9 @@ def test_sampling_multiple_grid_sizes(mode, ugridfactor): assert fieldset.U.grid is fieldset.V.grid pset.execute(AdvectionRK4, runtime=10, dt=1) assert np.isclose(pset[0].lon, 0.8) - assert np.all((pset[0].xi >= 0) & (pset[0].xi < xdim*ugridfactor)) + if mode == 'jit': + # assert np.all((pset[0].xi >= 0) & (pset[0].xi < xdim*ugridfactor)) + assert np.all(np.array([[p.xi[i] >= 0 & p.xi[i] < xdim*ugridfactor for i in range(0, len(p.xi))] for p in pset])) def test_multiple_grid_addlater_error(): From 29574a4f39c724bda4ad6363d0ef497e470be198 Mon Sep 17 00:00:00 2001 From: Christian Kehl Date: Tue, 1 Dec 2020 20:28:04 +0100 Subject: [PATCH 14/56] list_of_pset_array_trials - fixed issues in file writing --- parcels/particlefile.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/parcels/particlefile.py b/parcels/particlefile.py index ea4895b8a1..168285406c 100644 --- a/parcels/particlefile.py +++ b/parcels/particlefile.py @@ -240,7 +240,8 @@ def convert_pset_to_dict(self, pset, time, deleted_only=False): pset_towrite = [p for sublist in pset if sublist is not None - for p in sublist] + for p in sublist + if p.state == ErrorCode.Delete] # pset_towrite = pset # == commented due to git rebase to master 27 02 2020 == # # elif pset[0].dt > 0: From f418419e117fce5cf2913e7653234d569f8a8eb3 Mon Sep 17 00:00:00 2001 From: Christian Kehl Date: Wed, 9 Dec 2020 21:48:22 +0100 Subject: [PATCH 15/56] origin/list_of_pset_array_trials - added a version where the 1st-level structure is actually a deque. --- parcels/particleset_nodeloa.py | 1093 ++++++++++++++++++++++++++++++++ 1 file changed, 1093 insertions(+) create mode 100644 parcels/particleset_nodeloa.py diff --git a/parcels/particleset_nodeloa.py b/parcels/particleset_nodeloa.py new file mode 100644 index 0000000000..f209cb1c71 --- /dev/null +++ b/parcels/particleset_nodeloa.py @@ -0,0 +1,1093 @@ +import collections +import time as time_module +from datetime import date +from datetime import datetime +from datetime import timedelta as delta + +import numpy as np +import xarray as xr +import progressbar + +from parcels.compiler import GNUCompiler +from parcels.field import NestedField +from parcels.field import SummedField +from parcels.grid import GridCode +from parcels.kernel import Kernel +from parcels.kernels.advection import AdvectionRK4 +from parcels.particle import JITParticle +from parcels.particlefile import ParticleFile +from parcels.tools.loggers import logger +try: + from mpi4py import MPI +except: + MPI = None +if MPI: + try: + from sklearn.cluster import KMeans + except: + raise EnvironmentError('sklearn needs to be available if MPI is installed. ' + 'See http://oceanparcels.org/#parallel_install for more information') + +__all__ = ['ParticleSet'] + + +class ParticleSet(object): + """Container class for storing particle and executing kernel over them. + + Please note that this currently only supports fixed size particle sets. + + :param fieldset: :mod:`parcels.fieldset.FieldSet` object from which to sample velocity + :param pclass: Optional :mod:`parcels.particle.JITParticle` or + :mod:`parcels.particle.ScipyParticle` object that defines custom particle + :param lon: List of initial longitude values for particles + :param lat: List of initial latitude values for particles + :param depth: Optional list of initial depth values for particles. Default is 0m + :param time: Optional list of initial time values for particles. Default is fieldset.U.grid.time[0] + :param repeatdt: Optional interval (in seconds) on which to repeat the release of the ParticleSet + :param lonlatdepth_dtype: Floating precision for lon, lat, depth particle coordinates. + It is either np.float32 or np.float64. Default is np.float32 if fieldset.U.interp_method is 'linear' + and np.float64 if the interpolation method is 'cgrid_velocity' + :param partitions: List of cores on which to distribute the particles for MPI runs. Default: None, in which case particles + are distributed automatically on the processors + Other Variables can be initialised using further arguments (e.g. v=... for a Variable named 'v') + """ + + def __init__(self, fieldset, pclass=JITParticle, lon=None, lat=None, depth=None, time=None, repeatdt=None, lonlatdepth_dtype=None, pid_orig=None, **kwargs): + self.fieldset = fieldset + self.fieldset.check_complete() + partitions = kwargs.pop('partitions', None) + + def convert_to_array(var): + # Convert lists and single integers/floats to one-dimensional numpy arrays + if isinstance(var, np.ndarray): + return var.flatten() + elif isinstance(var, (int, float, np.float32, np.int32)): + return np.array([var]) + else: + return np.array(var) + + lon = np.empty(shape=0) if lon is None else convert_to_array(lon) + lat = np.empty(shape=0) if lat is None else convert_to_array(lat) + if pid_orig is None: + pid_orig = np.arange(lon.size) + pid = pid_orig + pclass.lastID + + if depth is None: + mindepth, _ = self.fieldset.gridset.dimrange('depth') + depth = np.ones(lon.size) * mindepth + else: + depth = convert_to_array(depth) + assert lon.size == lat.size and lon.size == depth.size, ( + 'lon, lat, depth don''t all have the same lenghts') + + time = convert_to_array(time) + time = np.repeat(time, lon.size) if time.size == 1 else time + if time.size > 0 and type(time[0]) in [datetime, date]: + time = np.array([np.datetime64(t) for t in time]) + self.time_origin = fieldset.time_origin + if time.size > 0 and isinstance(time[0], np.timedelta64) and not self.time_origin: + raise NotImplementedError('If fieldset.time_origin is not a date, time of a particle must be a double') + time = np.array([self.time_origin.reltime(t) if isinstance(t, np.datetime64) else t for t in time]) + assert lon.size == time.size, ( + 'time and positions (lon, lat, depth) don''t have the same lengths.') + + if partitions is not None and partitions is not False: + partitions = convert_to_array(partitions) + + for kwvar in kwargs: + kwargs[kwvar] = convert_to_array(kwargs[kwvar]) + assert lon.size == kwargs[kwvar].size, ( + '%s and positions (lon, lat, depth) don''t have the same lengths.' % kwargs[kwvar]) + + offset = np.max(pid) if len(pid) > 0 else -1 + if MPI: + mpi_comm = MPI.COMM_WORLD + mpi_rank = mpi_comm.Get_rank() + mpi_size = mpi_comm.Get_size() + + if lon.size < mpi_size and mpi_size > 1: + raise RuntimeError('Cannot initialise with fewer particles than MPI processors') + + if mpi_size > 1: + if partitions is not False: + if partitions is None: + if mpi_rank == 0: + coords = np.vstack((lon, lat)).transpose() + kmeans = KMeans(n_clusters=mpi_size, random_state=0).fit(coords) + partitions = kmeans.labels_ + else: + partitions = None + partitions = mpi_comm.bcast(partitions, root=0) + elif np.max(partitions >= mpi_rank): + raise RuntimeError('Particle partitions must vary between 0 and the number of mpi procs') + lon = lon[partitions == mpi_rank] + lat = lat[partitions == mpi_rank] + time = time[partitions == mpi_rank] + depth = depth[partitions == mpi_rank] + pid = pid[partitions == mpi_rank] + for kwvar in kwargs: + kwargs[kwvar] = kwargs[kwvar][partitions == mpi_rank] + offset = MPI.COMM_WORLD.allreduce(offset, op=MPI.MAX) + + self.repeatdt = repeatdt.total_seconds() if isinstance(repeatdt, delta) else repeatdt + if self.repeatdt: + if self.repeatdt <= 0: + raise('Repeatdt should be > 0') + if time[0] and not np.allclose(time, time[0]): + raise ('All Particle.time should be the same when repeatdt is not None') + self.repeat_starttime = time[0] + self.repeatlon = lon + self.repeatlat = lat + if not hasattr(self, 'repeatpid'): + self.repeatpid = pid - pclass.lastID + self.repeatdepth = depth + self.repeatpclass = pclass + self.partitions = partitions + self.repeatkwargs = kwargs + pclass.setLastID(offset+1) + + if lonlatdepth_dtype is None: + self.lonlatdepth_dtype = self.lonlatdepth_dtype_from_field_interp_method(fieldset.U) + else: + self.lonlatdepth_dtype = lonlatdepth_dtype + assert self.lonlatdepth_dtype in [np.float32, np.float64], \ + 'lon lat depth precision should be set to either np.float32 or np.float64' + JITParticle.set_lonlatdepth_dtype(self.lonlatdepth_dtype) + + # == this should be a list not a dict! == # + self._pid_mapping_bounds = [] # plist bracket_index -> (min_id, max_id, size_bracket) + # ======================================= # + self.nlist_limit = 4096 + # self.particles = np.empty(lon.size, dtype=pclass) + assert lon.shape[0] == lat.shape[0], ('Length of lon and lat do not match.') + # self._plist = [] + self._plist = collections.deque([]) + self._bracket_no = 0 + start_index = 0 + end_index = 0 + while end_index < lon.shape[0]: + end_index = min(lon.shape[0], start_index + self.nlist_limit) + self._plist.append(np.empty(end_index-start_index, dtype=pclass)) + self._pid_mapping_bounds.append((np.iinfo(np.int32).max, np.iinfo(np.int32).min, end_index-start_index)) + # bracket_index = len(self._plist)-1 + start_index = end_index + # end_index = min(lon.shape[0], end_index+self.nlist_limit) + # self._particle_data = None + self._plist_c = None + self._pclass = pclass + self._ptype = self.pclass.getPType() + self.kernel = None + + # ==== ANNOTATION ==== # + # the particles themselves (self.particles, dtype: pclass) are already separate memory fields from + # their c-pointers (self._particle_data, dtype: PType.dtype [ie. padded memory block]). Thus, it makes little + # sense to require them being static over the whole ParticleSet lifetime. In other words: there should be a mechanism + # to update the particles, reallocate the c-pointer field, and then just value-copy the particle data into the c-memory. + # nparticles = self.allocate_cptrs() + if self._ptype.uses_jit: + # self._plist_c = [] + self._plist_c = collections.deque([]) + for bracket_index in range(len(self._plist)): + self._plist_c.append(np.empty(self._pid_mapping_bounds[bracket_index][2], dtype=self.ptype.dtype)) + + if lon is not None and lat is not None: + # == Initialise from arrays of lon/lat coordinates == # + # assert self.particles.size == lon.size and self.particles.size == lat.size, ('Size of ParticleSet does not match length of lon and lat.') + + nbrackets = 0 + for i in range(lon.size): + bracket_index = int(float(i)/self.nlist_limit) + slot_index = int(i % self.nlist_limit) + # self._plist[bracket_index][slot_index] = pclass(lon[i], lat[i], pid[i], fieldset=fieldset, depth=depth[i], cptr=self.cptr(bracket_index, slot_index), time=time[i]) + self.current_p_node()[slot_index] = pclass(lon[i], lat[i], pid[i], fieldset=fieldset, depth=depth[i], cptr=self.current_c_node()(slot_index), time=time[i]) + bracket_info = self._pid_mapping_bounds[bracket_index] + self._pid_mapping_bounds[bracket_index] = (min(pid[i], bracket_info[0]), max(pid[i], bracket_info[1]), bracket_info[2]) + # self.particles[i] = pclass(lon[i], lat[i], pid[i], fieldset=fieldset, depth=depth[i], cptr=cptr(i), time=time[i]) + # == Set other Variables if provided == # + for kwvar in kwargs: + # if not hasattr(self._plist[bracket_index][slot_index], kwvar): + if not hasattr(self.current_p_node()[slot_index], kwvar): + raise RuntimeError('Particle class does not have Variable %s' % kwvar) + # setattr(self._plist[bracket_index][slot_index], kwvar, kwargs[kwvar][i]) + setattr(self.current_p_node()[slot_index], kwvar, kwargs[kwvar][i]) + # if not hasattr(self.particles[i], kwvar): + # raise RuntimeError('Particle class does not have Variable %s' % kwvar) + # setattr(self.particles[i], kwvar, kwargs[kwvar][i]) + if slot_index >= (self.nlist_limit-1): + self.next() + self.reset_current_node() + else: + raise ValueError("Latitude and longitude required for generating ParticleSet") + + def next(self): + self._bracket_no += 1 + self._plist.rotate(-1) + if self._plist_c is not None and self._ptype.uses_jit: + self._plist_c.rotate(-1) + + def prev(self): + self._bracket_no -= 1 + self._plist.rotate(1) + if self._plist_c is not None and self._ptype.uses_jit: + self._plist_c.rotate(1) + + def reset_current_node(self): + self._plist.rotate(self._bracket_no) + if self._plist_c is not None and self._ptype.uses_jit: + self._plist_c.rotate(self._bracket_no) + self._bracket_no = 0 + + @property + def current_p_node(self): + return self._plist[0] + + @current_p_node.setter + def current_p_node(self, value): + self._plist[0] = value + + @property + def current_c_node(self): + if self._plist_c is not None and self._ptype.uses_jit: + return self._plist_c[0] + return None + + @current_c_node.setter + def current_c_node(self, value): + if self._plist_c is not None and self._ptype.uses_jit: + self._plist_c[0] = value + + @property + def particles(self): + if len(self._plist) == 0: + return None + result = self._plist[0] + if len(self._plist) > 1: + # for bracket_index in range(0, len(self._plist)): + # result = np.concatenate((result, self._plist[bracket_index]), axis=0) + for bracket in self._plist: + result = np.concatenate((result, bracket), axis=0) + return result + # return self._plist + + @property + def particles_a(self): + return self._plist + + @property + def pid_mapping_bounds(self): + return self._pid_mapping_bounds + + @property + def particles_c(self): + return self._plist_c + + @property + def ptype(self): + return self._ptype + + @property + def pclass(self): + return self._pclass + + def cptr(self, bracket_index, slot_index): + """ + THIS IS EXTREMELY SLOW + """ + # if self._particle_data is None: + if self._plist_c is None: + return None + if self.ptype.uses_jit: + return self._plist_c[bracket_index][slot_index] + else: + return None + + # def allocate_cptrs(self): + # nparticles = 0 + # for sublist in self.plist: + # nparticles += sublist.shape[0] + # # == Allocate underlying data for C-allocated particles == # + # if self.ptype.uses_jit: + # self._particle_data = np.empty(nparticles, dtype=self.ptype.dtype) + # return nparticles + + def get_list_array_index(self, pdata): + """Searches for the list-array indices for a given particle. + + :param pdata: :mod:`int` global index or :mod:`parcels.particle._Particle` (or its subclass) object to search for + :return :mod:`tuple` in format (, ), to index + """ + if isinstance(pdata, int): + # == parameter is the global particle index == # + bracket_index = 0 + start_index = 0 + end_index = self._pid_mapping_bounds[bracket_index][2] + while pdata >= end_index: + start_index = end_index + end_index = start_index+self._pid_mapping_bounds[bracket_index][2] + bracket_index += 1 + slot_index = pdata-start_index + return (bracket_index, slot_index) + elif isinstance(pdata, self.pclass): + # == parameter is a particle itself == # + bracket_index = 0 + bracket_info = self.pid_mapping_bounds[bracket_index] + while bracket_index < len(self._plist) and ((pdata.id < bracket_info[0]) or (pdata.id > bracket_info[1])): + bracket_index += 1 + bracket_info = self.pid_mapping_bounds[bracket_index] + self.next() + if bracket_index >= len(self._plist): + logger.warn_once("ParticleSet.get_cptr_index() - requested particle ({}) not found.".format(pdata)) + self.reset_current_node() + return -1 + slot_index = 0 + while (slot_index < self.current_p_node().shape[0]) and (self.current_p_node()[slot_index].id != pdata.id): + slot_index += 1 + if slot_index >= self.current_p_node().shape[0]: + logger.warn_once("ParticleSet.get_cptr_index() - requested particle ({}) not found.".format(pdata)) + self.reset_current_node() + return -1 + self.reset_current_node() + return (bracket_index, slot_index) + + def get_list_array_index_by_PID(self, pdata): + # == parameter is the particle ID == # + located_particle = [(p, bracket_index, slot_index) for bracket_index, sublist in enumerate(self._plist) for slot_index, p in enumerate(sublist) if p.id == pdata] + if len(located_particle) < 1: + logger.warn_once("ParticleSet.get_cptr_index() - requested particle ID ({}) not existent.".format(pdata)) + return -1 + if len(located_particle) > 1: + logger.warn_once("ParticleSet.get_cptr_index() - requested particle ID ({}) ambiguous.".format(pdata)) + return -1 + located_particle = located_particle[0] + # nparticle = 0 + # for i in range(located_particle[1]): + # nparticle += self.pid_mapping_bounds[i][2] + # return nparticle+located_particle[2] + return (located_particle[1], located_particle[2]) + + def covert_static_array_to_list_array(self, property_array): + if not isinstance(property_array, np.ndarray): + property_array = np.array(property_array) + nparticles = 0 + for bracket_index in range(len(self._pid_mapping_bounds)): + nparticles += self._pid_mapping_bounds[bracket_index][2] + assert property_array.shape[0] == nparticles, ("Attaching the property array to the list of particle prohibited because num. of entries do not match.") + results = [] + bracket_index = 0 + start_index = 0 + end_index = self._pid_mapping_bounds[bracket_index][2] + while start_index < property_array.shape[0]: + results.append(np.array(property_array[start_index:min(end_index, property_array.shape[0])])) + start_index = end_index + end_index = start_index+self._pid_mapping_bounds[bracket_index][2] + bracket_index += 1 + return results + + @staticmethod + def convert_list_of_arrays_to_static_list(pset): + return [p + for sublist in pset + if sublist is not None + for p in sublist] + + def convert_pset_to_static_list(self): + return [p + for sublist in self.particles_a + if sublist is not None + for p in sublist] + + def convert_pset_to_static_array(self): + return np.array(self.convert_pset_to_static_list(), dtype=self._pclass) + + @classmethod + def from_list(cls, fieldset, pclass, lon, lat, depth=None, time=None, repeatdt=None, lonlatdepth_dtype=None, **kwargs): + """Initialise the ParticleSet from lists of lon and lat + + :param fieldset: :mod:`parcels.fieldset.FieldSet` object from which to sample velocity + :param pclass: mod:`parcels.particle.JITParticle` or :mod:`parcels.particle.ScipyParticle` + object that defines custom particle + :param lon: List of initial longitude values for particles + :param lat: List of initial latitude values for particles + :param depth: Optional list of initial depth values for particles. Default is 0m + :param time: Optional list of start time values for particles. Default is fieldset.U.time[0] + :param repeatdt: Optional interval (in seconds) on which to repeat the release of the ParticleSet + :param lonlatdepth_dtype: Floating precision for lon, lat, depth particle coordinates. + It is either np.float32 or np.float64. Default is np.float32 if fieldset.U.interp_method is 'linear' + and np.float64 if the interpolation method is 'cgrid_velocity' + Other Variables can be initialised using further arguments (e.g. v=... for a Variable named 'v') + """ + return cls(fieldset=fieldset, pclass=pclass, lon=lon, lat=lat, depth=depth, time=time, repeatdt=repeatdt, lonlatdepth_dtype=lonlatdepth_dtype, **kwargs) + + @classmethod + def from_line(cls, fieldset, pclass, start, finish, size, depth=None, time=None, repeatdt=None, lonlatdepth_dtype=None): + """Initialise the ParticleSet from start/finish coordinates with equidistant spacing + Note that this method uses simple numpy.linspace calls and does not take into account + great circles, so may not be a exact on a globe + + :param fieldset: :mod:`parcels.fieldset.FieldSet` object from which to sample velocity + :param pclass: mod:`parcels.particle.JITParticle` or :mod:`parcels.particle.ScipyParticle` + object that defines custom particle + :param start: Starting point for initialisation of particles on a straight line. + :param finish: End point for initialisation of particles on a straight line. + :param size: Initial size of particle set + :param depth: Optional list of initial depth values for particles. Default is 0m + :param time: Optional start time value for particles. Default is fieldset.U.time[0] + :param repeatdt: Optional interval (in seconds) on which to repeat the release of the ParticleSet + :param lonlatdepth_dtype: Floating precision for lon, lat, depth particle coordinates. + It is either np.float32 or np.float64. Default is np.float32 if fieldset.U.interp_method is 'linear' + and np.float64 if the interpolation method is 'cgrid_velocity' + """ + + lon = np.linspace(start[0], finish[0], size) + lat = np.linspace(start[1], finish[1], size) + if type(depth) in [int, float]: + depth = [depth] * size + return cls(fieldset=fieldset, pclass=pclass, lon=lon, lat=lat, depth=depth, time=time, repeatdt=repeatdt, lonlatdepth_dtype=lonlatdepth_dtype) + + @classmethod + def from_field(cls, fieldset, pclass, start_field, size, mode='monte_carlo', depth=None, time=None, repeatdt=None, lonlatdepth_dtype=None): + """Initialise the ParticleSet randomly drawn according to distribution from a field + + :param fieldset: :mod:`parcels.fieldset.FieldSet` object from which to sample velocity + :param pclass: mod:`parcels.particle.JITParticle` or :mod:`parcels.particle.ScipyParticle` + object that defines custom particle + :param start_field: Field for initialising particles stochastically (horizontally) according to the presented density field. + :param size: Initial size of particle set + :param mode: Type of random sampling. Currently only 'monte_carlo' is implemented + :param depth: Optional list of initial depth values for particles. Default is 0m + :param time: Optional start time value for particles. Default is fieldset.U.time[0] + :param repeatdt: Optional interval (in seconds) on which to repeat the release of the ParticleSet + :param lonlatdepth_dtype: Floating precision for lon, lat, depth particle coordinates. + It is either np.float32 or np.float64. Default is np.float32 if fieldset.U.interp_method is 'linear' + and np.float64 if the interpolation method is 'cgrid_velocity' + """ + + if mode == 'monte_carlo': + data = start_field.data if isinstance(start_field.data, np.ndarray) else np.array(start_field.data) + if start_field.interp_method == 'cgrid_tracer': + p_interior = np.squeeze(data[0, 1:, 1:]) + else: # if A-grid + d = data + p_interior = (d[0, :-1, :-1] + d[0, 1:, :-1] + d[0, :-1, 1:] + d[0, 1:, 1:])/4. + p_interior = np.where(d[0, :-1, :-1] == 0, 0, p_interior) + p_interior = np.where(d[0, 1:, :-1] == 0, 0, p_interior) + p_interior = np.where(d[0, 1:, 1:] == 0, 0, p_interior) + p_interior = np.where(d[0, :-1, 1:] == 0, 0, p_interior) + p = np.reshape(p_interior, (1, p_interior.size)) + inds = np.random.choice(p_interior.size, size, replace=True, p=p[0] / np.sum(p)) + xsi = np.random.uniform(size=len(inds)) + eta = np.random.uniform(size=len(inds)) + j, i = np.unravel_index(inds, p_interior.shape) + grid = start_field.grid + if grid.gtype in [GridCode.RectilinearZGrid, GridCode.RectilinearSGrid]: + lon = grid.lon[i] + xsi * (grid.lon[i + 1] - grid.lon[i]) + lat = grid.lat[j] + eta * (grid.lat[j + 1] - grid.lat[j]) + else: + lons = np.array([grid.lon[j, i], grid.lon[j, i+1], grid.lon[j+1, i+1], grid.lon[j+1, i]]) + if grid.mesh == 'spherical': + lons[1:] = np.where(lons[1:] - lons[0] > 180, lons[1:]-360, lons[1:]) + lons[1:] = np.where(-lons[1:] + lons[0] > 180, lons[1:]+360, lons[1:]) + lon = (1-xsi)*(1-eta) * lons[0] +\ + xsi*(1-eta) * lons[1] +\ + xsi*eta * lons[2] +\ + (1-xsi)*eta * lons[3] + lat = (1-xsi)*(1-eta) * grid.lat[j, i] +\ + xsi*(1-eta) * grid.lat[j, i+1] +\ + xsi*eta * grid.lat[j+1, i+1] +\ + (1-xsi)*eta * grid.lat[j+1, i] + else: + raise NotImplementedError('Mode %s not implemented. Please use "monte carlo" algorithm instead.' % mode) + + return cls(fieldset=fieldset, pclass=pclass, lon=lon, lat=lat, depth=depth, time=time, lonlatdepth_dtype=lonlatdepth_dtype, repeatdt=repeatdt) + + @classmethod + def from_particlefile(cls, fieldset, pclass, filename, restart=True, repeatdt=None, lonlatdepth_dtype=None): + """Initialise the ParticleSet from a netcdf ParticleFile. + This creates a new ParticleSet based on the last locations and time of all particles + in the netcdf ParticleFile. Particle IDs are preserved if restart=True + + :param fieldset: :mod:`parcels.fieldset.FieldSet` object from which to sample velocity + :param pclass: mod:`parcels.particle.JITParticle` or :mod:`parcels.particle.ScipyParticle` + object that defines custom particle + :param filename: Name of the particlefile from which to read initial conditions + :param restart: Boolean to signal if pset is used for a restart (default is True). + In that case, Particle IDs are preserved. + :param repeatdt: Optional interval (in seconds) on which to repeat the release of the ParticleSet + :param lonlatdepth_dtype: Floating precision for lon, lat, depth particle coordinates. + It is either np.float32 or np.float64. Default is np.float32 if fieldset.U.interp_method is 'linear' + and np.float64 if the interpolation method is 'cgrid_velocity' + """ + + pfile = xr.open_dataset(str(filename), decode_cf=True) + + lon = np.ma.filled(pfile.variables['lon'][:, -1], np.nan) + lat = np.ma.filled(pfile.variables['lat'][:, -1], np.nan) + depth = np.ma.filled(pfile.variables['z'][:, -1], np.nan) + time = np.ma.filled(pfile.variables['time'][:, -1], np.nan) + pid = np.ma.filled(pfile.variables['trajectory'][:, -1], np.nan) + if isinstance(time[0], np.timedelta64): + time = np.array([t/np.timedelta64(1, 's') for t in time]) + + inds = np.where(np.isfinite(lon))[0] + lon = lon[inds] + lat = lat[inds] + depth = depth[inds] + time = time[inds] + pid = pid[inds] if restart else None + + return cls(fieldset=fieldset, pclass=pclass, lon=lon, lat=lat, depth=depth, time=time, + pid_orig=pid, lonlatdepth_dtype=lonlatdepth_dtype, repeatdt=repeatdt) + + @staticmethod + def lonlatdepth_dtype_from_field_interp_method(field): + if type(field) in [SummedField, NestedField]: + for f in field: + if f.interp_method == 'cgrid_velocity': + return np.float64 + else: + if field.interp_method == 'cgrid_velocity': + return np.float64 + return np.float32 + + @property + def size(self): + nparticles = 0 + for bracket_index in range(len(self._pid_mapping_bounds)): + nparticles += self._pid_mapping_bounds[bracket_index][2] + return nparticles + # return self.particles.size + + def __repr__(self): + return "\n".join([str(p) for sublist in self._plist for p in sublist]) + + def __len__(self): + # return len(self._plist) + return self.size + + def __getitem__(self, key): + bracket_index = 0 + start_index = 0 + end_index = self._pid_mapping_bounds[bracket_index][2] + # while start_index= key: + while key >= end_index: + start_index = end_index + end_index = start_index+self._pid_mapping_bounds[bracket_index][2] + bracket_index += 1 + self.next() + slot_index = key - start_index + value = self.current_p_node[slot_index] + self.reset_current_node() + return value + # return self._plist[bracket_index][slot_index] + # return self.particles[key] + + def __setitem__(self, key, value): + assert isinstance(key, tuple) and len(key) == 2, ("Error: trying to set/address data item without (bracket_index, slot_index) key.") + # self._plist[key[0]][key[1]] = value + for i in range(0, key[0]): + self.next() + self.current_p_node[key[1]] = value + self.reset_current_node() + # def __setitem__(self, key, value): + # self.particles[key] = value + + def __iadd__(self, particles): + self.add(particles) + return self + + def _create_progressbar_(self, starttime, endtime): + pbar = None + try: + pbar = progressbar.ProgressBar(max_value=abs(endtime - starttime)).start() + except: # for old versions of progressbar + try: + pbar = progressbar.ProgressBar(maxvalue=abs(endtime - starttime)).start() + except: # for even older OR newer versions + pbar = progressbar.ProgressBar(maxval=abs(endtime - starttime)).start() + return pbar + + def add(self, particles): + """Method to add particles to the ParticleSet""" + if isinstance(particles, ParticleSet): + self._plist += particles.particles_a + self._pid_mapping_bounds += particles.pid_mapping_bounds + + if self.ptype.uses_jit: + for sublist in particles.particles_c: + self._plist_c.append(sublist) + else: + raise NotImplementedError('Only ParticleSets can be added to a ParticleSet') + + # if isinstance(particles, ParticleSet): + # particles = particles.particles + # else: + # raise NotImplementedError('Only ParticleSets can be added to a ParticleSet') + # self.particles = np.append(self.particles, particles) + # if self.ptype.uses_jit: + # particles_data = [p._cptr for p in particles] + # self._particle_data = np.append(self._particle_data, particles_data) + # # Update C-pointer on particles + # for p, pdata in zip(self.particles, self._particle_data): + # p._cptr = pdata + + def _merge_brackets_(self): + lw_bound_nlist = int(self.nlist_limit / 2) + nmerges = 1 + while nmerges > 0: + nmerges = 0 + trg_bracket = None + trg_bracket_c = None + trg_bracket_index = -1 + src_bracket = None + src_bracket_c = None + src_bracket_index = -1 + for bracket_index in range(len(self._plist)): + if self._pid_mapping_bounds[bracket_index][2] < lw_bound_nlist: + if trg_bracket is None: + trg_bracket = self._plist[bracket_index] + if self._plist_c is not None and self.ptype.uses_jit: + trg_bracket_c = self._plist_c[bracket_index] + trg_bracket_index = bracket_index + elif src_bracket is None: + src_bracket = self._plist[bracket_index] + if self._plist_c is not None and self.ptype.uses_jit: + src_bracket_c = self._plist_c[bracket_index] + src_bracket_index = bracket_index + if src_bracket is not None and trg_bracket is not None: + break + if src_bracket is not None and trg_bracket is not None: + for i in range(0, trg_bracket_index): + self.next() + # self._plist[trg_bracket_index] = np.append(trg_bracket, src_bracket) + self.current_p_node = np.append(trg_bracket, src_bracket) + # local_ids = np.array([p.id for p in self._plist[trg_bracket_index]]) + local_ids = np.array([p.id for p in self.current_p_node]) + self.pid_mapping_bounds[trg_bracket_index] = (np.min(local_ids), np.max(local_ids), local_ids.shape[0]) + if self.ptype.uses_jit: + # self._plist_c[trg_bracket_index] = np.append(self._plist_c[trg_bracket_index], self._plist_c[src_bracket_index]) + self.current_c_node = np.append(trg_bracket_c, src_bracket_c) + # for p, pdata in zip(self._plist[trg_bracket_index], self._plist_c[trg_bracket_index]): + for p, pdata in zip(self.current_p_node, self.current_c_node): + p._cptr = pdata + self._plist.pop(src_bracket_index) + # self._plist.remove(self._plist[src_bracket_index]) + # del self._plist[src_bracket_index] + if self.ptype.uses_jit: + self._plist_c.pop(src_bracket_index) + # self._plist_c.remove(self._plist_c[src_bracket_index]) + # del self._plist_c[src_bracket_index] + # self._pid_mapping_bounds.remove(self._pid_mapping_bounds[src_bracket_index]) + del self._pid_mapping_bounds[src_bracket_index] + nmerges += 1 + + def pop(self, indices=-1): + """ + Method to remove particles from the ParticleSet, based on their `indices`. + Returns numpy.ndarray of removed particles. + """ + lindices = indices + if isinstance(indices, np.ndarray) or not isinstance(indices, list): + ncounter = 0 + for bracket in self._pid_mapping_bounds: + ncounter += bracket[2] + if isinstance(indices, np.ndarray): + indices = indices.tolist() + if not isinstance(indices, list): + indices = [indices, ] + lindices = [] + for i in range(len(self._plist)): + lindices.append(None) + for gindex in indices: + while gindex < 0: + gindex = ncounter+gindex + (bracket_index, slot_index) = self.get_list_array_index(gindex) + if lindices[bracket_index] is None: + lindices[bracket_index] = [slot_index, ] + else: + lindices[bracket_index].append(slot_index) + assert len(lindices) == len(self._plist), ("ParticleSet.remove() - particle set length ({}) does not match index list length ({})".format( + len(self._plist), len(lindices))) + result = None + for bracket_index in range(len(lindices)): + local_indices = lindices[bracket_index] + if local_indices is None: + continue + if isinstance(local_indices, collections.Iterable): + local_indices = np.array(local_indices) + if result is None: + result = self._remove_(bracket_index, local_indices) + else: + result = np.concatenate((result, self._remove_(bracket_index, local_indices)), axis=0) + self._merge_brackets_() + if result is None: + return result + result = result.tolist() + if len(result) == 1: + return result[0] + return result + + def remove(self, indices): + """Method to remove particles from the ParticleSet, based on their `indices`""" + lindices = indices + if isinstance(indices, np.ndarray) or not isinstance(indices, list): + ncounter = 0 + for bracket in self._pid_mapping_bounds: + ncounter += bracket[2] + if isinstance(indices, np.ndarray): + indices = indices.tolist() + if not isinstance(indices, list): + indices = [indices, ] + lindices = [] + for i in range(len(self._plist)): + lindices.append(None) + for gindex in indices: + while gindex < 0: + gindex = ncounter+gindex + (bracket_index, slot_index) = self.get_list_array_index(gindex) + if lindices[bracket_index] is None: + lindices[bracket_index] = [slot_index, ] + else: + lindices[bracket_index].append(slot_index) + assert len(lindices) == len(self._plist), ("ParticleSet.remove() - particle set length ({}) does not match index list length ({})".format( + len(self._plist), len(lindices))) + + for bracket_index in range(len(lindices)): + local_indices = lindices[bracket_index] + if local_indices is None: + continue + if isinstance(local_indices, collections.Iterable): + local_indices = np.array(local_indices) + self._remove_(bracket_index, local_indices) + self._merge_brackets_() + + # if isinstance(indices, collections.Iterable): + # particles = [self.particles[i] for i in indices] + # else: + # particles = self.particles[indices] + # self.particles = np.delete(self.particles, indices) + # if self.ptype.uses_jit: + # self._particle_data = np.delete(self._particle_data, indices) + # # == Update C-pointer on particles == # + # for p, pdata in zip(self.particles, self._particle_data): + # p._cptr = pdata + # return particles + + def remove_local_particles(self, bracket_index, local_indices): + return self._remove_(bracket_index, local_indices) + + def _remove_(self, bracket_index, local_indices): + # if isinstance(indices, collections.Iterable): + # particles = [self.particles[i] for i in indices] + # else: + # particles = self.particles[indices] + if isinstance(local_indices, collections.Iterable): + local_indices = np.array(local_indices) + for i in range(0, bracket_index): + self.next() + # return_p = self._plist[bracket_index][local_indices] + return_p = self.current_p_node[local_indices] + # self._plist[bracket_index] = np.delete(self._plist[bracket_index], local_indices) + self.current_p_node = np.delete(self.current_p_node, local_indices) + if self._ptype.uses_jit: + # self._plist_c[bracket_index] = np.delete(self._plist_c[bracket_index], local_indices) + self.current_c_node = np.delete(self.current_c_node, local_indices) + # if self._plist[bracket_index] is not None and self._plist[bracket_index].shape[0] > 0: + if self.current_p_node is not None and self.current_p_node.shape[0] > 0: + # local_ids = np.array([p.id for p in self._plist[bracket_index]]) + local_ids = np.array([p.id for p in self.current_p_node]) + self.pid_mapping_bounds[bracket_index] = (np.min(local_ids), np.max(local_ids), local_ids.shape[0]) + if self._ptype.uses_jit: + # == Update C-pointer on particles == # + # -- OLD -- # + # for p, pdata in zip(self._plist[bracket_index], self._plist_c[bracket_index]): + # p._cptr = pdata + # -- END -- # + # for slot_index in range(len(self._plist[bracket_index])): + # self._plist[bracket_index][slot_index].update_cptr(self._plist_c[bracket_index][slot_index]) + for slot_index in range(len(self.current_p_node)): + self.current_p_node[slot_index].update_cptr(self.current_p_node[slot_index]) + else: + if len(self._plist) > 1: + # self._plist.remove(self._plist[bracket_index]) + self._plist.remove(self.current_p_node) + if self._ptype.uses_jit: + # self._plist_c.remove(self._plist_c[bracket_index]) + self._plist_c.remove(self.current_c_node) + self._pid_mapping_bounds.remove(self._pid_mapping_bounds[bracket_index]) + else: + self.pid_mapping_bounds[bracket_index] = (np.iinfo(np.int32).max, np.iinfo(np.int32).min, 0) + # return particles + return return_p + + def execute(self, pyfunc=AdvectionRK4, endtime=None, runtime=None, dt=1., + moviedt=None, recovery=None, output_file=None, movie_background_field=None, + verbose_progress=None, postIterationCallbacks=None, callbackdt=None): + """Execute a given kernel function over the particle set for + multiple timesteps. Optionally also provide sub-timestepping + for particle output. + + :param pyfunc: Kernel function to execute. This can be the name of a + defined Python function or a :class:`parcels.kernel.Kernel` object. + Kernels can be concatenated using the + operator + :param endtime: End time for the timestepping loop. + It is either a datetime object or a positive double. + :param runtime: Length of the timestepping loop. Use instead of endtime. + It is either a timedelta object or a positive double. + :param dt: Timestep interval to be passed to the kernel. + It is either a timedelta object or a double. + Use a negative value for a backward-in-time simulation. + :param moviedt: Interval for inner sub-timestepping (leap), which dictates + the update frequency of animation. + It is either a timedelta object or a positive double. + None value means no animation. + :param output_file: :mod:`parcels.particlefile.ParticleFile` object for particle output + :param recovery: Dictionary with additional `:mod:parcels.tools.error` + recovery kernels to allow custom recovery behaviour in case of + kernel errors. + :param movie_background_field: field plotted as background in the movie if moviedt is set. + 'vector' shows the velocity as a vector field. + :param verbose_progress: Boolean for providing a progress bar for the kernel execution loop. + :param postIterationCallbacks: (Optional) Array of functions that are to be called after each iteration (post-process, non-Kernel) + :param callbackdt: (Optional, in conjecture with 'postIterationCallbacks) timestep inverval to (latestly) interrupt the running kernel and invoke post-iteration callbacks from 'postIterationCallbacks' + """ + + # check if pyfunc has changed since last compile. If so, recompile + if self.kernel is None or (self.kernel.pyfunc is not pyfunc and self.kernel is not pyfunc): + # Generate and store Kernel + if isinstance(pyfunc, Kernel): + self.kernel = pyfunc + else: + self.kernel = self.Kernel(pyfunc) + # Prepare JIT kernel execution + if self._ptype.uses_jit: + self.kernel.remove_lib() + cppargs = ['-DDOUBLE_COORD_VARIABLES'] if self.lonlatdepth_dtype == np.float64 else None + self.kernel.compile(compiler=GNUCompiler(cppargs=cppargs)) + self.kernel.load_lib() + + # Convert all time variables to seconds + if isinstance(endtime, delta): + raise RuntimeError('endtime must be either a datetime or a double') + if isinstance(endtime, datetime): + endtime = np.datetime64(endtime) + if isinstance(endtime, np.datetime64): + if self.time_origin.calendar is None: + raise NotImplementedError('If fieldset.time_origin is not a date, execution endtime must be a double') + endtime = self.time_origin.reltime(endtime) + if isinstance(runtime, delta): + runtime = runtime.total_seconds() + if isinstance(dt, delta): + dt = dt.total_seconds() + outputdt = output_file.outputdt if output_file else np.infty + if isinstance(outputdt, delta): + outputdt = outputdt.total_seconds() + if isinstance(moviedt, delta): + moviedt = moviedt.total_seconds() + if isinstance(callbackdt, delta): + callbackdt = callbackdt.total_seconds() + + assert runtime is None or runtime >= 0, 'runtime must be positive' + assert outputdt is None or outputdt >= 0, 'outputdt must be positive' + assert moviedt is None or moviedt >= 0, 'moviedt must be positive' + + # Set particle.time defaults based on sign of dt, if not set at ParticleSet construction + for sublist in self._plist: + for p in sublist: + if np.isnan(p.time): + mintime, maxtime = self.fieldset.gridset.dimrange('time_full') + p.time = mintime if dt >= 0 else maxtime + + # Derive _starttime and endtime from arguments or fieldset defaults + if runtime is not None and endtime is not None: + raise RuntimeError('Only one of (endtime, runtime) can be specified') + # _starttime = min([p.time for p in self]) if dt >= 0 else max([p.time for p in self]) + _starttime = min([p.time for sublist in self._plist for p in sublist]) if dt >= 0 else max([p.time for sublist in self._plist for p in sublist]) + if self.repeatdt is not None and self.repeat_starttime is None: + self.repeat_starttime = _starttime + if runtime is not None: + endtime = _starttime + runtime * np.sign(dt) + elif endtime is None: + mintime, maxtime = self.fieldset.gridset.dimrange('time_full') + endtime = maxtime if dt >= 0 else mintime + + execute_once = False + if abs(endtime-_starttime) < 1e-5 or dt == 0 or runtime == 0: + dt = 0 + runtime = 0 + endtime = _starttime + logger.warning_once("dt or runtime are zero, or endtime is equal to Particle.time. " + "The kernels will be executed once, without incrementing time") + execute_once = True + + # Initialise particle timestepping + for sublist in self._plist: + for p in sublist: + p.dt = dt + + # First write output_file, because particles could have been added + if output_file: + output_file.write(self, _starttime) + if moviedt: + self.show(field=movie_background_field, show_time=_starttime, animation=True) + + if moviedt is None: + moviedt = np.infty + if callbackdt is None: + interupt_dts = [np.infty, moviedt, outputdt] + if self.repeatdt is not None: + interupt_dts.append(self.repeatdt) + callbackdt = np.min(np.array(interupt_dts)) + time = _starttime + if self.repeatdt: + next_prelease = self.repeat_starttime + (abs(time - self.repeat_starttime) // self.repeatdt + 1) * self.repeatdt * np.sign(dt) + else: + next_prelease = np.infty if dt > 0 else - np.infty + next_output = time + outputdt if dt > 0 else time - outputdt + next_movie = time + moviedt if dt > 0 else time - moviedt + next_callback = time + callbackdt if dt > 0 else time - callbackdt + next_input = self.fieldset.computeTimeChunk(time, np.sign(dt)) + + tol = 1e-12 + if verbose_progress is None: + walltime_start = time_module.time() + if verbose_progress: + pbar = self._create_progressbar_(_starttime, endtime) + while (time < endtime and dt > 0) or (time > endtime and dt < 0) or dt == 0: + pbar = None + if verbose_progress is None and time_module.time() - walltime_start > 10: + # Showing progressbar if runtime > 10 seconds + if output_file: + logger.info('Temporary output files are stored in %s.' % output_file.tempwritedir_base) + logger.info('You can use "parcels_convert_npydir_to_netcdf %s" to convert these ' + 'to a NetCDF file during the run.' % output_file.tempwritedir_base) + pbar = self._create_progressbar_(_starttime, endtime) + verbose_progress = True + if dt > 0: + time = min(next_prelease, next_input, next_output, next_movie, next_callback, endtime) + else: + time = max(next_prelease, next_input, next_output, next_movie, next_callback, endtime) + self.kernel.execute(self, endtime=time, dt=dt, recovery=recovery, output_file=output_file, execute_once=execute_once) + logger.info("Computed time = {}".format(time)) + if abs(time-next_prelease) < tol: + pset_new = ParticleSet(fieldset=self.fieldset, time=time, lon=self.repeatlon, + lat=self.repeatlat, depth=self.repeatdepth, + pclass=self.repeatpclass, lonlatdepth_dtype=self.lonlatdepth_dtype, + partitions=False, pid_orig=self.repeatpid, **self.repeatkwargs) + for sublist in pset_new.particles_a: + for p in sublist: + p.dt = dt + self.add(pset_new) + next_prelease += self.repeatdt * np.sign(dt) + if abs(time-next_output) < tol: + if output_file: + output_file.write(self, time) + next_output += outputdt * np.sign(dt) + if abs(time-next_movie) < tol: + self.show(field=movie_background_field, show_time=time, animation=True) + next_movie += moviedt * np.sign(dt) + # ==== insert post-process here to also allow for memory clean-up via external func ==== # + if abs(time-next_callback) < tol: + if postIterationCallbacks is not None: + for extFunc in postIterationCallbacks: + extFunc() + next_callback += callbackdt * np.sign(dt) + if time != endtime: + next_input = self.fieldset.computeTimeChunk(time, dt) + if dt == 0: + break + if verbose_progress: + pbar.update(abs(time - _starttime)) + + if output_file: + output_file.write(self, time) + if verbose_progress: + pbar.finish() + + def show(self, with_particles=True, show_time=None, field=None, domain=None, projection=None, + land=True, vmin=None, vmax=None, savefile=None, animation=False, **kwargs): + """Method to 'show' a Parcels ParticleSet + + :param with_particles: Boolean whether to show particles + :param show_time: Time at which to show the ParticleSet + :param field: Field to plot under particles (either None, a Field object, or 'vector') + :param domain: dictionary (with keys 'N', 'S', 'E', 'W') defining domain to show + :param projection: type of cartopy projection to use (default PlateCarree) + :param land: Boolean whether to show land. This is ignored for flat meshes + :param vmin: minimum colour scale (only in single-plot mode) + :param vmax: maximum colour scale (only in single-plot mode) + :param savefile: Name of a file to save the plot to + :param animation: Boolean whether result is a single plot, or an animation + """ + from parcels.plotting import plotparticles + plotparticles(particles=self, with_particles=with_particles, show_time=show_time, field=field, domain=domain, + projection=projection, land=land, vmin=vmin, vmax=vmax, savefile=savefile, animation=animation, **kwargs) + + def density(self, field=None, particle_val=None, relative=False, area_scale=False): + """Method to calculate the density of particles in a ParticleSet from their locations, + through a 2D histogram. + + :param field: Optional :mod:`parcels.field.Field` object to calculate the histogram + on. Default is `fieldset.U` + :param particle_val: Optional numpy-array of values to weigh each particle with, + or string name of particle variable to use weigh particles with. + Default is None, resulting in a value of 1 for each particle + :param relative: Boolean to control whether the density is scaled by the total + weight of all particles. Default is False + :param area_scale: Boolean to control whether the density is scaled by the area + (in m^2) of each grid cell. Default is False + """ + + field = field if field else self.fieldset.U + if isinstance(particle_val, str): + particle_val = [np.array([getattr(p, particle_val) for p in sublist]) for sublist in self._plist] + else: + particle_val = particle_val if particle_val else np.ones(self.size) + particle_val = self.covert_static_array_to_list_array(particle_val) + density = np.zeros((field.grid.lat.size, field.grid.lon.size), dtype=np.float32) +# =================================================================================================================== # + for bracket_index, subfield in enumerate(self._plist): + for slot_index, p in enumerate(subfield): + try: # breaks if either p.xi, p.yi, p.zi, p.ti do not exist (in scipy) or field not in fieldset + if p.ti[field.igrid] < 0: # xi, yi, zi, ti, not initialised + raise('error') + xi = p.xi[field.igrid] + yi = p.yi[field.igrid] + except: + _, _, _, xi, yi, _ = field.search_indices(p.lon, p.lat, p.depth, 0, 0, search2D=True) + density[yi, xi] += particle_val[bracket_index][slot_index] + + # for pi, p in enumerate(self.particles): + # try: # breaks if either p.xi, p.yi, p.zi, p.ti do not exist (in scipy) or field not in fieldset + # if p.ti[field.igrid] < 0: # xi, yi, zi, ti, not initialised + # raise('error') + # xi = p.xi[field.igrid] + # yi = p.yi[field.igrid] + # except: + # _, _, _, xi, yi, _ = field.search_indices(p.lon, p.lat, p.depth, 0, 0, search2D=True) + # density[yi, xi] += particle_val[pi] + + if relative: + psum = 0 + for bracket_index, subfield in enumerate(particle_val): + psum += np.sum(particle_val[bracket_index]) + density /= psum + # density /= np.sum(particle_val) + + if area_scale: + density /= field.cell_areas() + + return density + + def Kernel(self, pyfunc, c_include="", delete_cfiles=True): + """Wrapper method to convert a `pyfunc` into a :class:`parcels.kernel.Kernel` object + based on `fieldset` and `ptype` of the ParticleSet + :param delete_cfiles: Boolean whether to delete the C-files after compilation in JIT mode (default is True) + """ + return Kernel(self.fieldset, self.ptype, pyfunc=pyfunc, c_include=c_include, + delete_cfiles=delete_cfiles) + + def ParticleFile(self, *args, **kwargs): + """Wrapper method to initialise a :class:`parcels.particlefile.ParticleFile` + object from the ParticleSet""" + return ParticleFile(*args, particleset=self, **kwargs) From 7a626a377df84645b172e4fa45bf0477a9dfd67d Mon Sep 17 00:00:00 2001 From: Christian Kehl Date: Wed, 9 Dec 2020 21:49:27 +0100 Subject: [PATCH 16/56] list_of_pset_array_trials - finalised the correct memory consumption measurement --- parcels/particleset_benchmark.py | 8 +++++--- parcels/tools/performance_logger.py | 15 ++++++++++----- 2 files changed, 15 insertions(+), 8 deletions(-) diff --git a/parcels/particleset_benchmark.py b/parcels/particleset_benchmark.py index 6e5c53c504..24face7e5c 100644 --- a/parcels/particleset_benchmark.py +++ b/parcels/particleset_benchmark.py @@ -319,7 +319,7 @@ def execute(self, pyfunc=AdvectionRK4, endtime=None, runtime=None, dt=1., # else: # mem_B_used_total = self.process.memory_info().rss # mem_B_used_total = self.process.memory_info().rss - mem_B_used_total = 0 + # mem_B_used_total = 0 if USE_RUSE_SYNC_MEMLOG: mem_B_used_total = measure_mem_usage() else: @@ -350,7 +350,7 @@ def execute(self, pyfunc=AdvectionRK4, endtime=None, runtime=None, dt=1., # self.nparticle_log.advance_iteration(self.size) # self.compute_log.advance_iteration() # self.io_log.advance_iteration() - # self.mem_log.advance_iteration(self.process.memory_info().rss) + # self.mem_log.advance_iteration() # self.mem_io_log.advance_iteration() # self.plot_log.advance_iteration() # self.total_log.advance_iteration() @@ -419,6 +419,8 @@ def plot_and_log(self, total_times = None, compute_times = None, io_times = None do_drawt_plot = False do_mem_plot = True do_mem_plot_async = True + if not USE_ASYNC_MEMLOG: + do_mem_plot_async = False do_npart_plot = True assert (len(plot_t) == len(plot_ct)) if len(plot_t) != len(plot_iot): @@ -480,7 +482,7 @@ def plot_and_log(self, total_times = None, compute_times = None, io_times = None max_mem_sync = 0 if memory_used is not None and len(memory_used) > 1: memory_used = np.floor(memory_used / (1024*1024)) - memory_used = memory_used.astype(dtype=np.uint32) + memory_used = memory_used.astype(dtype=np.int64) max_mem_sync = memory_used.max() max_mem_async = 0 if USE_ASYNC_MEMLOG: diff --git a/parcels/tools/performance_logger.py b/parcels/tools/performance_logger.py index cd13634545..5fadb01a1b 100644 --- a/parcels/tools/performance_logger.py +++ b/parcels/tools/performance_logger.py @@ -248,19 +248,21 @@ def async_run(self): if self.differential_measurement: self.async_run_diff_measurement() else: - pass + self.async_run_measurement() def async_run_diff_measurement(self): if self._measure_start_value is None: self._measure_start_value = self._measure_func() self._measure_partial_values.append(0) while not self._event.is_set(): - self._measure_partial_values.append( self._measure_func()-self._measure_start_value ) + value = self._measure_func()-self._measure_start_value + self._measure_partial_values.append( value ) sleep(self._measure_interval) def async_run_measurement(self): while not self._event.is_set(): - self._measure_partial_values.append( self.measure_func() ) + value = self.measure_func() + self._measure_partial_values.append( value ) sleep(self.measure_interval) def start_partial_measurement(self): @@ -270,7 +272,7 @@ def start_partial_measurement(self): del self._measure_partial_values[:] self._measure_partial_values = [] self._event = Event() - self._thread = Thread(target=self.async_run_measurement) + self._thread = Thread(target=self.async_run) self._thread.start() def stop_partial_measurement(self): @@ -284,9 +286,12 @@ def stop_partial_measurement(self): sleep(self._measure_interval) del self._thread self._thread = None + self._event.clear() + del self._event + self._event = None self._measure_start_value = None # param_partial_mean = np.array(self._measure_partial_values).mean() - param_partial_mean = np.array(self._measure_partial_values).max() + param_partial_mean = np.array(self._measure_partial_values).max(initial=0) self.advance_iteration(param_partial_mean) def get_params(self): From 2f165b622442e3424f9b01fe97439e8c7ae0bcfb Mon Sep 17 00:00:00 2001 From: Christian Kehl Date: Tue, 7 Jul 2020 20:02:48 +0200 Subject: [PATCH 17/56] list_of_pset_array_trials - initial addenda for list-of-particles version --- parcels/particleset.py | 1 - 1 file changed, 1 deletion(-) diff --git a/parcels/particleset.py b/parcels/particleset.py index bba5488c6d..fcc9288ebd 100644 --- a/parcels/particleset.py +++ b/parcels/particleset.py @@ -30,7 +30,6 @@ __all__ = ['ParticleSet'] - class ParticleSet(object): """Container class for storing particle and executing kernel over them. From cdf4a52c08b80e9f210127cd8295d6658b95326d Mon Sep 17 00:00:00 2001 From: Christian Kehl Date: Wed, 8 Jul 2020 16:10:52 +0200 Subject: [PATCH 18/56] list_of_pset_array_trials - wrapped up list-of-array adaptations in the ParticleSet class. --- parcels/particleset.py | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/parcels/particleset.py b/parcels/particleset.py index fcc9288ebd..47cf6f05e8 100644 --- a/parcels/particleset.py +++ b/parcels/particleset.py @@ -276,9 +276,11 @@ def get_list_array_index(self, pdata): return (bracket_index, slot_index) elif isinstance(pdata, self.pclass): # == parameter is a particle itself == # + # nparticle = 0 bracket_index = 0 bracket_info = self.pid_mapping_bounds[bracket_index] while bracket_index < len(self._plist) and ((pdata.id < bracket_info[0]) or (pdata.id > bracket_info[1])): + # nparticle += bracket_info[2] bracket_index += 1 bracket_info = self.pid_mapping_bounds[bracket_index] if bracket_index >= len(self._plist): @@ -287,9 +289,11 @@ def get_list_array_index(self, pdata): slot_index = 0 while (slot_index < self._plist[bracket_index].shape[0]) and (self._plist[bracket_index][slot_index].id != pdata.id): slot_index += 1 + # nparticle += 1 if slot_index >= self._plist[bracket_index].shape[0]: logger.warn_once("ParticleSet.get_cptr_index() - requested particle ({}) not found.".format(pdata)) return -1 + # return nparticle return (bracket_index, slot_index) def get_list_array_index_by_PID(self, pdata): From 084c6e23d82649edc255183c54e922277237b243 Mon Sep 17 00:00:00 2001 From: Christian Kehl Date: Tue, 14 Jul 2020 18:20:29 +0200 Subject: [PATCH 19/56] list_of_pset_array_trials - fixing linter errors --- parcels/particleset.py | 3 +- parcels/particleset_benchmark.py | 77 +++++++++++++++++++++++++++++++- 2 files changed, 77 insertions(+), 3 deletions(-) diff --git a/parcels/particleset.py b/parcels/particleset.py index 47cf6f05e8..f1787c1f21 100644 --- a/parcels/particleset.py +++ b/parcels/particleset.py @@ -30,6 +30,7 @@ __all__ = ['ParticleSet'] + class ParticleSet(object): """Container class for storing particle and executing kernel over them. @@ -276,7 +277,6 @@ def get_list_array_index(self, pdata): return (bracket_index, slot_index) elif isinstance(pdata, self.pclass): # == parameter is a particle itself == # - # nparticle = 0 bracket_index = 0 bracket_info = self.pid_mapping_bounds[bracket_index] while bracket_index < len(self._plist) and ((pdata.id < bracket_info[0]) or (pdata.id > bracket_info[1])): @@ -293,7 +293,6 @@ def get_list_array_index(self, pdata): if slot_index >= self._plist[bracket_index].shape[0]: logger.warn_once("ParticleSet.get_cptr_index() - requested particle ({}) not found.".format(pdata)) return -1 - # return nparticle return (bracket_index, slot_index) def get_list_array_index_by_PID(self, pdata): diff --git a/parcels/particleset_benchmark.py b/parcels/particleset_benchmark.py index 24face7e5c..e81fdf9945 100644 --- a/parcels/particleset_benchmark.py +++ b/parcels/particleset_benchmark.py @@ -51,6 +51,82 @@ def measure_mem_usage(): USE_ASYNC_MEMLOG = False USE_RUSE_SYNC_MEMLOG = False # can be faulty +class ParticleSet_TimingLog(): + stime = 0 + etime = 0 + mtime = 0 + samples = [] + times_steps = [] + _iter = 0 + + def start_timing(self): + if MPI: + mpi_comm = MPI.COMM_WORLD + mpi_rank = mpi_comm.Get_rank() + if mpi_rank == 0: + # self.stime = MPI.Wtime() + # self.stime = time_module.perf_counter() + self.stime = time_module.process_time() + else: + self.stime = time_module.perf_counter() + + def stop_timing(self): + if MPI: + mpi_comm = MPI.COMM_WORLD + mpi_rank = mpi_comm.Get_rank() + if mpi_rank == 0: + # self.etime = MPI.Wtime() + # self.etime = time_module.perf_counter() + self.etime = time_module.process_time() + else: + self.etime = time_module.perf_counter() + + def accumulate_timing(self): + if MPI: + mpi_comm = MPI.COMM_WORLD + mpi_rank = mpi_comm.Get_rank() + if mpi_rank == 0: + self.mtime += (self.etime-self.stime) + else: + self.mtime = 0 + else: + self.mtime += (self.etime-self.stime) + + def advance_iteration(self): + if MPI: + mpi_comm = MPI.COMM_WORLD + mpi_rank = mpi_comm.Get_rank() + if mpi_rank == 0: + self.times_steps.append(self.mtime) + self.samples.append(self._iter) + self._iter += 1 + self.mtime = 0 + else: + self.times_steps.append(self.mtime) + self.samples.append(self._iter) + self._iter += 1 + self.mtime = 0 + + +class ParticleSet_ParamLogging(): + samples = [] + params = [] + _iter = 0 + + def advance_iteration(self, param): + if MPI: + mpi_comm = MPI.COMM_WORLD + mpi_rank = mpi_comm.Get_rank() + if mpi_rank == 0: + self.params.append(param) + self.samples.append(self._iter) + self._iter += 1 + else: + self.params.append(param) + self.samples.append(self._iter) + self._iter += 1 + + class ParticleSet_Benchmark(ParticleSet): def __init__(self, fieldset, pclass=JITParticle, lon=None, lat=None, depth=None, time=None, repeatdt=None, @@ -69,7 +145,6 @@ def __init__(self, fieldset, pclass=JITParticle, lon=None, lat=None, depth=None, def set_async_memlog_interval(self, interval): self.async_mem_log.measure_interval = interval - # @profile def execute(self, pyfunc=AdvectionRK4, endtime=None, runtime=None, dt=1., moviedt=None, recovery=None, output_file=None, movie_background_field=None, verbose_progress=None, postIterationCallbacks=None, callbackdt=None): From d16560b7ca10abcb967fe28ef7d58b3b4b85c74a Mon Sep 17 00:00:00 2001 From: Christian Kehl Date: Thu, 16 Jul 2020 17:07:21 +0200 Subject: [PATCH 20/56] list_of_pset_array_trials - added benchmark refactoring. tested and working (stommel). --- parcels/particleset_benchmark.py | 71 ++++++++++++++++++++++++++------ performance/benchmark_stommel.py | 3 -- 2 files changed, 59 insertions(+), 15 deletions(-) diff --git a/parcels/particleset_benchmark.py b/parcels/particleset_benchmark.py index e81fdf9945..c96f4de9a3 100644 --- a/parcels/particleset_benchmark.py +++ b/parcels/particleset_benchmark.py @@ -55,10 +55,35 @@ class ParticleSet_TimingLog(): stime = 0 etime = 0 mtime = 0 - samples = [] - times_steps = [] + _samples = [] + _times_steps = [] _iter = 0 + def __init__(self): + self.stime = 0 + self.etime = 0 + self.mtime = 0 + self._samples = [] + self._times_steps = [] + self._iter = 0 + + @property + def timing(self): + return self._times_steps + + @property + def samples(self): + return self._samples + + def __len__(self): + return len(self._samples) + + def get_values(self): + return self._times_steps + + def get_value(self, index): + return self._times_steps[index] + def start_timing(self): if MPI: mpi_comm = MPI.COMM_WORLD @@ -97,33 +122,55 @@ def advance_iteration(self): mpi_comm = MPI.COMM_WORLD mpi_rank = mpi_comm.Get_rank() if mpi_rank == 0: - self.times_steps.append(self.mtime) - self.samples.append(self._iter) + self._times_steps.append(self.mtime) + self._samples.append(self._iter) self._iter += 1 self.mtime = 0 else: - self.times_steps.append(self.mtime) - self.samples.append(self._iter) + self._times_steps.append(self.mtime) + self._samples.append(self._iter) self._iter += 1 self.mtime = 0 class ParticleSet_ParamLogging(): - samples = [] - params = [] + _samples = [] + _params = [] _iter = 0 + def __init__(self): + self._samples = [] + self._params = [] + self._iter = 0 + + @property + def samples(self): + return self._samples + + @property + def params(self): + return self._params + + def get_params(self): + return self._params + + def get_param(self, index): + return self._params[index] + + def __len__(self): + return len(self._samples) + def advance_iteration(self, param): if MPI: mpi_comm = MPI.COMM_WORLD mpi_rank = mpi_comm.Get_rank() if mpi_rank == 0: - self.params.append(param) - self.samples.append(self._iter) + self._params.append(param) + self._samples.append(self._iter) self._iter += 1 else: - self.params.append(param) - self.samples.append(self._iter) + self._params.append(param) + self._samples.append(self._iter) self._iter += 1 diff --git a/performance/benchmark_stommel.py b/performance/benchmark_stommel.py index 6bf63192c4..eb0e125b64 100644 --- a/performance/benchmark_stommel.py +++ b/performance/benchmark_stommel.py @@ -426,11 +426,8 @@ def Age(particle, fieldset, time): mpi_comm = MPI.COMM_WORLD mpi_rank = mpi_comm.Get_rank() if mpi_rank==0: - # endtime = ostime.time() - # endtime = MPI.Wtime() endtime = ostime.process_time() else: - #endtime = ostime.time() endtime = ostime.process_time() if args.write_out: From 2bb9ad5f962d717cd9a1c3702e224962a9add89f Mon Sep 17 00:00:00 2001 From: Christian Kehl Date: Mon, 20 Jul 2020 23:06:58 +0200 Subject: [PATCH 21/56] list_of_pset_array_trials - rebase to benchmarking branch rev 54e4c18 --- parcels/particleset_benchmark.py | 123 ------------------------------- performance/benchmark_stommel.py | 3 + 2 files changed, 3 insertions(+), 123 deletions(-) diff --git a/parcels/particleset_benchmark.py b/parcels/particleset_benchmark.py index c96f4de9a3..3cc9288151 100644 --- a/parcels/particleset_benchmark.py +++ b/parcels/particleset_benchmark.py @@ -51,129 +51,6 @@ def measure_mem_usage(): USE_ASYNC_MEMLOG = False USE_RUSE_SYNC_MEMLOG = False # can be faulty -class ParticleSet_TimingLog(): - stime = 0 - etime = 0 - mtime = 0 - _samples = [] - _times_steps = [] - _iter = 0 - - def __init__(self): - self.stime = 0 - self.etime = 0 - self.mtime = 0 - self._samples = [] - self._times_steps = [] - self._iter = 0 - - @property - def timing(self): - return self._times_steps - - @property - def samples(self): - return self._samples - - def __len__(self): - return len(self._samples) - - def get_values(self): - return self._times_steps - - def get_value(self, index): - return self._times_steps[index] - - def start_timing(self): - if MPI: - mpi_comm = MPI.COMM_WORLD - mpi_rank = mpi_comm.Get_rank() - if mpi_rank == 0: - # self.stime = MPI.Wtime() - # self.stime = time_module.perf_counter() - self.stime = time_module.process_time() - else: - self.stime = time_module.perf_counter() - - def stop_timing(self): - if MPI: - mpi_comm = MPI.COMM_WORLD - mpi_rank = mpi_comm.Get_rank() - if mpi_rank == 0: - # self.etime = MPI.Wtime() - # self.etime = time_module.perf_counter() - self.etime = time_module.process_time() - else: - self.etime = time_module.perf_counter() - - def accumulate_timing(self): - if MPI: - mpi_comm = MPI.COMM_WORLD - mpi_rank = mpi_comm.Get_rank() - if mpi_rank == 0: - self.mtime += (self.etime-self.stime) - else: - self.mtime = 0 - else: - self.mtime += (self.etime-self.stime) - - def advance_iteration(self): - if MPI: - mpi_comm = MPI.COMM_WORLD - mpi_rank = mpi_comm.Get_rank() - if mpi_rank == 0: - self._times_steps.append(self.mtime) - self._samples.append(self._iter) - self._iter += 1 - self.mtime = 0 - else: - self._times_steps.append(self.mtime) - self._samples.append(self._iter) - self._iter += 1 - self.mtime = 0 - - -class ParticleSet_ParamLogging(): - _samples = [] - _params = [] - _iter = 0 - - def __init__(self): - self._samples = [] - self._params = [] - self._iter = 0 - - @property - def samples(self): - return self._samples - - @property - def params(self): - return self._params - - def get_params(self): - return self._params - - def get_param(self, index): - return self._params[index] - - def __len__(self): - return len(self._samples) - - def advance_iteration(self, param): - if MPI: - mpi_comm = MPI.COMM_WORLD - mpi_rank = mpi_comm.Get_rank() - if mpi_rank == 0: - self._params.append(param) - self._samples.append(self._iter) - self._iter += 1 - else: - self._params.append(param) - self._samples.append(self._iter) - self._iter += 1 - - class ParticleSet_Benchmark(ParticleSet): def __init__(self, fieldset, pclass=JITParticle, lon=None, lat=None, depth=None, time=None, repeatdt=None, diff --git a/performance/benchmark_stommel.py b/performance/benchmark_stommel.py index eb0e125b64..6bf63192c4 100644 --- a/performance/benchmark_stommel.py +++ b/performance/benchmark_stommel.py @@ -426,8 +426,11 @@ def Age(particle, fieldset, time): mpi_comm = MPI.COMM_WORLD mpi_rank = mpi_comm.Get_rank() if mpi_rank==0: + # endtime = ostime.time() + # endtime = MPI.Wtime() endtime = ostime.process_time() else: + #endtime = ostime.time() endtime = ostime.process_time() if args.write_out: From 77d057baaf736ee3586dff45a0b861b6044d0a62 Mon Sep 17 00:00:00 2001 From: Christian Kehl Date: Wed, 9 Dec 2020 21:49:27 +0100 Subject: [PATCH 22/56] list_of_pset_array_trials - finalised the correct memory consumption measurement --- parcels/particleset_benchmark.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/parcels/particleset_benchmark.py b/parcels/particleset_benchmark.py index 3cc9288151..ece52e391b 100644 --- a/parcels/particleset_benchmark.py +++ b/parcels/particleset_benchmark.py @@ -477,7 +477,7 @@ def plot_and_log(self, total_times = None, compute_times = None, io_times = None header_string = "target_N, start_N, final_N, avg_N, ncores, avg_kt_total[s], avg_kt_compute[s], avg_kt_io[s], avg_kt_plot[s], cum_t_total[s], cum_t_compute[s], com_t_io[s], cum_t_plot[s], max_mem[MB]\n" f.write(header_string) data_string = "{}, {}, {}, {}, {}, ".format(target_N, nparticles_t0, nparticles_tN, nparticles.mean(), ncores) - data_string += "{:2.10f}, {:2.10f}, {:2.10f}, {:2.10f}, ".format(total_times.mean(), compute_times.mean(), io_times.mean(), plot_times.mean()) + data_string+= "{:2.10f}, {:2.10f}, {:2.10f}, {:2.10f}, ".format(total_times.mean(), compute_times.mean(), io_times.mean(), plot_times.mean()) max_mem_sync = 0 if memory_used is not None and len(memory_used) > 1: memory_used = np.floor(memory_used / (1024*1024)) @@ -490,6 +490,6 @@ def plot_and_log(self, total_times = None, compute_times = None, io_times = None memory_used_async = memory_used_async.astype(dtype=np.int64) max_mem_async = memory_used_async.max() max_mem = max(max_mem_sync, max_mem_async) - data_string += "{:10.4f}, {:10.4f}, {:10.4f}, {:10.4f}, {}".format(total_times.sum(), compute_times.sum(), io_times.sum(), plot_times.sum(), max_mem) + data_string+= "{:10.4f}, {:10.4f}, {:10.4f}, {:10.4f}, {}".format(total_times.sum(), compute_times.sum(), io_times.sum(), plot_times.sum(), max_mem) f.write(data_string) From b78fcd2b5c64ad39ff6812dd8eca4797666d2802 Mon Sep 17 00:00:00 2001 From: Christian Kehl Date: Tue, 7 Jul 2020 20:02:48 +0200 Subject: [PATCH 23/56] list_of_pset_array_trials - initial addenda for list-of-particles version --- parcels/particleset.py | 1 - 1 file changed, 1 deletion(-) diff --git a/parcels/particleset.py b/parcels/particleset.py index f1787c1f21..667c1d2285 100644 --- a/parcels/particleset.py +++ b/parcels/particleset.py @@ -30,7 +30,6 @@ __all__ = ['ParticleSet'] - class ParticleSet(object): """Container class for storing particle and executing kernel over them. From 44985e7fbe40afe90a2c0662e08979bc3920f93e Mon Sep 17 00:00:00 2001 From: Christian Kehl Date: Wed, 8 Jul 2020 16:10:52 +0200 Subject: [PATCH 24/56] list_of_pset_array_trials - wrapped up list-of-array adaptations in the ParticleSet class. --- parcels/particleset.py | 1 + 1 file changed, 1 insertion(+) diff --git a/parcels/particleset.py b/parcels/particleset.py index 667c1d2285..59e4b45153 100644 --- a/parcels/particleset.py +++ b/parcels/particleset.py @@ -276,6 +276,7 @@ def get_list_array_index(self, pdata): return (bracket_index, slot_index) elif isinstance(pdata, self.pclass): # == parameter is a particle itself == # + # nparticle = 0 bracket_index = 0 bracket_info = self.pid_mapping_bounds[bracket_index] while bracket_index < len(self._plist) and ((pdata.id < bracket_info[0]) or (pdata.id > bracket_info[1])): From 6c04661fd7e0a86583575ff56c6602983fbb69b0 Mon Sep 17 00:00:00 2001 From: Christian Kehl Date: Tue, 14 Jul 2020 18:20:29 +0200 Subject: [PATCH 25/56] list_of_pset_array_trials - fixing linter errors --- parcels/particleset.py | 2 +- parcels/particleset_benchmark.py | 76 ++++++++++++++++++++++++++++++++ 2 files changed, 77 insertions(+), 1 deletion(-) diff --git a/parcels/particleset.py b/parcels/particleset.py index 59e4b45153..f1787c1f21 100644 --- a/parcels/particleset.py +++ b/parcels/particleset.py @@ -30,6 +30,7 @@ __all__ = ['ParticleSet'] + class ParticleSet(object): """Container class for storing particle and executing kernel over them. @@ -276,7 +277,6 @@ def get_list_array_index(self, pdata): return (bracket_index, slot_index) elif isinstance(pdata, self.pclass): # == parameter is a particle itself == # - # nparticle = 0 bracket_index = 0 bracket_info = self.pid_mapping_bounds[bracket_index] while bracket_index < len(self._plist) and ((pdata.id < bracket_info[0]) or (pdata.id > bracket_info[1])): diff --git a/parcels/particleset_benchmark.py b/parcels/particleset_benchmark.py index ece52e391b..d71943c5ac 100644 --- a/parcels/particleset_benchmark.py +++ b/parcels/particleset_benchmark.py @@ -51,6 +51,82 @@ def measure_mem_usage(): USE_ASYNC_MEMLOG = False USE_RUSE_SYNC_MEMLOG = False # can be faulty +class ParticleSet_TimingLog(): + stime = 0 + etime = 0 + mtime = 0 + samples = [] + times_steps = [] + _iter = 0 + + def start_timing(self): + if MPI: + mpi_comm = MPI.COMM_WORLD + mpi_rank = mpi_comm.Get_rank() + if mpi_rank == 0: + # self.stime = MPI.Wtime() + # self.stime = time_module.perf_counter() + self.stime = time_module.process_time() + else: + self.stime = time_module.perf_counter() + + def stop_timing(self): + if MPI: + mpi_comm = MPI.COMM_WORLD + mpi_rank = mpi_comm.Get_rank() + if mpi_rank == 0: + # self.etime = MPI.Wtime() + # self.etime = time_module.perf_counter() + self.etime = time_module.process_time() + else: + self.etime = time_module.perf_counter() + + def accumulate_timing(self): + if MPI: + mpi_comm = MPI.COMM_WORLD + mpi_rank = mpi_comm.Get_rank() + if mpi_rank == 0: + self.mtime += (self.etime-self.stime) + else: + self.mtime = 0 + else: + self.mtime += (self.etime-self.stime) + + def advance_iteration(self): + if MPI: + mpi_comm = MPI.COMM_WORLD + mpi_rank = mpi_comm.Get_rank() + if mpi_rank == 0: + self.times_steps.append(self.mtime) + self.samples.append(self._iter) + self._iter += 1 + self.mtime = 0 + else: + self.times_steps.append(self.mtime) + self.samples.append(self._iter) + self._iter += 1 + self.mtime = 0 + + +class ParticleSet_ParamLogging(): + samples = [] + params = [] + _iter = 0 + + def advance_iteration(self, param): + if MPI: + mpi_comm = MPI.COMM_WORLD + mpi_rank = mpi_comm.Get_rank() + if mpi_rank == 0: + self.params.append(param) + self.samples.append(self._iter) + self._iter += 1 + else: + self.params.append(param) + self.samples.append(self._iter) + self._iter += 1 + + class ParticleSet_Benchmark(ParticleSet): def __init__(self, fieldset, pclass=JITParticle, lon=None, lat=None, depth=None, time=None, repeatdt=None, From 84b05f1925ffbe7218c03a5d7591ae50d8ff789b Mon Sep 17 00:00:00 2001 From: Christian Kehl Date: Thu, 16 Jul 2020 17:07:21 +0200 Subject: [PATCH 26/56] list_of_pset_array_trials - added benchmark refactoring. tested and working (stommel). --- parcels/particleset_benchmark.py | 71 ++++++++++++++++++++++++++------ performance/benchmark_stommel.py | 3 -- 2 files changed, 59 insertions(+), 15 deletions(-) diff --git a/parcels/particleset_benchmark.py b/parcels/particleset_benchmark.py index d71943c5ac..999f02a402 100644 --- a/parcels/particleset_benchmark.py +++ b/parcels/particleset_benchmark.py @@ -55,10 +55,35 @@ class ParticleSet_TimingLog(): stime = 0 etime = 0 mtime = 0 - samples = [] - times_steps = [] + _samples = [] + _times_steps = [] _iter = 0 + def __init__(self): + self.stime = 0 + self.etime = 0 + self.mtime = 0 + self._samples = [] + self._times_steps = [] + self._iter = 0 + + @property + def timing(self): + return self._times_steps + + @property + def samples(self): + return self._samples + + def __len__(self): + return len(self._samples) + + def get_values(self): + return self._times_steps + + def get_value(self, index): + return self._times_steps[index] + def start_timing(self): if MPI: mpi_comm = MPI.COMM_WORLD @@ -97,33 +122,55 @@ def advance_iteration(self): mpi_comm = MPI.COMM_WORLD mpi_rank = mpi_comm.Get_rank() if mpi_rank == 0: - self.times_steps.append(self.mtime) - self.samples.append(self._iter) + self._times_steps.append(self.mtime) + self._samples.append(self._iter) self._iter += 1 self.mtime = 0 else: - self.times_steps.append(self.mtime) - self.samples.append(self._iter) + self._times_steps.append(self.mtime) + self._samples.append(self._iter) self._iter += 1 self.mtime = 0 class ParticleSet_ParamLogging(): - samples = [] - params = [] + _samples = [] + _params = [] _iter = 0 + def __init__(self): + self._samples = [] + self._params = [] + self._iter = 0 + + @property + def samples(self): + return self._samples + + @property + def params(self): + return self._params + + def get_params(self): + return self._params + + def get_param(self, index): + return self._params[index] + + def __len__(self): + return len(self._samples) + def advance_iteration(self, param): if MPI: mpi_comm = MPI.COMM_WORLD mpi_rank = mpi_comm.Get_rank() if mpi_rank == 0: - self.params.append(param) - self.samples.append(self._iter) + self._params.append(param) + self._samples.append(self._iter) self._iter += 1 else: - self.params.append(param) - self.samples.append(self._iter) + self._params.append(param) + self._samples.append(self._iter) self._iter += 1 diff --git a/performance/benchmark_stommel.py b/performance/benchmark_stommel.py index 6bf63192c4..eb0e125b64 100644 --- a/performance/benchmark_stommel.py +++ b/performance/benchmark_stommel.py @@ -426,11 +426,8 @@ def Age(particle, fieldset, time): mpi_comm = MPI.COMM_WORLD mpi_rank = mpi_comm.Get_rank() if mpi_rank==0: - # endtime = ostime.time() - # endtime = MPI.Wtime() endtime = ostime.process_time() else: - #endtime = ostime.time() endtime = ostime.process_time() if args.write_out: From 9b3a1c32153a857a45428c8fd22ded613e4707e1 Mon Sep 17 00:00:00 2001 From: Christian Kehl Date: Mon, 20 Jul 2020 23:06:58 +0200 Subject: [PATCH 27/56] list_of_pset_array_trials - rebase to benchmarking branch rev 54e4c18 --- parcels/particleset_benchmark.py | 123 ------------------------------- performance/benchmark_stommel.py | 3 + 2 files changed, 3 insertions(+), 123 deletions(-) diff --git a/parcels/particleset_benchmark.py b/parcels/particleset_benchmark.py index 999f02a402..ece52e391b 100644 --- a/parcels/particleset_benchmark.py +++ b/parcels/particleset_benchmark.py @@ -51,129 +51,6 @@ def measure_mem_usage(): USE_ASYNC_MEMLOG = False USE_RUSE_SYNC_MEMLOG = False # can be faulty -class ParticleSet_TimingLog(): - stime = 0 - etime = 0 - mtime = 0 - _samples = [] - _times_steps = [] - _iter = 0 - - def __init__(self): - self.stime = 0 - self.etime = 0 - self.mtime = 0 - self._samples = [] - self._times_steps = [] - self._iter = 0 - - @property - def timing(self): - return self._times_steps - - @property - def samples(self): - return self._samples - - def __len__(self): - return len(self._samples) - - def get_values(self): - return self._times_steps - - def get_value(self, index): - return self._times_steps[index] - - def start_timing(self): - if MPI: - mpi_comm = MPI.COMM_WORLD - mpi_rank = mpi_comm.Get_rank() - if mpi_rank == 0: - # self.stime = MPI.Wtime() - # self.stime = time_module.perf_counter() - self.stime = time_module.process_time() - else: - self.stime = time_module.perf_counter() - - def stop_timing(self): - if MPI: - mpi_comm = MPI.COMM_WORLD - mpi_rank = mpi_comm.Get_rank() - if mpi_rank == 0: - # self.etime = MPI.Wtime() - # self.etime = time_module.perf_counter() - self.etime = time_module.process_time() - else: - self.etime = time_module.perf_counter() - - def accumulate_timing(self): - if MPI: - mpi_comm = MPI.COMM_WORLD - mpi_rank = mpi_comm.Get_rank() - if mpi_rank == 0: - self.mtime += (self.etime-self.stime) - else: - self.mtime = 0 - else: - self.mtime += (self.etime-self.stime) - - def advance_iteration(self): - if MPI: - mpi_comm = MPI.COMM_WORLD - mpi_rank = mpi_comm.Get_rank() - if mpi_rank == 0: - self._times_steps.append(self.mtime) - self._samples.append(self._iter) - self._iter += 1 - self.mtime = 0 - else: - self._times_steps.append(self.mtime) - self._samples.append(self._iter) - self._iter += 1 - self.mtime = 0 - - -class ParticleSet_ParamLogging(): - _samples = [] - _params = [] - _iter = 0 - - def __init__(self): - self._samples = [] - self._params = [] - self._iter = 0 - - @property - def samples(self): - return self._samples - - @property - def params(self): - return self._params - - def get_params(self): - return self._params - - def get_param(self, index): - return self._params[index] - - def __len__(self): - return len(self._samples) - - def advance_iteration(self, param): - if MPI: - mpi_comm = MPI.COMM_WORLD - mpi_rank = mpi_comm.Get_rank() - if mpi_rank == 0: - self._params.append(param) - self._samples.append(self._iter) - self._iter += 1 - else: - self._params.append(param) - self._samples.append(self._iter) - self._iter += 1 - - class ParticleSet_Benchmark(ParticleSet): def __init__(self, fieldset, pclass=JITParticle, lon=None, lat=None, depth=None, time=None, repeatdt=None, diff --git a/performance/benchmark_stommel.py b/performance/benchmark_stommel.py index eb0e125b64..6bf63192c4 100644 --- a/performance/benchmark_stommel.py +++ b/performance/benchmark_stommel.py @@ -426,8 +426,11 @@ def Age(particle, fieldset, time): mpi_comm = MPI.COMM_WORLD mpi_rank = mpi_comm.Get_rank() if mpi_rank==0: + # endtime = ostime.time() + # endtime = MPI.Wtime() endtime = ostime.process_time() else: + #endtime = ostime.time() endtime = ostime.process_time() if args.write_out: From 60ee34c87203518da4d82e13ff10ce284af2e9d8 Mon Sep 17 00:00:00 2001 From: Christian Kehl Date: Tue, 7 Jul 2020 20:02:48 +0200 Subject: [PATCH 28/56] list_of_pset_array_trials - initial addenda for list-of-particles version --- parcels/particleset.py | 1 - 1 file changed, 1 deletion(-) diff --git a/parcels/particleset.py b/parcels/particleset.py index f1787c1f21..667c1d2285 100644 --- a/parcels/particleset.py +++ b/parcels/particleset.py @@ -30,7 +30,6 @@ __all__ = ['ParticleSet'] - class ParticleSet(object): """Container class for storing particle and executing kernel over them. From c5ea2744d9d5f26cf3cb50fd67247f476fe30982 Mon Sep 17 00:00:00 2001 From: Christian Kehl Date: Wed, 8 Jul 2020 17:55:10 +0200 Subject: [PATCH 29/56] list_of_pset_array_trials - adapted the particle file to work with the new pset (basically: change the iterator). structure to be wrapped up and tested now. --- parcels/particlefile.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/parcels/particlefile.py b/parcels/particlefile.py index 168285406c..02fa4cc39e 100644 --- a/parcels/particlefile.py +++ b/parcels/particlefile.py @@ -260,6 +260,7 @@ def convert_pset_to_dict(self, pset, time, deleted_only=False): data_dict[var] = np.array([getattr(p, var) for p in pset_towrite]) self.maxid_written = np.max([self.maxid_written, np.max(data_dict['id'])]) + # pset_errs = [p for sublist in pset_towrite.particles for p in sublist if p.state != ErrorCode.Delete and abs(time - p.time) > 1e-3] pset_errs = [p for p in pset_towrite if p.state != ErrorCode.Delete and abs(time-p.time) > 1e-3] # pset_errs = [p for sublist in pset_towrite.particles_a for p in sublist if p.state != ErrorCode.Delete and abs(time - p.time) > 1e-3] for p in pset_errs: @@ -271,6 +272,7 @@ def convert_pset_to_dict(self, pset, time, deleted_only=False): if len(self.var_names_once) > 0: # first_write = [p for p in pset if (p.id not in self.written_once) and _is_particle_started_yet(p, time)] + # first_write = [p for sublist in pset.particles for p in sublist if (p.id not in self.written_once) and _is_particle_started_yet(p, time)] first_write = [p for sublist in pset if sublist is not None From 722d44c19d28c1f5c86c20ce8164c439d68227b6 Mon Sep 17 00:00:00 2001 From: Christian Kehl Date: Tue, 14 Jul 2020 18:20:29 +0200 Subject: [PATCH 30/56] list_of_pset_array_trials - fixing linter errors --- parcels/particleset.py | 1 + parcels/particleset_benchmark.py | 76 ++++++++++++++++++++++++++++++++ 2 files changed, 77 insertions(+) diff --git a/parcels/particleset.py b/parcels/particleset.py index 667c1d2285..f1787c1f21 100644 --- a/parcels/particleset.py +++ b/parcels/particleset.py @@ -30,6 +30,7 @@ __all__ = ['ParticleSet'] + class ParticleSet(object): """Container class for storing particle and executing kernel over them. diff --git a/parcels/particleset_benchmark.py b/parcels/particleset_benchmark.py index ece52e391b..d71943c5ac 100644 --- a/parcels/particleset_benchmark.py +++ b/parcels/particleset_benchmark.py @@ -51,6 +51,82 @@ def measure_mem_usage(): USE_ASYNC_MEMLOG = False USE_RUSE_SYNC_MEMLOG = False # can be faulty +class ParticleSet_TimingLog(): + stime = 0 + etime = 0 + mtime = 0 + samples = [] + times_steps = [] + _iter = 0 + + def start_timing(self): + if MPI: + mpi_comm = MPI.COMM_WORLD + mpi_rank = mpi_comm.Get_rank() + if mpi_rank == 0: + # self.stime = MPI.Wtime() + # self.stime = time_module.perf_counter() + self.stime = time_module.process_time() + else: + self.stime = time_module.perf_counter() + + def stop_timing(self): + if MPI: + mpi_comm = MPI.COMM_WORLD + mpi_rank = mpi_comm.Get_rank() + if mpi_rank == 0: + # self.etime = MPI.Wtime() + # self.etime = time_module.perf_counter() + self.etime = time_module.process_time() + else: + self.etime = time_module.perf_counter() + + def accumulate_timing(self): + if MPI: + mpi_comm = MPI.COMM_WORLD + mpi_rank = mpi_comm.Get_rank() + if mpi_rank == 0: + self.mtime += (self.etime-self.stime) + else: + self.mtime = 0 + else: + self.mtime += (self.etime-self.stime) + + def advance_iteration(self): + if MPI: + mpi_comm = MPI.COMM_WORLD + mpi_rank = mpi_comm.Get_rank() + if mpi_rank == 0: + self.times_steps.append(self.mtime) + self.samples.append(self._iter) + self._iter += 1 + self.mtime = 0 + else: + self.times_steps.append(self.mtime) + self.samples.append(self._iter) + self._iter += 1 + self.mtime = 0 + + +class ParticleSet_ParamLogging(): + samples = [] + params = [] + _iter = 0 + + def advance_iteration(self, param): + if MPI: + mpi_comm = MPI.COMM_WORLD + mpi_rank = mpi_comm.Get_rank() + if mpi_rank == 0: + self.params.append(param) + self.samples.append(self._iter) + self._iter += 1 + else: + self.params.append(param) + self.samples.append(self._iter) + self._iter += 1 + + class ParticleSet_Benchmark(ParticleSet): def __init__(self, fieldset, pclass=JITParticle, lon=None, lat=None, depth=None, time=None, repeatdt=None, From 081f888fb10d273cd4273d5e424ed822a231a575 Mon Sep 17 00:00:00 2001 From: Christian Kehl Date: Thu, 16 Jul 2020 17:07:21 +0200 Subject: [PATCH 31/56] list_of_pset_array_trials - added benchmark refactoring. tested and working (stommel). --- parcels/particleset_benchmark.py | 71 ++++++++++++++++++++++++++------ performance/benchmark_stommel.py | 3 -- 2 files changed, 59 insertions(+), 15 deletions(-) diff --git a/parcels/particleset_benchmark.py b/parcels/particleset_benchmark.py index d71943c5ac..999f02a402 100644 --- a/parcels/particleset_benchmark.py +++ b/parcels/particleset_benchmark.py @@ -55,10 +55,35 @@ class ParticleSet_TimingLog(): stime = 0 etime = 0 mtime = 0 - samples = [] - times_steps = [] + _samples = [] + _times_steps = [] _iter = 0 + def __init__(self): + self.stime = 0 + self.etime = 0 + self.mtime = 0 + self._samples = [] + self._times_steps = [] + self._iter = 0 + + @property + def timing(self): + return self._times_steps + + @property + def samples(self): + return self._samples + + def __len__(self): + return len(self._samples) + + def get_values(self): + return self._times_steps + + def get_value(self, index): + return self._times_steps[index] + def start_timing(self): if MPI: mpi_comm = MPI.COMM_WORLD @@ -97,33 +122,55 @@ def advance_iteration(self): mpi_comm = MPI.COMM_WORLD mpi_rank = mpi_comm.Get_rank() if mpi_rank == 0: - self.times_steps.append(self.mtime) - self.samples.append(self._iter) + self._times_steps.append(self.mtime) + self._samples.append(self._iter) self._iter += 1 self.mtime = 0 else: - self.times_steps.append(self.mtime) - self.samples.append(self._iter) + self._times_steps.append(self.mtime) + self._samples.append(self._iter) self._iter += 1 self.mtime = 0 class ParticleSet_ParamLogging(): - samples = [] - params = [] + _samples = [] + _params = [] _iter = 0 + def __init__(self): + self._samples = [] + self._params = [] + self._iter = 0 + + @property + def samples(self): + return self._samples + + @property + def params(self): + return self._params + + def get_params(self): + return self._params + + def get_param(self, index): + return self._params[index] + + def __len__(self): + return len(self._samples) + def advance_iteration(self, param): if MPI: mpi_comm = MPI.COMM_WORLD mpi_rank = mpi_comm.Get_rank() if mpi_rank == 0: - self.params.append(param) - self.samples.append(self._iter) + self._params.append(param) + self._samples.append(self._iter) self._iter += 1 else: - self.params.append(param) - self.samples.append(self._iter) + self._params.append(param) + self._samples.append(self._iter) self._iter += 1 diff --git a/performance/benchmark_stommel.py b/performance/benchmark_stommel.py index 6bf63192c4..eb0e125b64 100644 --- a/performance/benchmark_stommel.py +++ b/performance/benchmark_stommel.py @@ -426,11 +426,8 @@ def Age(particle, fieldset, time): mpi_comm = MPI.COMM_WORLD mpi_rank = mpi_comm.Get_rank() if mpi_rank==0: - # endtime = ostime.time() - # endtime = MPI.Wtime() endtime = ostime.process_time() else: - #endtime = ostime.time() endtime = ostime.process_time() if args.write_out: From 8601514e0b7566c30f899a0e85bfbf623ad47171 Mon Sep 17 00:00:00 2001 From: Christian Kehl Date: Mon, 20 Jul 2020 23:06:58 +0200 Subject: [PATCH 32/56] list_of_pset_array_trials - rebase to benchmarking branch rev 54e4c18 --- parcels/particleset_benchmark.py | 123 ------------------------------- performance/benchmark_stommel.py | 3 + 2 files changed, 3 insertions(+), 123 deletions(-) diff --git a/parcels/particleset_benchmark.py b/parcels/particleset_benchmark.py index 999f02a402..ece52e391b 100644 --- a/parcels/particleset_benchmark.py +++ b/parcels/particleset_benchmark.py @@ -51,129 +51,6 @@ def measure_mem_usage(): USE_ASYNC_MEMLOG = False USE_RUSE_SYNC_MEMLOG = False # can be faulty -class ParticleSet_TimingLog(): - stime = 0 - etime = 0 - mtime = 0 - _samples = [] - _times_steps = [] - _iter = 0 - - def __init__(self): - self.stime = 0 - self.etime = 0 - self.mtime = 0 - self._samples = [] - self._times_steps = [] - self._iter = 0 - - @property - def timing(self): - return self._times_steps - - @property - def samples(self): - return self._samples - - def __len__(self): - return len(self._samples) - - def get_values(self): - return self._times_steps - - def get_value(self, index): - return self._times_steps[index] - - def start_timing(self): - if MPI: - mpi_comm = MPI.COMM_WORLD - mpi_rank = mpi_comm.Get_rank() - if mpi_rank == 0: - # self.stime = MPI.Wtime() - # self.stime = time_module.perf_counter() - self.stime = time_module.process_time() - else: - self.stime = time_module.perf_counter() - - def stop_timing(self): - if MPI: - mpi_comm = MPI.COMM_WORLD - mpi_rank = mpi_comm.Get_rank() - if mpi_rank == 0: - # self.etime = MPI.Wtime() - # self.etime = time_module.perf_counter() - self.etime = time_module.process_time() - else: - self.etime = time_module.perf_counter() - - def accumulate_timing(self): - if MPI: - mpi_comm = MPI.COMM_WORLD - mpi_rank = mpi_comm.Get_rank() - if mpi_rank == 0: - self.mtime += (self.etime-self.stime) - else: - self.mtime = 0 - else: - self.mtime += (self.etime-self.stime) - - def advance_iteration(self): - if MPI: - mpi_comm = MPI.COMM_WORLD - mpi_rank = mpi_comm.Get_rank() - if mpi_rank == 0: - self._times_steps.append(self.mtime) - self._samples.append(self._iter) - self._iter += 1 - self.mtime = 0 - else: - self._times_steps.append(self.mtime) - self._samples.append(self._iter) - self._iter += 1 - self.mtime = 0 - - -class ParticleSet_ParamLogging(): - _samples = [] - _params = [] - _iter = 0 - - def __init__(self): - self._samples = [] - self._params = [] - self._iter = 0 - - @property - def samples(self): - return self._samples - - @property - def params(self): - return self._params - - def get_params(self): - return self._params - - def get_param(self, index): - return self._params[index] - - def __len__(self): - return len(self._samples) - - def advance_iteration(self, param): - if MPI: - mpi_comm = MPI.COMM_WORLD - mpi_rank = mpi_comm.Get_rank() - if mpi_rank == 0: - self._params.append(param) - self._samples.append(self._iter) - self._iter += 1 - else: - self._params.append(param) - self._samples.append(self._iter) - self._iter += 1 - - class ParticleSet_Benchmark(ParticleSet): def __init__(self, fieldset, pclass=JITParticle, lon=None, lat=None, depth=None, time=None, repeatdt=None, diff --git a/performance/benchmark_stommel.py b/performance/benchmark_stommel.py index eb0e125b64..6bf63192c4 100644 --- a/performance/benchmark_stommel.py +++ b/performance/benchmark_stommel.py @@ -426,8 +426,11 @@ def Age(particle, fieldset, time): mpi_comm = MPI.COMM_WORLD mpi_rank = mpi_comm.Get_rank() if mpi_rank==0: + # endtime = ostime.time() + # endtime = MPI.Wtime() endtime = ostime.process_time() else: + #endtime = ostime.time() endtime = ostime.process_time() if args.write_out: From d98d9fc0bbfe07de207efaa27232b67900ef6e45 Mon Sep 17 00:00:00 2001 From: Christian Kehl Date: Tue, 7 Jul 2020 20:02:48 +0200 Subject: [PATCH 33/56] list_of_pset_array_trials - initial addenda for list-of-particles version --- parcels/particleset.py | 1 - 1 file changed, 1 deletion(-) diff --git a/parcels/particleset.py b/parcels/particleset.py index f1787c1f21..667c1d2285 100644 --- a/parcels/particleset.py +++ b/parcels/particleset.py @@ -30,7 +30,6 @@ __all__ = ['ParticleSet'] - class ParticleSet(object): """Container class for storing particle and executing kernel over them. From 791f0670dea3b58474e79be1f7776cf6d7d3bfb9 Mon Sep 17 00:00:00 2001 From: Christian Kehl Date: Wed, 8 Jul 2020 16:10:52 +0200 Subject: [PATCH 34/56] list_of_pset_array_trials - wrapped up list-of-array adaptations in the ParticleSet class. --- parcels/particleset.py | 1 + 1 file changed, 1 insertion(+) diff --git a/parcels/particleset.py b/parcels/particleset.py index 667c1d2285..59e4b45153 100644 --- a/parcels/particleset.py +++ b/parcels/particleset.py @@ -276,6 +276,7 @@ def get_list_array_index(self, pdata): return (bracket_index, slot_index) elif isinstance(pdata, self.pclass): # == parameter is a particle itself == # + # nparticle = 0 bracket_index = 0 bracket_info = self.pid_mapping_bounds[bracket_index] while bracket_index < len(self._plist) and ((pdata.id < bracket_info[0]) or (pdata.id > bracket_info[1])): From e71e4660089a52f582d0543ccdce4ca4ca77f5f1 Mon Sep 17 00:00:00 2001 From: Christian Kehl Date: Tue, 14 Jul 2020 18:20:29 +0200 Subject: [PATCH 35/56] list_of_pset_array_trials - fixing linter errors --- parcels/particleset.py | 2 +- parcels/particleset_benchmark.py | 76 ++++++++++++++++++++++++++++++++ 2 files changed, 77 insertions(+), 1 deletion(-) diff --git a/parcels/particleset.py b/parcels/particleset.py index 59e4b45153..f1787c1f21 100644 --- a/parcels/particleset.py +++ b/parcels/particleset.py @@ -30,6 +30,7 @@ __all__ = ['ParticleSet'] + class ParticleSet(object): """Container class for storing particle and executing kernel over them. @@ -276,7 +277,6 @@ def get_list_array_index(self, pdata): return (bracket_index, slot_index) elif isinstance(pdata, self.pclass): # == parameter is a particle itself == # - # nparticle = 0 bracket_index = 0 bracket_info = self.pid_mapping_bounds[bracket_index] while bracket_index < len(self._plist) and ((pdata.id < bracket_info[0]) or (pdata.id > bracket_info[1])): diff --git a/parcels/particleset_benchmark.py b/parcels/particleset_benchmark.py index ece52e391b..d71943c5ac 100644 --- a/parcels/particleset_benchmark.py +++ b/parcels/particleset_benchmark.py @@ -51,6 +51,82 @@ def measure_mem_usage(): USE_ASYNC_MEMLOG = False USE_RUSE_SYNC_MEMLOG = False # can be faulty +class ParticleSet_TimingLog(): + stime = 0 + etime = 0 + mtime = 0 + samples = [] + times_steps = [] + _iter = 0 + + def start_timing(self): + if MPI: + mpi_comm = MPI.COMM_WORLD + mpi_rank = mpi_comm.Get_rank() + if mpi_rank == 0: + # self.stime = MPI.Wtime() + # self.stime = time_module.perf_counter() + self.stime = time_module.process_time() + else: + self.stime = time_module.perf_counter() + + def stop_timing(self): + if MPI: + mpi_comm = MPI.COMM_WORLD + mpi_rank = mpi_comm.Get_rank() + if mpi_rank == 0: + # self.etime = MPI.Wtime() + # self.etime = time_module.perf_counter() + self.etime = time_module.process_time() + else: + self.etime = time_module.perf_counter() + + def accumulate_timing(self): + if MPI: + mpi_comm = MPI.COMM_WORLD + mpi_rank = mpi_comm.Get_rank() + if mpi_rank == 0: + self.mtime += (self.etime-self.stime) + else: + self.mtime = 0 + else: + self.mtime += (self.etime-self.stime) + + def advance_iteration(self): + if MPI: + mpi_comm = MPI.COMM_WORLD + mpi_rank = mpi_comm.Get_rank() + if mpi_rank == 0: + self.times_steps.append(self.mtime) + self.samples.append(self._iter) + self._iter += 1 + self.mtime = 0 + else: + self.times_steps.append(self.mtime) + self.samples.append(self._iter) + self._iter += 1 + self.mtime = 0 + + +class ParticleSet_ParamLogging(): + samples = [] + params = [] + _iter = 0 + + def advance_iteration(self, param): + if MPI: + mpi_comm = MPI.COMM_WORLD + mpi_rank = mpi_comm.Get_rank() + if mpi_rank == 0: + self.params.append(param) + self.samples.append(self._iter) + self._iter += 1 + else: + self.params.append(param) + self.samples.append(self._iter) + self._iter += 1 + + class ParticleSet_Benchmark(ParticleSet): def __init__(self, fieldset, pclass=JITParticle, lon=None, lat=None, depth=None, time=None, repeatdt=None, From 14827a393503720f269d727278a36e88da11241c Mon Sep 17 00:00:00 2001 From: Christian Kehl Date: Thu, 16 Jul 2020 17:07:21 +0200 Subject: [PATCH 36/56] list_of_pset_array_trials - added benchmark refactoring. tested and working (stommel). --- parcels/particleset_benchmark.py | 71 ++++++++++++++++++++++++++------ performance/benchmark_stommel.py | 3 -- 2 files changed, 59 insertions(+), 15 deletions(-) diff --git a/parcels/particleset_benchmark.py b/parcels/particleset_benchmark.py index d71943c5ac..999f02a402 100644 --- a/parcels/particleset_benchmark.py +++ b/parcels/particleset_benchmark.py @@ -55,10 +55,35 @@ class ParticleSet_TimingLog(): stime = 0 etime = 0 mtime = 0 - samples = [] - times_steps = [] + _samples = [] + _times_steps = [] _iter = 0 + def __init__(self): + self.stime = 0 + self.etime = 0 + self.mtime = 0 + self._samples = [] + self._times_steps = [] + self._iter = 0 + + @property + def timing(self): + return self._times_steps + + @property + def samples(self): + return self._samples + + def __len__(self): + return len(self._samples) + + def get_values(self): + return self._times_steps + + def get_value(self, index): + return self._times_steps[index] + def start_timing(self): if MPI: mpi_comm = MPI.COMM_WORLD @@ -97,33 +122,55 @@ def advance_iteration(self): mpi_comm = MPI.COMM_WORLD mpi_rank = mpi_comm.Get_rank() if mpi_rank == 0: - self.times_steps.append(self.mtime) - self.samples.append(self._iter) + self._times_steps.append(self.mtime) + self._samples.append(self._iter) self._iter += 1 self.mtime = 0 else: - self.times_steps.append(self.mtime) - self.samples.append(self._iter) + self._times_steps.append(self.mtime) + self._samples.append(self._iter) self._iter += 1 self.mtime = 0 class ParticleSet_ParamLogging(): - samples = [] - params = [] + _samples = [] + _params = [] _iter = 0 + def __init__(self): + self._samples = [] + self._params = [] + self._iter = 0 + + @property + def samples(self): + return self._samples + + @property + def params(self): + return self._params + + def get_params(self): + return self._params + + def get_param(self, index): + return self._params[index] + + def __len__(self): + return len(self._samples) + def advance_iteration(self, param): if MPI: mpi_comm = MPI.COMM_WORLD mpi_rank = mpi_comm.Get_rank() if mpi_rank == 0: - self.params.append(param) - self.samples.append(self._iter) + self._params.append(param) + self._samples.append(self._iter) self._iter += 1 else: - self.params.append(param) - self.samples.append(self._iter) + self._params.append(param) + self._samples.append(self._iter) self._iter += 1 diff --git a/performance/benchmark_stommel.py b/performance/benchmark_stommel.py index 6bf63192c4..eb0e125b64 100644 --- a/performance/benchmark_stommel.py +++ b/performance/benchmark_stommel.py @@ -426,11 +426,8 @@ def Age(particle, fieldset, time): mpi_comm = MPI.COMM_WORLD mpi_rank = mpi_comm.Get_rank() if mpi_rank==0: - # endtime = ostime.time() - # endtime = MPI.Wtime() endtime = ostime.process_time() else: - #endtime = ostime.time() endtime = ostime.process_time() if args.write_out: From afacbd3d68ceff4e76973a8abc4f76f6da07d8d8 Mon Sep 17 00:00:00 2001 From: Christian Kehl Date: Mon, 20 Jul 2020 23:06:58 +0200 Subject: [PATCH 37/56] list_of_pset_array_trials - rebase to benchmarking branch rev 54e4c18 --- parcels/particleset_benchmark.py | 123 ------------------------------- performance/benchmark_perlin.py | 9 --- performance/benchmark_stommel.py | 3 + 3 files changed, 3 insertions(+), 132 deletions(-) diff --git a/parcels/particleset_benchmark.py b/parcels/particleset_benchmark.py index 999f02a402..ece52e391b 100644 --- a/parcels/particleset_benchmark.py +++ b/parcels/particleset_benchmark.py @@ -51,129 +51,6 @@ def measure_mem_usage(): USE_ASYNC_MEMLOG = False USE_RUSE_SYNC_MEMLOG = False # can be faulty -class ParticleSet_TimingLog(): - stime = 0 - etime = 0 - mtime = 0 - _samples = [] - _times_steps = [] - _iter = 0 - - def __init__(self): - self.stime = 0 - self.etime = 0 - self.mtime = 0 - self._samples = [] - self._times_steps = [] - self._iter = 0 - - @property - def timing(self): - return self._times_steps - - @property - def samples(self): - return self._samples - - def __len__(self): - return len(self._samples) - - def get_values(self): - return self._times_steps - - def get_value(self, index): - return self._times_steps[index] - - def start_timing(self): - if MPI: - mpi_comm = MPI.COMM_WORLD - mpi_rank = mpi_comm.Get_rank() - if mpi_rank == 0: - # self.stime = MPI.Wtime() - # self.stime = time_module.perf_counter() - self.stime = time_module.process_time() - else: - self.stime = time_module.perf_counter() - - def stop_timing(self): - if MPI: - mpi_comm = MPI.COMM_WORLD - mpi_rank = mpi_comm.Get_rank() - if mpi_rank == 0: - # self.etime = MPI.Wtime() - # self.etime = time_module.perf_counter() - self.etime = time_module.process_time() - else: - self.etime = time_module.perf_counter() - - def accumulate_timing(self): - if MPI: - mpi_comm = MPI.COMM_WORLD - mpi_rank = mpi_comm.Get_rank() - if mpi_rank == 0: - self.mtime += (self.etime-self.stime) - else: - self.mtime = 0 - else: - self.mtime += (self.etime-self.stime) - - def advance_iteration(self): - if MPI: - mpi_comm = MPI.COMM_WORLD - mpi_rank = mpi_comm.Get_rank() - if mpi_rank == 0: - self._times_steps.append(self.mtime) - self._samples.append(self._iter) - self._iter += 1 - self.mtime = 0 - else: - self._times_steps.append(self.mtime) - self._samples.append(self._iter) - self._iter += 1 - self.mtime = 0 - - -class ParticleSet_ParamLogging(): - _samples = [] - _params = [] - _iter = 0 - - def __init__(self): - self._samples = [] - self._params = [] - self._iter = 0 - - @property - def samples(self): - return self._samples - - @property - def params(self): - return self._params - - def get_params(self): - return self._params - - def get_param(self, index): - return self._params[index] - - def __len__(self): - return len(self._samples) - - def advance_iteration(self, param): - if MPI: - mpi_comm = MPI.COMM_WORLD - mpi_rank = mpi_comm.Get_rank() - if mpi_rank == 0: - self._params.append(param) - self._samples.append(self._iter) - self._iter += 1 - else: - self._params.append(param) - self._samples.append(self._iter) - self._iter += 1 - - class ParticleSet_Benchmark(ParticleSet): def __init__(self, fieldset, pclass=JITParticle, lon=None, lat=None, depth=None, time=None, repeatdt=None, diff --git a/performance/benchmark_perlin.py b/performance/benchmark_perlin.py index e0f62c2062..e2a8246fc9 100644 --- a/performance/benchmark_perlin.py +++ b/performance/benchmark_perlin.py @@ -100,15 +100,6 @@ def perlin_fieldset_from_numpy(periodic_wrap=False, write_out=False): U = np.transpose(U, (0,2,1)) V = perlin2d.generate_fractal_noise_temporal2d(img_shape, tsteps, (perlinres[1], perlinres[2]), noctaves, perlin_persistence, max_shift=((-1, 2), (-1, 2))) V = np.transpose(V, (0,2,1)) - # V = np.swapaxes(V, 1, 2) - # print("V-statistics - min: {:10.7f}; max: {:10.7f}; avg. {:10.7f}; std_dev: {:10.7f}".format(V.min(initial=0), V.max(initial=0), V.mean(), V.std())) - - # U = perlin3d.generate_fractal_noise_3d(img_shape, perlinres, noctaves, perlin_persistence) * scalefac - # U = np.transpose(U, (0,2,1)) - # sys.stdout.write("U field shape: {} - [tdim][ydim][xdim]=[{}][{}][{}]\n".format(U.shape, time.shape[0], lat.shape[0], lon.shape[0])) - # V = perlin3d.generate_fractal_noise_3d(img_shape, perlinres, noctaves, perlin_persistence) * scalefac - # V = np.transpose(V, (0,2,1)) - # sys.stdout.write("V field shape: {} - [tdim][ydim][xdim]=[{}][{}][{}]\n".format(V.shape, time.shape[0], lat.shape[0], lon.shape[0])) U *= scalefac V *= scalefac diff --git a/performance/benchmark_stommel.py b/performance/benchmark_stommel.py index eb0e125b64..6bf63192c4 100644 --- a/performance/benchmark_stommel.py +++ b/performance/benchmark_stommel.py @@ -426,8 +426,11 @@ def Age(particle, fieldset, time): mpi_comm = MPI.COMM_WORLD mpi_rank = mpi_comm.Get_rank() if mpi_rank==0: + # endtime = ostime.time() + # endtime = MPI.Wtime() endtime = ostime.process_time() else: + #endtime = ostime.time() endtime = ostime.process_time() if args.write_out: From e7dd4f652e39f6acf997b68e8aa11bf59e98919f Mon Sep 17 00:00:00 2001 From: Christian Kehl Date: Tue, 7 Jul 2020 20:02:48 +0200 Subject: [PATCH 38/56] list_of_pset_array_trials - initial addenda for list-of-particles version --- parcels/particle.py | 2 +- parcels/particleset.py | 40 +--------------------------------------- 2 files changed, 2 insertions(+), 40 deletions(-) diff --git a/parcels/particle.py b/parcels/particle.py index 0762b46749..1800b52b77 100644 --- a/parcels/particle.py +++ b/parcels/particle.py @@ -260,5 +260,5 @@ def __init__(self, *args, **kwargs): setattr(self, 'c'+index, getattr(self, index+'p').value) def update_cptr(self, new_ptr): - # np.copyto(new_ptr, self._cptr) + np.copyto(new_ptr, self._cptr) self._cptr = new_ptr diff --git a/parcels/particleset.py b/parcels/particleset.py index f1787c1f21..45553ac86c 100644 --- a/parcels/particleset.py +++ b/parcels/particleset.py @@ -30,7 +30,6 @@ __all__ = ['ParticleSet'] - class ParticleSet(object): """Container class for storing particle and executing kernel over them. @@ -154,9 +153,7 @@ def convert_to_array(var): 'lon lat depth precision should be set to either np.float32 or np.float64' JITParticle.set_lonlatdepth_dtype(self.lonlatdepth_dtype) - # == this should be a list not a dict! == # - self._pid_mapping_bounds = [] # plist bracket_index -> (min_id, max_id, size_bracket) - # ======================================= # + self._pid_mapping_bounds = [] # plist bracket_index -> (min_id, max_id, size_bracket) self.nlist_limit = 4096 # self.particles = np.empty(lon.size, dtype=pclass) assert lon.shape[0] == lat.shape[0], ('Length of lon and lat do not match.') @@ -501,34 +498,27 @@ def size(self): for bracket_index in range(len(self._pid_mapping_bounds)): nparticles += self._pid_mapping_bounds[bracket_index][2] return nparticles - # return self.particles.size def __repr__(self): return "\n".join([str(p) for sublist in self._plist for p in sublist]) def __len__(self): - # return len(self._plist) return self.size def __getitem__(self, key): bracket_index = 0 start_index = 0 end_index = self._pid_mapping_bounds[bracket_index][2] - # while start_index= key: while key >= end_index: start_index = end_index end_index = start_index+self._pid_mapping_bounds[bracket_index][2] bracket_index += 1 slot_index = key - start_index return self._plist[bracket_index][slot_index] - # return self.particles[key] def __setitem__(self, key, value): assert isinstance(key, tuple) and len(key) == 2, ("Error: trying to set/address data item without (bracket_index, slot_index) key.") self._plist[key[0]][key[1]] = value - # def __setitem__(self, key, value): - # self.particles[key] = value def __iadd__(self, particles): self.add(particles) @@ -557,18 +547,6 @@ def add(self, particles): else: raise NotImplementedError('Only ParticleSets can be added to a ParticleSet') - # if isinstance(particles, ParticleSet): - # particles = particles.particles - # else: - # raise NotImplementedError('Only ParticleSets can be added to a ParticleSet') - # self.particles = np.append(self.particles, particles) - # if self.ptype.uses_jit: - # particles_data = [p._cptr for p in particles] - # self._particle_data = np.append(self._particle_data, particles_data) - # # Update C-pointer on particles - # for p, pdata in zip(self.particles, self._particle_data): - # p._cptr = pdata - def _merge_brackets_(self): lw_bound_nlist = int(self.nlist_limit / 2) nmerges = 1 @@ -685,26 +663,10 @@ def remove(self, indices): self._remove_(bracket_index, local_indices) self._merge_brackets_() - # if isinstance(indices, collections.Iterable): - # particles = [self.particles[i] for i in indices] - # else: - # particles = self.particles[indices] - # self.particles = np.delete(self.particles, indices) - # if self.ptype.uses_jit: - # self._particle_data = np.delete(self._particle_data, indices) - # # == Update C-pointer on particles == # - # for p, pdata in zip(self.particles, self._particle_data): - # p._cptr = pdata - # return particles - def remove_local_particles(self, bracket_index, local_indices): return self._remove_(bracket_index, local_indices) def _remove_(self, bracket_index, local_indices): - # if isinstance(indices, collections.Iterable): - # particles = [self.particles[i] for i in indices] - # else: - # particles = self.particles[indices] if isinstance(local_indices, collections.Iterable): local_indices = np.array(local_indices) return_p = self._plist[bracket_index][local_indices] From 2080637ed53b53b7307635d393671328dbea44cd Mon Sep 17 00:00:00 2001 From: Christian Kehl Date: Wed, 8 Jul 2020 16:10:52 +0200 Subject: [PATCH 39/56] list_of_pset_array_trials - wrapped up list-of-array adaptations in the ParticleSet class. --- parcels/particleset.py | 10 ---------- 1 file changed, 10 deletions(-) diff --git a/parcels/particleset.py b/parcels/particleset.py index 45553ac86c..1cd035a209 100644 --- a/parcels/particleset.py +++ b/parcels/particleset.py @@ -932,16 +932,6 @@ def density(self, field=None, particle_val=None, relative=False, area_scale=Fals _, _, _, xi, yi, _ = field.search_indices(p.lon, p.lat, p.depth, 0, 0, search2D=True) density[yi, xi] += particle_val[bracket_index][slot_index] - # for pi, p in enumerate(self.particles): - # try: # breaks if either p.xi, p.yi, p.zi, p.ti do not exist (in scipy) or field not in fieldset - # if p.ti[field.igrid] < 0: # xi, yi, zi, ti, not initialised - # raise('error') - # xi = p.xi[field.igrid] - # yi = p.yi[field.igrid] - # except: - # _, _, _, xi, yi, _ = field.search_indices(p.lon, p.lat, p.depth, 0, 0, search2D=True) - # density[yi, xi] += particle_val[pi] - if relative: psum = 0 for bracket_index, subfield in enumerate(particle_val): From c7005dfcd8c64b776c6f0d537066f44cfcab1f19 Mon Sep 17 00:00:00 2001 From: Christian Kehl Date: Wed, 8 Jul 2020 16:48:14 +0200 Subject: [PATCH 40/56] list_of_pset_array_trials - adapted the kernel to work with list-of-array particles. TODO: particle file. --- parcels/kernel.py | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/parcels/kernel.py b/parcels/kernel.py index e98e122540..3193532168 100644 --- a/parcels/kernel.py +++ b/parcels/kernel.py @@ -415,8 +415,7 @@ def execute(self, pset, endtime, dt, recovery=None, output_file=None, execute_on # Remove all particles that signalled deletion self.remove_deleted(pset, output_file=output_file, endtime=endtime) - # == Identify particles that threw errors == # - # error_particles = [p for p in pset.particles if p.state not in [ErrorCode.Success, ErrorCode.Evaluate]] + # == Identify particles that threw errors == #\ error_particles = [p for subfield in pset.particles_a for p in subfield if p.state not in [ErrorCode.Success, ErrorCode.Evaluate]] while len(error_particles) > 0: @@ -448,7 +447,6 @@ def execute(self, pset, endtime, dt, recovery=None, output_file=None, execute_on else: self.execute_python(pset, endtime, dt) - # error_particles = [p for p in pset.particles if p.state not in [ErrorCode.Success, ErrorCode.Evaluate]] error_particles = [p for subfield in pset.particles_a for p in subfield if p.state not in [ErrorCode.Success, ErrorCode.Evaluate]] def merge(self, kernel): From fee95d1c2631b3e656bbca9b0e64a95237c73bac Mon Sep 17 00:00:00 2001 From: Christian Kehl Date: Wed, 8 Jul 2020 17:55:10 +0200 Subject: [PATCH 41/56] list_of_pset_array_trials - adapted the particle file to work with the new pset (basically: change the iterator). structure to be wrapped up and tested now. --- parcels/particlefile.py | 1 - 1 file changed, 1 deletion(-) diff --git a/parcels/particlefile.py b/parcels/particlefile.py index 02fa4cc39e..df136bed92 100644 --- a/parcels/particlefile.py +++ b/parcels/particlefile.py @@ -232,7 +232,6 @@ def convert_pset_to_dict(self, pset, time, deleted_only=False): if not isinstance(pset, list): pset = pset.particles_a - # if pset.size == 0: if len(pset) == 0: logger.warning("ParticleSet is empty on writing as array at time %g" % time) else: From 4d61245f5375cd9e51f4bc86adc8ff8d0f40a29e Mon Sep 17 00:00:00 2001 From: Christian Kehl Date: Tue, 14 Jul 2020 17:38:01 +0200 Subject: [PATCH 42/56] list_of_pset_array_trials - fixed errors emerging from tests --- parcels/kernel.py | 2 +- parcels/particle.py | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/parcels/kernel.py b/parcels/kernel.py index 3193532168..a9b8dec583 100644 --- a/parcels/kernel.py +++ b/parcels/kernel.py @@ -415,7 +415,7 @@ def execute(self, pset, endtime, dt, recovery=None, output_file=None, execute_on # Remove all particles that signalled deletion self.remove_deleted(pset, output_file=output_file, endtime=endtime) - # == Identify particles that threw errors == #\ + # == Identify particles that threw errors == # error_particles = [p for subfield in pset.particles_a for p in subfield if p.state not in [ErrorCode.Success, ErrorCode.Evaluate]] while len(error_particles) > 0: diff --git a/parcels/particle.py b/parcels/particle.py index 1800b52b77..0762b46749 100644 --- a/parcels/particle.py +++ b/parcels/particle.py @@ -260,5 +260,5 @@ def __init__(self, *args, **kwargs): setattr(self, 'c'+index, getattr(self, index+'p').value) def update_cptr(self, new_ptr): - np.copyto(new_ptr, self._cptr) + # np.copyto(new_ptr, self._cptr) self._cptr = new_ptr From 0fe9417628bdeffb4d1ece5e6988b7f77ef18f5e Mon Sep 17 00:00:00 2001 From: Christian Kehl Date: Tue, 14 Jul 2020 18:20:29 +0200 Subject: [PATCH 43/56] list_of_pset_array_trials - fixing linter errors --- parcels/particleset.py | 3 +- parcels/particleset_benchmark.py | 76 ++++++++++++++++++++++++++++++++ 2 files changed, 78 insertions(+), 1 deletion(-) diff --git a/parcels/particleset.py b/parcels/particleset.py index 1cd035a209..c07e7431eb 100644 --- a/parcels/particleset.py +++ b/parcels/particleset.py @@ -30,6 +30,7 @@ __all__ = ['ParticleSet'] + class ParticleSet(object): """Container class for storing particle and executing kernel over them. @@ -153,7 +154,7 @@ def convert_to_array(var): 'lon lat depth precision should be set to either np.float32 or np.float64' JITParticle.set_lonlatdepth_dtype(self.lonlatdepth_dtype) - self._pid_mapping_bounds = [] # plist bracket_index -> (min_id, max_id, size_bracket) + self._pid_mapping_bounds = [] # plist bracket_index -> (min_id, max_id, size_bracket) self.nlist_limit = 4096 # self.particles = np.empty(lon.size, dtype=pclass) assert lon.shape[0] == lat.shape[0], ('Length of lon and lat do not match.') diff --git a/parcels/particleset_benchmark.py b/parcels/particleset_benchmark.py index ece52e391b..d71943c5ac 100644 --- a/parcels/particleset_benchmark.py +++ b/parcels/particleset_benchmark.py @@ -51,6 +51,82 @@ def measure_mem_usage(): USE_ASYNC_MEMLOG = False USE_RUSE_SYNC_MEMLOG = False # can be faulty +class ParticleSet_TimingLog(): + stime = 0 + etime = 0 + mtime = 0 + samples = [] + times_steps = [] + _iter = 0 + + def start_timing(self): + if MPI: + mpi_comm = MPI.COMM_WORLD + mpi_rank = mpi_comm.Get_rank() + if mpi_rank == 0: + # self.stime = MPI.Wtime() + # self.stime = time_module.perf_counter() + self.stime = time_module.process_time() + else: + self.stime = time_module.perf_counter() + + def stop_timing(self): + if MPI: + mpi_comm = MPI.COMM_WORLD + mpi_rank = mpi_comm.Get_rank() + if mpi_rank == 0: + # self.etime = MPI.Wtime() + # self.etime = time_module.perf_counter() + self.etime = time_module.process_time() + else: + self.etime = time_module.perf_counter() + + def accumulate_timing(self): + if MPI: + mpi_comm = MPI.COMM_WORLD + mpi_rank = mpi_comm.Get_rank() + if mpi_rank == 0: + self.mtime += (self.etime-self.stime) + else: + self.mtime = 0 + else: + self.mtime += (self.etime-self.stime) + + def advance_iteration(self): + if MPI: + mpi_comm = MPI.COMM_WORLD + mpi_rank = mpi_comm.Get_rank() + if mpi_rank == 0: + self.times_steps.append(self.mtime) + self.samples.append(self._iter) + self._iter += 1 + self.mtime = 0 + else: + self.times_steps.append(self.mtime) + self.samples.append(self._iter) + self._iter += 1 + self.mtime = 0 + + +class ParticleSet_ParamLogging(): + samples = [] + params = [] + _iter = 0 + + def advance_iteration(self, param): + if MPI: + mpi_comm = MPI.COMM_WORLD + mpi_rank = mpi_comm.Get_rank() + if mpi_rank == 0: + self.params.append(param) + self.samples.append(self._iter) + self._iter += 1 + else: + self.params.append(param) + self.samples.append(self._iter) + self._iter += 1 + + class ParticleSet_Benchmark(ParticleSet): def __init__(self, fieldset, pclass=JITParticle, lon=None, lat=None, depth=None, time=None, repeatdt=None, From 16a8903299533bd62e716ba99da4225285e67e39 Mon Sep 17 00:00:00 2001 From: Christian Kehl Date: Thu, 16 Jul 2020 17:07:21 +0200 Subject: [PATCH 44/56] list_of_pset_array_trials - added benchmark refactoring. tested and working (stommel). --- parcels/particleset_benchmark.py | 71 ++++++++++++++++++++++++++------ performance/benchmark_stommel.py | 3 -- 2 files changed, 59 insertions(+), 15 deletions(-) diff --git a/parcels/particleset_benchmark.py b/parcels/particleset_benchmark.py index d71943c5ac..999f02a402 100644 --- a/parcels/particleset_benchmark.py +++ b/parcels/particleset_benchmark.py @@ -55,10 +55,35 @@ class ParticleSet_TimingLog(): stime = 0 etime = 0 mtime = 0 - samples = [] - times_steps = [] + _samples = [] + _times_steps = [] _iter = 0 + def __init__(self): + self.stime = 0 + self.etime = 0 + self.mtime = 0 + self._samples = [] + self._times_steps = [] + self._iter = 0 + + @property + def timing(self): + return self._times_steps + + @property + def samples(self): + return self._samples + + def __len__(self): + return len(self._samples) + + def get_values(self): + return self._times_steps + + def get_value(self, index): + return self._times_steps[index] + def start_timing(self): if MPI: mpi_comm = MPI.COMM_WORLD @@ -97,33 +122,55 @@ def advance_iteration(self): mpi_comm = MPI.COMM_WORLD mpi_rank = mpi_comm.Get_rank() if mpi_rank == 0: - self.times_steps.append(self.mtime) - self.samples.append(self._iter) + self._times_steps.append(self.mtime) + self._samples.append(self._iter) self._iter += 1 self.mtime = 0 else: - self.times_steps.append(self.mtime) - self.samples.append(self._iter) + self._times_steps.append(self.mtime) + self._samples.append(self._iter) self._iter += 1 self.mtime = 0 class ParticleSet_ParamLogging(): - samples = [] - params = [] + _samples = [] + _params = [] _iter = 0 + def __init__(self): + self._samples = [] + self._params = [] + self._iter = 0 + + @property + def samples(self): + return self._samples + + @property + def params(self): + return self._params + + def get_params(self): + return self._params + + def get_param(self, index): + return self._params[index] + + def __len__(self): + return len(self._samples) + def advance_iteration(self, param): if MPI: mpi_comm = MPI.COMM_WORLD mpi_rank = mpi_comm.Get_rank() if mpi_rank == 0: - self.params.append(param) - self.samples.append(self._iter) + self._params.append(param) + self._samples.append(self._iter) self._iter += 1 else: - self.params.append(param) - self.samples.append(self._iter) + self._params.append(param) + self._samples.append(self._iter) self._iter += 1 diff --git a/performance/benchmark_stommel.py b/performance/benchmark_stommel.py index 6bf63192c4..eb0e125b64 100644 --- a/performance/benchmark_stommel.py +++ b/performance/benchmark_stommel.py @@ -426,11 +426,8 @@ def Age(particle, fieldset, time): mpi_comm = MPI.COMM_WORLD mpi_rank = mpi_comm.Get_rank() if mpi_rank==0: - # endtime = ostime.time() - # endtime = MPI.Wtime() endtime = ostime.process_time() else: - #endtime = ostime.time() endtime = ostime.process_time() if args.write_out: From e25e8708da20cfdec9ef9ea7a8f2e5f73c2ef778 Mon Sep 17 00:00:00 2001 From: Christian Kehl Date: Mon, 20 Jul 2020 23:06:58 +0200 Subject: [PATCH 45/56] list_of_pset_array_trials - rebase to benchmarking branch rev 54e4c18 --- parcels/particleset_benchmark.py | 123 ------------------------------- performance/benchmark_stommel.py | 3 + 2 files changed, 3 insertions(+), 123 deletions(-) diff --git a/parcels/particleset_benchmark.py b/parcels/particleset_benchmark.py index 999f02a402..ece52e391b 100644 --- a/parcels/particleset_benchmark.py +++ b/parcels/particleset_benchmark.py @@ -51,129 +51,6 @@ def measure_mem_usage(): USE_ASYNC_MEMLOG = False USE_RUSE_SYNC_MEMLOG = False # can be faulty -class ParticleSet_TimingLog(): - stime = 0 - etime = 0 - mtime = 0 - _samples = [] - _times_steps = [] - _iter = 0 - - def __init__(self): - self.stime = 0 - self.etime = 0 - self.mtime = 0 - self._samples = [] - self._times_steps = [] - self._iter = 0 - - @property - def timing(self): - return self._times_steps - - @property - def samples(self): - return self._samples - - def __len__(self): - return len(self._samples) - - def get_values(self): - return self._times_steps - - def get_value(self, index): - return self._times_steps[index] - - def start_timing(self): - if MPI: - mpi_comm = MPI.COMM_WORLD - mpi_rank = mpi_comm.Get_rank() - if mpi_rank == 0: - # self.stime = MPI.Wtime() - # self.stime = time_module.perf_counter() - self.stime = time_module.process_time() - else: - self.stime = time_module.perf_counter() - - def stop_timing(self): - if MPI: - mpi_comm = MPI.COMM_WORLD - mpi_rank = mpi_comm.Get_rank() - if mpi_rank == 0: - # self.etime = MPI.Wtime() - # self.etime = time_module.perf_counter() - self.etime = time_module.process_time() - else: - self.etime = time_module.perf_counter() - - def accumulate_timing(self): - if MPI: - mpi_comm = MPI.COMM_WORLD - mpi_rank = mpi_comm.Get_rank() - if mpi_rank == 0: - self.mtime += (self.etime-self.stime) - else: - self.mtime = 0 - else: - self.mtime += (self.etime-self.stime) - - def advance_iteration(self): - if MPI: - mpi_comm = MPI.COMM_WORLD - mpi_rank = mpi_comm.Get_rank() - if mpi_rank == 0: - self._times_steps.append(self.mtime) - self._samples.append(self._iter) - self._iter += 1 - self.mtime = 0 - else: - self._times_steps.append(self.mtime) - self._samples.append(self._iter) - self._iter += 1 - self.mtime = 0 - - -class ParticleSet_ParamLogging(): - _samples = [] - _params = [] - _iter = 0 - - def __init__(self): - self._samples = [] - self._params = [] - self._iter = 0 - - @property - def samples(self): - return self._samples - - @property - def params(self): - return self._params - - def get_params(self): - return self._params - - def get_param(self, index): - return self._params[index] - - def __len__(self): - return len(self._samples) - - def advance_iteration(self, param): - if MPI: - mpi_comm = MPI.COMM_WORLD - mpi_rank = mpi_comm.Get_rank() - if mpi_rank == 0: - self._params.append(param) - self._samples.append(self._iter) - self._iter += 1 - else: - self._params.append(param) - self._samples.append(self._iter) - self._iter += 1 - - class ParticleSet_Benchmark(ParticleSet): def __init__(self, fieldset, pclass=JITParticle, lon=None, lat=None, depth=None, time=None, repeatdt=None, diff --git a/performance/benchmark_stommel.py b/performance/benchmark_stommel.py index eb0e125b64..6bf63192c4 100644 --- a/performance/benchmark_stommel.py +++ b/performance/benchmark_stommel.py @@ -426,8 +426,11 @@ def Age(particle, fieldset, time): mpi_comm = MPI.COMM_WORLD mpi_rank = mpi_comm.Get_rank() if mpi_rank==0: + # endtime = ostime.time() + # endtime = MPI.Wtime() endtime = ostime.process_time() else: + #endtime = ostime.time() endtime = ostime.process_time() if args.write_out: From 135986a8c390bc008d2a277ed479d8e42ccb7c95 Mon Sep 17 00:00:00 2001 From: Christian Kehl Date: Tue, 7 Jul 2020 20:02:48 +0200 Subject: [PATCH 46/56] list_of_pset_array_trials - initial addenda for list-of-particles version --- parcels/particle.py | 2 +- parcels/particleset.py | 18 ++++++------------ 2 files changed, 7 insertions(+), 13 deletions(-) diff --git a/parcels/particle.py b/parcels/particle.py index 0762b46749..1800b52b77 100644 --- a/parcels/particle.py +++ b/parcels/particle.py @@ -260,5 +260,5 @@ def __init__(self, *args, **kwargs): setattr(self, 'c'+index, getattr(self, index+'p').value) def update_cptr(self, new_ptr): - # np.copyto(new_ptr, self._cptr) + np.copyto(new_ptr, self._cptr) self._cptr = new_ptr diff --git a/parcels/particleset.py b/parcels/particleset.py index c07e7431eb..d2443f1db9 100644 --- a/parcels/particleset.py +++ b/parcels/particleset.py @@ -30,7 +30,6 @@ __all__ = ['ParticleSet'] - class ParticleSet(object): """Container class for storing particle and executing kernel over them. @@ -247,15 +246,6 @@ def cptr(self, bracket_index, slot_index): else: return None - # def allocate_cptrs(self): - # nparticles = 0 - # for sublist in self.plist: - # nparticles += sublist.shape[0] - # # == Allocate underlying data for C-allocated particles == # - # if self.ptype.uses_jit: - # self._particle_data = np.empty(nparticles, dtype=self.ptype.dtype) - # return nparticles - def get_list_array_index(self, pdata): """Searches for the list-array indices for a given particle. @@ -503,8 +493,9 @@ def size(self): def __repr__(self): return "\n".join([str(p) for sublist in self._plist for p in sublist]) - def __len__(self): - return self.size + # def __len__(self): + # return len(self._plist) + # # return self.size def __getitem__(self, key): bracket_index = 0 @@ -548,6 +539,9 @@ def add(self, particles): else: raise NotImplementedError('Only ParticleSets can be added to a ParticleSet') + # for sublist in particles.particles_c: + # self._plist_c.append(sublist) + def _merge_brackets_(self): lw_bound_nlist = int(self.nlist_limit / 2) nmerges = 1 From e0f995dcf5f9e80ea99f2cd435bca0d2e427969e Mon Sep 17 00:00:00 2001 From: Christian Kehl Date: Wed, 8 Jul 2020 20:26:39 +0200 Subject: [PATCH 47/56] list_of_pset_array_trials - atted __getitem__(self, key) func to ParticleSet for element iteration - slow, but working. --- parcels/particleset.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/parcels/particleset.py b/parcels/particleset.py index d2443f1db9..47f0664b4d 100644 --- a/parcels/particleset.py +++ b/parcels/particleset.py @@ -493,9 +493,9 @@ def size(self): def __repr__(self): return "\n".join([str(p) for sublist in self._plist for p in sublist]) - # def __len__(self): - # return len(self._plist) - # # return self.size + def __len__(self): + # return len(self._plist) + return self.size def __getitem__(self, key): bracket_index = 0 From 95959939f93169e2f9999666b506f2465eb9b349 Mon Sep 17 00:00:00 2001 From: Christian Kehl Date: Tue, 14 Jul 2020 17:38:01 +0200 Subject: [PATCH 48/56] list_of_pset_array_trials - fixed errors emerging from tests --- parcels/particle.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/parcels/particle.py b/parcels/particle.py index 1800b52b77..0762b46749 100644 --- a/parcels/particle.py +++ b/parcels/particle.py @@ -260,5 +260,5 @@ def __init__(self, *args, **kwargs): setattr(self, 'c'+index, getattr(self, index+'p').value) def update_cptr(self, new_ptr): - np.copyto(new_ptr, self._cptr) + # np.copyto(new_ptr, self._cptr) self._cptr = new_ptr From 613bda53e762f1c65d3f4009d12a2e3cc4519faf Mon Sep 17 00:00:00 2001 From: Christian Kehl Date: Tue, 14 Jul 2020 18:20:29 +0200 Subject: [PATCH 49/56] list_of_pset_array_trials - fixing linter errors --- parcels/particleset.py | 1 + parcels/particleset_benchmark.py | 76 ++++++++++++++++++++++++++++++++ 2 files changed, 77 insertions(+) diff --git a/parcels/particleset.py b/parcels/particleset.py index 47f0664b4d..ce44edd3cf 100644 --- a/parcels/particleset.py +++ b/parcels/particleset.py @@ -30,6 +30,7 @@ __all__ = ['ParticleSet'] + class ParticleSet(object): """Container class for storing particle and executing kernel over them. diff --git a/parcels/particleset_benchmark.py b/parcels/particleset_benchmark.py index ece52e391b..d71943c5ac 100644 --- a/parcels/particleset_benchmark.py +++ b/parcels/particleset_benchmark.py @@ -51,6 +51,82 @@ def measure_mem_usage(): USE_ASYNC_MEMLOG = False USE_RUSE_SYNC_MEMLOG = False # can be faulty +class ParticleSet_TimingLog(): + stime = 0 + etime = 0 + mtime = 0 + samples = [] + times_steps = [] + _iter = 0 + + def start_timing(self): + if MPI: + mpi_comm = MPI.COMM_WORLD + mpi_rank = mpi_comm.Get_rank() + if mpi_rank == 0: + # self.stime = MPI.Wtime() + # self.stime = time_module.perf_counter() + self.stime = time_module.process_time() + else: + self.stime = time_module.perf_counter() + + def stop_timing(self): + if MPI: + mpi_comm = MPI.COMM_WORLD + mpi_rank = mpi_comm.Get_rank() + if mpi_rank == 0: + # self.etime = MPI.Wtime() + # self.etime = time_module.perf_counter() + self.etime = time_module.process_time() + else: + self.etime = time_module.perf_counter() + + def accumulate_timing(self): + if MPI: + mpi_comm = MPI.COMM_WORLD + mpi_rank = mpi_comm.Get_rank() + if mpi_rank == 0: + self.mtime += (self.etime-self.stime) + else: + self.mtime = 0 + else: + self.mtime += (self.etime-self.stime) + + def advance_iteration(self): + if MPI: + mpi_comm = MPI.COMM_WORLD + mpi_rank = mpi_comm.Get_rank() + if mpi_rank == 0: + self.times_steps.append(self.mtime) + self.samples.append(self._iter) + self._iter += 1 + self.mtime = 0 + else: + self.times_steps.append(self.mtime) + self.samples.append(self._iter) + self._iter += 1 + self.mtime = 0 + + +class ParticleSet_ParamLogging(): + samples = [] + params = [] + _iter = 0 + + def advance_iteration(self, param): + if MPI: + mpi_comm = MPI.COMM_WORLD + mpi_rank = mpi_comm.Get_rank() + if mpi_rank == 0: + self.params.append(param) + self.samples.append(self._iter) + self._iter += 1 + else: + self.params.append(param) + self.samples.append(self._iter) + self._iter += 1 + + class ParticleSet_Benchmark(ParticleSet): def __init__(self, fieldset, pclass=JITParticle, lon=None, lat=None, depth=None, time=None, repeatdt=None, From 11f80275e23a1377a43456f8de762ab5cf72ec23 Mon Sep 17 00:00:00 2001 From: Christian Kehl Date: Thu, 16 Jul 2020 17:07:21 +0200 Subject: [PATCH 50/56] list_of_pset_array_trials - added benchmark refactoring. tested and working (stommel). --- parcels/particleset_benchmark.py | 71 ++++++++++++++++++++++++++------ performance/benchmark_stommel.py | 3 -- 2 files changed, 59 insertions(+), 15 deletions(-) diff --git a/parcels/particleset_benchmark.py b/parcels/particleset_benchmark.py index d71943c5ac..999f02a402 100644 --- a/parcels/particleset_benchmark.py +++ b/parcels/particleset_benchmark.py @@ -55,10 +55,35 @@ class ParticleSet_TimingLog(): stime = 0 etime = 0 mtime = 0 - samples = [] - times_steps = [] + _samples = [] + _times_steps = [] _iter = 0 + def __init__(self): + self.stime = 0 + self.etime = 0 + self.mtime = 0 + self._samples = [] + self._times_steps = [] + self._iter = 0 + + @property + def timing(self): + return self._times_steps + + @property + def samples(self): + return self._samples + + def __len__(self): + return len(self._samples) + + def get_values(self): + return self._times_steps + + def get_value(self, index): + return self._times_steps[index] + def start_timing(self): if MPI: mpi_comm = MPI.COMM_WORLD @@ -97,33 +122,55 @@ def advance_iteration(self): mpi_comm = MPI.COMM_WORLD mpi_rank = mpi_comm.Get_rank() if mpi_rank == 0: - self.times_steps.append(self.mtime) - self.samples.append(self._iter) + self._times_steps.append(self.mtime) + self._samples.append(self._iter) self._iter += 1 self.mtime = 0 else: - self.times_steps.append(self.mtime) - self.samples.append(self._iter) + self._times_steps.append(self.mtime) + self._samples.append(self._iter) self._iter += 1 self.mtime = 0 class ParticleSet_ParamLogging(): - samples = [] - params = [] + _samples = [] + _params = [] _iter = 0 + def __init__(self): + self._samples = [] + self._params = [] + self._iter = 0 + + @property + def samples(self): + return self._samples + + @property + def params(self): + return self._params + + def get_params(self): + return self._params + + def get_param(self, index): + return self._params[index] + + def __len__(self): + return len(self._samples) + def advance_iteration(self, param): if MPI: mpi_comm = MPI.COMM_WORLD mpi_rank = mpi_comm.Get_rank() if mpi_rank == 0: - self.params.append(param) - self.samples.append(self._iter) + self._params.append(param) + self._samples.append(self._iter) self._iter += 1 else: - self.params.append(param) - self.samples.append(self._iter) + self._params.append(param) + self._samples.append(self._iter) self._iter += 1 diff --git a/performance/benchmark_stommel.py b/performance/benchmark_stommel.py index 6bf63192c4..eb0e125b64 100644 --- a/performance/benchmark_stommel.py +++ b/performance/benchmark_stommel.py @@ -426,11 +426,8 @@ def Age(particle, fieldset, time): mpi_comm = MPI.COMM_WORLD mpi_rank = mpi_comm.Get_rank() if mpi_rank==0: - # endtime = ostime.time() - # endtime = MPI.Wtime() endtime = ostime.process_time() else: - #endtime = ostime.time() endtime = ostime.process_time() if args.write_out: From 8737b2f89b7383641ab357ffce0e170c820a2642 Mon Sep 17 00:00:00 2001 From: Christian Kehl Date: Mon, 20 Jul 2020 23:06:58 +0200 Subject: [PATCH 51/56] list_of_pset_array_trials - rebase to benchmarking branch rev 54e4c18 --- parcels/particleset_benchmark.py | 123 ------------------------------- performance/benchmark_stommel.py | 3 + 2 files changed, 3 insertions(+), 123 deletions(-) diff --git a/parcels/particleset_benchmark.py b/parcels/particleset_benchmark.py index 999f02a402..ece52e391b 100644 --- a/parcels/particleset_benchmark.py +++ b/parcels/particleset_benchmark.py @@ -51,129 +51,6 @@ def measure_mem_usage(): USE_ASYNC_MEMLOG = False USE_RUSE_SYNC_MEMLOG = False # can be faulty -class ParticleSet_TimingLog(): - stime = 0 - etime = 0 - mtime = 0 - _samples = [] - _times_steps = [] - _iter = 0 - - def __init__(self): - self.stime = 0 - self.etime = 0 - self.mtime = 0 - self._samples = [] - self._times_steps = [] - self._iter = 0 - - @property - def timing(self): - return self._times_steps - - @property - def samples(self): - return self._samples - - def __len__(self): - return len(self._samples) - - def get_values(self): - return self._times_steps - - def get_value(self, index): - return self._times_steps[index] - - def start_timing(self): - if MPI: - mpi_comm = MPI.COMM_WORLD - mpi_rank = mpi_comm.Get_rank() - if mpi_rank == 0: - # self.stime = MPI.Wtime() - # self.stime = time_module.perf_counter() - self.stime = time_module.process_time() - else: - self.stime = time_module.perf_counter() - - def stop_timing(self): - if MPI: - mpi_comm = MPI.COMM_WORLD - mpi_rank = mpi_comm.Get_rank() - if mpi_rank == 0: - # self.etime = MPI.Wtime() - # self.etime = time_module.perf_counter() - self.etime = time_module.process_time() - else: - self.etime = time_module.perf_counter() - - def accumulate_timing(self): - if MPI: - mpi_comm = MPI.COMM_WORLD - mpi_rank = mpi_comm.Get_rank() - if mpi_rank == 0: - self.mtime += (self.etime-self.stime) - else: - self.mtime = 0 - else: - self.mtime += (self.etime-self.stime) - - def advance_iteration(self): - if MPI: - mpi_comm = MPI.COMM_WORLD - mpi_rank = mpi_comm.Get_rank() - if mpi_rank == 0: - self._times_steps.append(self.mtime) - self._samples.append(self._iter) - self._iter += 1 - self.mtime = 0 - else: - self._times_steps.append(self.mtime) - self._samples.append(self._iter) - self._iter += 1 - self.mtime = 0 - - -class ParticleSet_ParamLogging(): - _samples = [] - _params = [] - _iter = 0 - - def __init__(self): - self._samples = [] - self._params = [] - self._iter = 0 - - @property - def samples(self): - return self._samples - - @property - def params(self): - return self._params - - def get_params(self): - return self._params - - def get_param(self, index): - return self._params[index] - - def __len__(self): - return len(self._samples) - - def advance_iteration(self, param): - if MPI: - mpi_comm = MPI.COMM_WORLD - mpi_rank = mpi_comm.Get_rank() - if mpi_rank == 0: - self._params.append(param) - self._samples.append(self._iter) - self._iter += 1 - else: - self._params.append(param) - self._samples.append(self._iter) - self._iter += 1 - - class ParticleSet_Benchmark(ParticleSet): def __init__(self, fieldset, pclass=JITParticle, lon=None, lat=None, depth=None, time=None, repeatdt=None, diff --git a/performance/benchmark_stommel.py b/performance/benchmark_stommel.py index eb0e125b64..6bf63192c4 100644 --- a/performance/benchmark_stommel.py +++ b/performance/benchmark_stommel.py @@ -426,8 +426,11 @@ def Age(particle, fieldset, time): mpi_comm = MPI.COMM_WORLD mpi_rank = mpi_comm.Get_rank() if mpi_rank==0: + # endtime = ostime.time() + # endtime = MPI.Wtime() endtime = ostime.process_time() else: + #endtime = ostime.time() endtime = ostime.process_time() if args.write_out: From f743a8aee8b2ff4db993edccad6f61788e820ad2 Mon Sep 17 00:00:00 2001 From: Christian Kehl Date: Tue, 7 Jul 2020 20:02:48 +0200 Subject: [PATCH 52/56] list_of_pset_array_trials - initial addenda for list-of-particles version --- parcels/particleset.py | 1 - 1 file changed, 1 deletion(-) diff --git a/parcels/particleset.py b/parcels/particleset.py index ce44edd3cf..47f0664b4d 100644 --- a/parcels/particleset.py +++ b/parcels/particleset.py @@ -30,7 +30,6 @@ __all__ = ['ParticleSet'] - class ParticleSet(object): """Container class for storing particle and executing kernel over them. From 75de5e4f0d3981eecf7c7f67848615fad25ee787 Mon Sep 17 00:00:00 2001 From: Christian Kehl Date: Tue, 14 Jul 2020 18:20:29 +0200 Subject: [PATCH 53/56] list_of_pset_array_trials - fixing linter errors --- parcels/particleset.py | 1 + parcels/particleset_benchmark.py | 76 ++++++++++++++++++++++++++++++++ 2 files changed, 77 insertions(+) diff --git a/parcels/particleset.py b/parcels/particleset.py index 47f0664b4d..ce44edd3cf 100644 --- a/parcels/particleset.py +++ b/parcels/particleset.py @@ -30,6 +30,7 @@ __all__ = ['ParticleSet'] + class ParticleSet(object): """Container class for storing particle and executing kernel over them. diff --git a/parcels/particleset_benchmark.py b/parcels/particleset_benchmark.py index ece52e391b..d71943c5ac 100644 --- a/parcels/particleset_benchmark.py +++ b/parcels/particleset_benchmark.py @@ -51,6 +51,82 @@ def measure_mem_usage(): USE_ASYNC_MEMLOG = False USE_RUSE_SYNC_MEMLOG = False # can be faulty +class ParticleSet_TimingLog(): + stime = 0 + etime = 0 + mtime = 0 + samples = [] + times_steps = [] + _iter = 0 + + def start_timing(self): + if MPI: + mpi_comm = MPI.COMM_WORLD + mpi_rank = mpi_comm.Get_rank() + if mpi_rank == 0: + # self.stime = MPI.Wtime() + # self.stime = time_module.perf_counter() + self.stime = time_module.process_time() + else: + self.stime = time_module.perf_counter() + + def stop_timing(self): + if MPI: + mpi_comm = MPI.COMM_WORLD + mpi_rank = mpi_comm.Get_rank() + if mpi_rank == 0: + # self.etime = MPI.Wtime() + # self.etime = time_module.perf_counter() + self.etime = time_module.process_time() + else: + self.etime = time_module.perf_counter() + + def accumulate_timing(self): + if MPI: + mpi_comm = MPI.COMM_WORLD + mpi_rank = mpi_comm.Get_rank() + if mpi_rank == 0: + self.mtime += (self.etime-self.stime) + else: + self.mtime = 0 + else: + self.mtime += (self.etime-self.stime) + + def advance_iteration(self): + if MPI: + mpi_comm = MPI.COMM_WORLD + mpi_rank = mpi_comm.Get_rank() + if mpi_rank == 0: + self.times_steps.append(self.mtime) + self.samples.append(self._iter) + self._iter += 1 + self.mtime = 0 + else: + self.times_steps.append(self.mtime) + self.samples.append(self._iter) + self._iter += 1 + self.mtime = 0 + + +class ParticleSet_ParamLogging(): + samples = [] + params = [] + _iter = 0 + + def advance_iteration(self, param): + if MPI: + mpi_comm = MPI.COMM_WORLD + mpi_rank = mpi_comm.Get_rank() + if mpi_rank == 0: + self.params.append(param) + self.samples.append(self._iter) + self._iter += 1 + else: + self.params.append(param) + self.samples.append(self._iter) + self._iter += 1 + + class ParticleSet_Benchmark(ParticleSet): def __init__(self, fieldset, pclass=JITParticle, lon=None, lat=None, depth=None, time=None, repeatdt=None, From 1981d446ebb6992a4d874175165e309f86d63b8a Mon Sep 17 00:00:00 2001 From: Christian Kehl Date: Thu, 16 Jul 2020 17:07:21 +0200 Subject: [PATCH 54/56] list_of_pset_array_trials - added benchmark refactoring. tested and working (stommel). --- parcels/particleset_benchmark.py | 71 ++++++++++++++++++++++++++------ performance/benchmark_stommel.py | 3 -- 2 files changed, 59 insertions(+), 15 deletions(-) diff --git a/parcels/particleset_benchmark.py b/parcels/particleset_benchmark.py index d71943c5ac..999f02a402 100644 --- a/parcels/particleset_benchmark.py +++ b/parcels/particleset_benchmark.py @@ -55,10 +55,35 @@ class ParticleSet_TimingLog(): stime = 0 etime = 0 mtime = 0 - samples = [] - times_steps = [] + _samples = [] + _times_steps = [] _iter = 0 + def __init__(self): + self.stime = 0 + self.etime = 0 + self.mtime = 0 + self._samples = [] + self._times_steps = [] + self._iter = 0 + + @property + def timing(self): + return self._times_steps + + @property + def samples(self): + return self._samples + + def __len__(self): + return len(self._samples) + + def get_values(self): + return self._times_steps + + def get_value(self, index): + return self._times_steps[index] + def start_timing(self): if MPI: mpi_comm = MPI.COMM_WORLD @@ -97,33 +122,55 @@ def advance_iteration(self): mpi_comm = MPI.COMM_WORLD mpi_rank = mpi_comm.Get_rank() if mpi_rank == 0: - self.times_steps.append(self.mtime) - self.samples.append(self._iter) + self._times_steps.append(self.mtime) + self._samples.append(self._iter) self._iter += 1 self.mtime = 0 else: - self.times_steps.append(self.mtime) - self.samples.append(self._iter) + self._times_steps.append(self.mtime) + self._samples.append(self._iter) self._iter += 1 self.mtime = 0 class ParticleSet_ParamLogging(): - samples = [] - params = [] + _samples = [] + _params = [] _iter = 0 + def __init__(self): + self._samples = [] + self._params = [] + self._iter = 0 + + @property + def samples(self): + return self._samples + + @property + def params(self): + return self._params + + def get_params(self): + return self._params + + def get_param(self, index): + return self._params[index] + + def __len__(self): + return len(self._samples) + def advance_iteration(self, param): if MPI: mpi_comm = MPI.COMM_WORLD mpi_rank = mpi_comm.Get_rank() if mpi_rank == 0: - self.params.append(param) - self.samples.append(self._iter) + self._params.append(param) + self._samples.append(self._iter) self._iter += 1 else: - self.params.append(param) - self.samples.append(self._iter) + self._params.append(param) + self._samples.append(self._iter) self._iter += 1 diff --git a/performance/benchmark_stommel.py b/performance/benchmark_stommel.py index 6bf63192c4..eb0e125b64 100644 --- a/performance/benchmark_stommel.py +++ b/performance/benchmark_stommel.py @@ -426,11 +426,8 @@ def Age(particle, fieldset, time): mpi_comm = MPI.COMM_WORLD mpi_rank = mpi_comm.Get_rank() if mpi_rank==0: - # endtime = ostime.time() - # endtime = MPI.Wtime() endtime = ostime.process_time() else: - #endtime = ostime.time() endtime = ostime.process_time() if args.write_out: From eb329ff965d672c0522e850246b07349b779f248 Mon Sep 17 00:00:00 2001 From: Christian Kehl Date: Mon, 20 Jul 2020 23:06:58 +0200 Subject: [PATCH 55/56] list_of_pset_array_trials - rebase to benchmarking branch rev 54e4c18 --- parcels/particleset_benchmark.py | 123 ------------------------------- performance/benchmark_stommel.py | 3 + 2 files changed, 3 insertions(+), 123 deletions(-) diff --git a/parcels/particleset_benchmark.py b/parcels/particleset_benchmark.py index 999f02a402..ece52e391b 100644 --- a/parcels/particleset_benchmark.py +++ b/parcels/particleset_benchmark.py @@ -51,129 +51,6 @@ def measure_mem_usage(): USE_ASYNC_MEMLOG = False USE_RUSE_SYNC_MEMLOG = False # can be faulty -class ParticleSet_TimingLog(): - stime = 0 - etime = 0 - mtime = 0 - _samples = [] - _times_steps = [] - _iter = 0 - - def __init__(self): - self.stime = 0 - self.etime = 0 - self.mtime = 0 - self._samples = [] - self._times_steps = [] - self._iter = 0 - - @property - def timing(self): - return self._times_steps - - @property - def samples(self): - return self._samples - - def __len__(self): - return len(self._samples) - - def get_values(self): - return self._times_steps - - def get_value(self, index): - return self._times_steps[index] - - def start_timing(self): - if MPI: - mpi_comm = MPI.COMM_WORLD - mpi_rank = mpi_comm.Get_rank() - if mpi_rank == 0: - # self.stime = MPI.Wtime() - # self.stime = time_module.perf_counter() - self.stime = time_module.process_time() - else: - self.stime = time_module.perf_counter() - - def stop_timing(self): - if MPI: - mpi_comm = MPI.COMM_WORLD - mpi_rank = mpi_comm.Get_rank() - if mpi_rank == 0: - # self.etime = MPI.Wtime() - # self.etime = time_module.perf_counter() - self.etime = time_module.process_time() - else: - self.etime = time_module.perf_counter() - - def accumulate_timing(self): - if MPI: - mpi_comm = MPI.COMM_WORLD - mpi_rank = mpi_comm.Get_rank() - if mpi_rank == 0: - self.mtime += (self.etime-self.stime) - else: - self.mtime = 0 - else: - self.mtime += (self.etime-self.stime) - - def advance_iteration(self): - if MPI: - mpi_comm = MPI.COMM_WORLD - mpi_rank = mpi_comm.Get_rank() - if mpi_rank == 0: - self._times_steps.append(self.mtime) - self._samples.append(self._iter) - self._iter += 1 - self.mtime = 0 - else: - self._times_steps.append(self.mtime) - self._samples.append(self._iter) - self._iter += 1 - self.mtime = 0 - - -class ParticleSet_ParamLogging(): - _samples = [] - _params = [] - _iter = 0 - - def __init__(self): - self._samples = [] - self._params = [] - self._iter = 0 - - @property - def samples(self): - return self._samples - - @property - def params(self): - return self._params - - def get_params(self): - return self._params - - def get_param(self, index): - return self._params[index] - - def __len__(self): - return len(self._samples) - - def advance_iteration(self, param): - if MPI: - mpi_comm = MPI.COMM_WORLD - mpi_rank = mpi_comm.Get_rank() - if mpi_rank == 0: - self._params.append(param) - self._samples.append(self._iter) - self._iter += 1 - else: - self._params.append(param) - self._samples.append(self._iter) - self._iter += 1 - - class ParticleSet_Benchmark(ParticleSet): def __init__(self, fieldset, pclass=JITParticle, lon=None, lat=None, depth=None, time=None, repeatdt=None, diff --git a/performance/benchmark_stommel.py b/performance/benchmark_stommel.py index eb0e125b64..6bf63192c4 100644 --- a/performance/benchmark_stommel.py +++ b/performance/benchmark_stommel.py @@ -426,8 +426,11 @@ def Age(particle, fieldset, time): mpi_comm = MPI.COMM_WORLD mpi_rank = mpi_comm.Get_rank() if mpi_rank==0: + # endtime = ostime.time() + # endtime = MPI.Wtime() endtime = ostime.process_time() else: + #endtime = ostime.time() endtime = ostime.process_time() if args.write_out: From ae1e349ec9ddc5e1b834e238717f8280067ae52a Mon Sep 17 00:00:00 2001 From: Christian Kehl Date: Thu, 4 Feb 2021 22:16:48 +0100 Subject: [PATCH 56/56] list_of_pset_array_trials - changed the branch names in the performance scripts --- performance/benchmark_CMEMS.py | 2 +- performance/benchmark_bickleyjet.py | 2 +- performance/benchmark_deep_migration_NPacific.py | 2 +- performance/benchmark_doublegyre.py | 2 +- performance/benchmark_galapagos_backwards.py | 2 +- performance/benchmark_palaeo_Y2K.py | 2 +- performance/benchmark_perlin.py | 2 +- performance/benchmark_stommel.py | 2 +- 8 files changed, 8 insertions(+), 8 deletions(-) diff --git a/performance/benchmark_CMEMS.py b/performance/benchmark_CMEMS.py index af6af774b0..a1e6e6c4c9 100644 --- a/performance/benchmark_CMEMS.py +++ b/performance/benchmark_CMEMS.py @@ -153,7 +153,7 @@ def perIterGC(): nowtime = datetime.now() random.seed(nowtime.microsecond) - branch = "benchmarking" + branch = "LoA" computer_env = "local/unspecified" scenario = "CMEMS" headdir = "" diff --git a/performance/benchmark_bickleyjet.py b/performance/benchmark_bickleyjet.py index 7d04d36d38..0a3718231e 100644 --- a/performance/benchmark_bickleyjet.py +++ b/performance/benchmark_bickleyjet.py @@ -230,7 +230,7 @@ def Age(particle, fieldset, time): nowtime = datetime.datetime.now() random.seed(nowtime.microsecond) - branch = "benchmarking" + branch = "LoA" computer_env = "local/unspecified" scenario = "bickleyjet" odir = "" diff --git a/performance/benchmark_deep_migration_NPacific.py b/performance/benchmark_deep_migration_NPacific.py index 2edff5bb07..8c3c763063 100644 --- a/performance/benchmark_deep_migration_NPacific.py +++ b/performance/benchmark_deep_migration_NPacific.py @@ -240,7 +240,7 @@ def Profiles(particle, fieldset, time): time_in_days = int(float(eval(args.time_in_days))) with_GC = args.useGC - branch = "benchmarking" + branch = "LoA" computer_env = "local/unspecified" scenario = "deep_migration" headdir = "" diff --git a/performance/benchmark_doublegyre.py b/performance/benchmark_doublegyre.py index 8733e6a02f..185fd06409 100644 --- a/performance/benchmark_doublegyre.py +++ b/performance/benchmark_doublegyre.py @@ -224,7 +224,7 @@ def Age(particle, fieldset, time): nowtime = datetime.datetime.now() random.seed(nowtime.microsecond) - branch = "benchmarking" + branch = "LoA" computer_env = "local/unspecified" scenario = "doublegyre" odir = "" diff --git a/performance/benchmark_galapagos_backwards.py b/performance/benchmark_galapagos_backwards.py index be53489051..ba6ea91f68 100644 --- a/performance/benchmark_galapagos_backwards.py +++ b/performance/benchmark_galapagos_backwards.py @@ -107,7 +107,7 @@ def periodicBC(particle, fieldSet, time): time_in_days = int(float(eval(args.time_in_days))) with_GC = args.useGC - branch = "benchmarking" + branch = "LoA" computer_env = "local/unspecified" scenario = "galapagos_backwards" headdir = "" diff --git a/performance/benchmark_palaeo_Y2K.py b/performance/benchmark_palaeo_Y2K.py index 981b4b9bd8..ed138738c4 100644 --- a/performance/benchmark_palaeo_Y2K.py +++ b/performance/benchmark_palaeo_Y2K.py @@ -263,7 +263,7 @@ class DinoParticle(JITParticle): time_in_years = int(time_in_days/366.0) with_GC = args.useGC - branch = "benchmarking" + branch = "LoA" computer_env = "local/unspecified" scenario = "palaeo-parcels" headdir = "" diff --git a/performance/benchmark_perlin.py b/performance/benchmark_perlin.py index e2a8246fc9..0d6eb96d2f 100644 --- a/performance/benchmark_perlin.py +++ b/performance/benchmark_perlin.py @@ -239,7 +239,7 @@ def Age(particle, fieldset, time): nowtime = datetime.datetime.now() random.seed(nowtime.microsecond) - branch = "benchmarking" + branch = "LoA" computer_env = "local/unspecified" scenario = "perlin" odir = "" diff --git a/performance/benchmark_stommel.py b/performance/benchmark_stommel.py index 6bf63192c4..d2ecccc158 100644 --- a/performance/benchmark_stommel.py +++ b/performance/benchmark_stommel.py @@ -257,7 +257,7 @@ def Age(particle, fieldset, time): nowtime = datetime.datetime.now() random.seed(nowtime.microsecond) - branch = "benchmarking" + branch = "LoA" computer_env = "local/unspecified" scenario = "stommel" odir = ""