diff --git a/parcels/codegenerator.py b/parcels/codegenerator.py index d22bfd478d..d7785d468e 100644 --- a/parcels/codegenerator.py +++ b/parcels/codegenerator.py @@ -12,6 +12,7 @@ from parcels.field import SummedField from parcels.field import VectorField from parcels.grid import Grid +from parcels.particle import JITParticle from parcels.tools.loggers import logger @@ -218,7 +219,7 @@ class IntrinsicTransformer(ast.NodeTransformer): names, such as 'particle' or 'fieldset', inserts placeholder objects and propagates attribute access.""" - def __init__(self, fieldset, ptype): + def __init__(self, fieldset=None, ptype=JITParticle): self.fieldset = fieldset self.ptype = ptype @@ -237,7 +238,7 @@ def get_tmp(self): def visit_Name(self, node): """Inject IntrinsicNode objects into the tree according to keyword""" - if node.id == 'fieldset': + if node.id == 'fieldset' and self.fieldset is not None: node = FieldSetNode(self.fieldset, ccode='fset') elif node.id == 'particle': node = ParticleNode(self.ptype, ccode='particles') @@ -382,7 +383,7 @@ class KernelGenerator(ast.NodeVisitor): kernel_vars = ['particle', 'fieldset', 'time', 'output_time', 'tol'] array_vars = [] - def __init__(self, fieldset, ptype): + def __init__(self, fieldset=None, ptype=JITParticle): self.fieldset = fieldset self.ptype = ptype self.field_args = collections.OrderedDict() @@ -901,7 +902,7 @@ def generate(self, funcname, field_args, const_args, kernel_ast, c_include): ccode += [str(c.Include("math.h", system=False))] ccode += [str(c.Assign('double _next_dt', '0'))] ccode += [str(c.Assign('size_t _next_dt_set', '0'))] - ccode += [str(c.Assign('const int ngrid', str(self.fieldset.gridset.size)))] + ccode += [str(c.Assign('const int ngrid', str(self.fieldset.gridset.size if self.fieldset is not None else 1)))] # ==== Generate type definition for particle type ==== # vdeclp = [c.Pointer(c.POD(v.dtype, v.name)) for v in self.ptype.variables] diff --git a/parcels/kernel.py b/parcels/kernel.py index a8c67bc1c1..491be7e51b 100644 --- a/parcels/kernel.py +++ b/parcels/kernel.py @@ -230,35 +230,36 @@ def load_lib(self): def execute_jit(self, pset, endtime, dt): """Invokes JIT engine to perform the core update loop""" - if len(pset) > 0 and pset.particle_data['xi'].ndim == 2: + if len(pset) > 0 and pset.particle_data['xi'].ndim == 2 and pset.fieldset is not None: assert pset.fieldset.gridset.size == pset.particle_data['xi'].shape[1], \ 'FieldSet has different number 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 - # C-contiguous memory layout for JIT mode. - for f in pset.fieldset.get_fields(): - if type(f) in [VectorField, NestedField, SummedField]: - continue - if f in self.field_args.values(): - f.chunk_data() - else: - for block_id in range(len(f.data_chunks)): - f.data_chunks[block_id] = None - f.c_data_chunks[block_id] = None - - for g in pset.fieldset.gridset.grids: - g.load_chunk = np.where(g.load_chunk == 1, 2, g.load_chunk) - if len(g.load_chunk) > 0: # not the case if a field in not called in the kernel - if not g.load_chunk.flags.c_contiguous: - g.load_chunk = g.load_chunk.copy() - if not g.depth.flags.c_contiguous: - g.depth = g.depth.copy() - if not g.lon.flags.c_contiguous: - g.lon = g.lon.copy() - if not g.lat.flags.c_contiguous: - g.lat = g.lat.copy() + if pset.fieldset is not None: + 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 + # C-contiguous memory layout for JIT mode. + for f in pset.fieldset.get_fields(): + if type(f) in [VectorField, NestedField, SummedField]: + continue + if f in self.field_args.values(): + f.chunk_data() + else: + for block_id in range(len(f.data_chunks)): + f.data_chunks[block_id] = None + f.c_data_chunks[block_id] = None + + for g in pset.fieldset.gridset.grids: + g.load_chunk = np.where(g.load_chunk == 1, 2, g.load_chunk) + if len(g.load_chunk) > 0: # not the case if a field in not called in the kernel + if not g.load_chunk.flags.c_contiguous: + g.load_chunk = g.load_chunk.copy() + if not g.depth.flags.c_contiguous: + g.depth = g.depth.copy() + if not g.lon.flags.c_contiguous: + g.lon = g.lon.copy() + if not g.lat.flags.c_contiguous: + g.lat = g.lat.copy() fargs = [byref(f.ctypes_struct) for f in self.field_args.values()] fargs += [c_double(f) for f in self.const_args.values()] @@ -283,10 +284,11 @@ def execute_python(self, pset, endtime, dt): # back up variables in case of ErrorCode.Repeat p_var_back = {} - for f in self.fieldset.get_fields(): - if type(f) in [VectorField, NestedField, SummedField]: - continue - f.data = np.array(f.data) + if self.fieldset is not None: + for f in self.fieldset.get_fields(): + if type(f) in [VectorField, NestedField, SummedField]: + continue + f.data = np.array(f.data) for p in range(pset.size): particles.set_index(p) @@ -390,9 +392,10 @@ def remove_deleted(pset): recovery_map = recovery_base_map.copy() recovery_map.update(recovery) - for g in pset.fieldset.gridset.grids: - if len(g.load_chunk) > 0: # not the case if a field in not called in the kernel - g.load_chunk = np.where(g.load_chunk == 2, 3, g.load_chunk) + if pset.fieldset is not None: + for g in pset.fieldset.gridset.grids: + if len(g.load_chunk) > 0: # not the case if a field in not called in the kernel + g.load_chunk = np.where(g.load_chunk == 2, 3, g.load_chunk) # Execute the kernel over the particle set if self.ptype.uses_jit: diff --git a/parcels/particleset.py b/parcels/particleset.py index eec7cf78d7..259b1b939e 100644 --- a/parcels/particleset.py +++ b/parcels/particleset.py @@ -96,7 +96,8 @@ class ParticleSet(object): Please note that this currently only supports fixed size particle sets. - :param fieldset: :mod:`parcels.fieldset.FieldSet` object from which to sample velocity + :param fieldset: :mod:`parcels.fieldset.FieldSet` object from which to sample velocity. + While fieldset=None is supported, this will throw a warning as it breaks most Parcels functionality :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 @@ -112,9 +113,13 @@ class ParticleSet(object): 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): + def __init__(self, fieldset=None, 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() + if self.fieldset is None: + logger.warning_once("No FieldSet provided in ParticleSet generation. " + "This breaks most Parcels functionality") + else: + self.fieldset.check_complete() partitions = kwargs.pop('partitions', None) def convert_to_array(var): @@ -133,7 +138,7 @@ def convert_to_array(var): pid = pid_orig + pclass.lastID if depth is None: - mindepth, _ = self.fieldset.gridset.dimrange('depth') + mindepth = self.fieldset.gridset.dimrange('depth')[0] if self.fieldset is not None else 0 depth = np.ones(lon.size) * mindepth else: depth = convert_to_array(depth) @@ -150,7 +155,7 @@ def _convert_to_reltime(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 + self.time_origin = fieldset.time_origin if self.fieldset is not None else 0 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 _convert_to_reltime(t) else t for t in time]) @@ -212,7 +217,10 @@ def _convert_to_reltime(time): pclass.setLastID(offset+1) if lonlatdepth_dtype is None: - self.lonlatdepth_dtype = self.lonlatdepth_dtype_from_field_interp_method(fieldset.U) + if fieldset is not None: + self.lonlatdepth_dtype = self.lonlatdepth_dtype_from_field_interp_method(fieldset.U) + else: + self.lonlatdepth_dtype = np.float32 else: self.lonlatdepth_dtype = lonlatdepth_dtype assert self.lonlatdepth_dtype in [np.float32, np.float64], \ @@ -227,7 +235,8 @@ def _convert_to_reltime(time): initialised = set() for v in self.ptype.variables: if v.name in ['xi', 'yi', 'zi', 'ti']: - self.particle_data[v.name] = np.empty((len(lon), fieldset.gridset.size), dtype=v.dtype) + ngrid = fieldset.gridset.size if fieldset is not None else 1 + self.particle_data[v.name] = np.empty((len(lon), ngrid), dtype=v.dtype) else: self.particle_data[v.name] = np.empty(len(lon), dtype=v.dtype) @@ -608,7 +617,7 @@ def execute(self, pyfunc=AdvectionRK4, endtime=None, runtime=None, dt=1., assert outputdt is None or outputdt >= 0, 'outputdt must be positive' assert moviedt is None or moviedt >= 0, 'moviedt must be positive' - mintime, maxtime = self.fieldset.gridset.dimrange('time_full') + mintime, maxtime = self.fieldset.gridset.dimrange('time_full') if self.fieldset is not None else (0, 1) if np.any(np.isnan(self.particle_data['time'])): self.particle_data['time'][np.isnan(self.particle_data['time'])] = mintime if dt >= 0 else maxtime @@ -656,7 +665,7 @@ def execute(self, pyfunc=AdvectionRK4, endtime=None, runtime=None, dt=1., 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)) + next_input = self.fieldset.computeTimeChunk(time, np.sign(dt)) if self.fieldset is not None else np.inf tol = 1e-12 if verbose_progress is None: diff --git a/tests/test_kernel_language.py b/tests/test_kernel_language.py index 374aeee6ad..45cb4b42c0 100644 --- a/tests/test_kernel_language.py +++ b/tests/test_kernel_language.py @@ -42,11 +42,11 @@ def fieldset_fixture(xdim=20, ydim=20): ('Mul', '3 * 5', 15), ('Div', '24 / 4', 6), ]) -def test_expression_int(fieldset, mode, name, expr, result, npart=10): +def test_expression_int(mode, name, expr, result, npart=10): """ Test basic arithmetic expressions """ class TestParticle(ptype[mode]): p = Variable('p', dtype=np.float32) - pset = ParticleSet(fieldset, pclass=TestParticle, + pset = ParticleSet(None, pclass=TestParticle, lon=np.linspace(0., 1., npart), lat=np.zeros(npart) + 0.5) pset.execute(expr_kernel('Test%s' % name, pset, expr), endtime=1., dt=1.) @@ -61,11 +61,11 @@ class TestParticle(ptype[mode]): ('Div', '24. / 4.', 6), ('Pow', '2 ** 3', 8), ]) -def test_expression_float(fieldset, mode, name, expr, result, npart=10): +def test_expression_float(mode, name, expr, result, npart=10): """ Test basic arithmetic expressions """ class TestParticle(ptype[mode]): p = Variable('p', dtype=np.float32) - pset = ParticleSet(fieldset, pclass=TestParticle, + pset = ParticleSet(None, pclass=TestParticle, lon=np.linspace(0., 1., npart), lat=np.zeros(npart) + 0.5) pset.execute(expr_kernel('Test%s' % name, pset, expr), endtime=1., dt=1.) @@ -84,11 +84,11 @@ class TestParticle(ptype[mode]): ('Greater', '4 > 2', True), ('GreaterEq', '2 >= 4', False), ]) -def test_expression_bool(fieldset, mode, name, expr, result, npart=10): +def test_expression_bool(mode, name, expr, result, npart=10): """ Test basic arithmetic expressions """ class TestParticle(ptype[mode]): p = Variable('p', dtype=np.float32) - pset = ParticleSet(fieldset, pclass=TestParticle, + pset = ParticleSet(None, pclass=TestParticle, lon=np.linspace(0., 1., npart), lat=np.zeros(npart) + 0.5) pset.execute(expr_kernel('Test%s' % name, pset, expr), endtime=1., dt=1.) @@ -99,11 +99,11 @@ class TestParticle(ptype[mode]): @pytest.mark.parametrize('mode', ['scipy', 'jit']) -def test_while_if_break(fieldset, mode): +def test_while_if_break(mode): """Test while, if and break commands""" class TestParticle(ptype[mode]): p = Variable('p', dtype=np.float32, initial=0.) - pset = ParticleSet(fieldset, pclass=TestParticle, lon=[0], lat=[0]) + pset = ParticleSet(pclass=TestParticle, lon=[0], lat=[0]) def kernel(particle, fieldset, time): while particle.p < 30: @@ -117,12 +117,12 @@ def kernel(particle, fieldset, time): @pytest.mark.parametrize('mode', ['scipy', 'jit']) -def test_nested_if(fieldset, mode): +def test_nested_if(mode): """Test nested if commands""" class TestParticle(ptype[mode]): p0 = Variable('p0', dtype=np.int32, initial=0) p1 = Variable('p1', dtype=np.int32, initial=1) - pset = ParticleSet(fieldset, pclass=TestParticle, lon=0, lat=0) + pset = ParticleSet(pclass=TestParticle, lon=0, lat=0) def kernel(particle, fieldset, time): if particle.p1 >= particle.p0: @@ -135,11 +135,11 @@ def kernel(particle, fieldset, time): @pytest.mark.parametrize('mode', ['scipy', 'jit']) -def test_pass(fieldset, mode): +def test_pass(mode): """Test pass commands""" class TestParticle(ptype[mode]): p = Variable('p', dtype=np.int32, initial=0) - pset = ParticleSet(fieldset, pclass=TestParticle, lon=0, lat=0) + pset = ParticleSet(pclass=TestParticle, lon=0, lat=0) def kernel(particle, fieldset, time): particle.p = -1 @@ -150,8 +150,8 @@ def kernel(particle, fieldset, time): @pytest.mark.parametrize('mode', ['scipy', 'jit']) -def test_dt_as_variable_in_kernel(fieldset, mode): - pset = ParticleSet(fieldset, pclass=ptype[mode], lon=0, lat=0) +def test_dt_as_variable_in_kernel(mode): + pset = ParticleSet(pclass=ptype[mode], lon=0, lat=0) def kernel(particle, fieldset, time): dt = 1. # noqa @@ -159,10 +159,10 @@ def kernel(particle, fieldset, time): pset.execute(kernel, endtime=10, dt=1.) -def test_parcels_tmpvar_in_kernel(fieldset): - """Tests for error thrown if vartiable with 'tmp' defined in custom kernel""" +def test_parcels_tmpvar_in_kernel(): + """Tests for error thrown if variable with 'tmp' defined in custom kernel""" error_thrown = False - pset = ParticleSet(fieldset, pclass=JITParticle, lon=0, lat=0) + pset = ParticleSet(pclass=JITParticle, lon=0, lat=0) def kernel_tmpvar(particle, fieldset, time): parcels_tmpvar0 = 0 # noqa @@ -277,11 +277,11 @@ def random_series(npart, rngfunc, rngargs, mode): ('uniform', [0., 20.]), ('randint', [0, 20]), ]) -def test_random_float(fieldset, mode, rngfunc, rngargs, npart=10): +def test_random_float(mode, rngfunc, rngargs, npart=10): """ Test basic random number generation """ class TestParticle(ptype[mode]): p = Variable('p', dtype=np.float32 if rngfunc == 'randint' else np.float32) - pset = ParticleSet(fieldset, pclass=TestParticle, + pset = ParticleSet(pclass=TestParticle, lon=np.linspace(0., 1., npart), lat=np.zeros(npart) + 0.5) series = random_series(npart, rngfunc, rngargs, mode) @@ -331,10 +331,10 @@ def pykernel(particle, fieldset, time): @pytest.mark.parametrize('mode', ['scipy', 'jit']) -def test_dt_modif_by_kernel(fieldset, mode): +def test_dt_modif_by_kernel(mode): class TestParticle(ptype[mode]): age = Variable('age', dtype=np.float32) - pset = ParticleSet(fieldset, pclass=TestParticle, lon=[0.5], lat=[0]) + pset = ParticleSet(pclass=TestParticle, lon=[0.5], lat=[0]) def modif_dt(particle, fieldset, time): particle.age += particle.dt @@ -347,8 +347,8 @@ def modif_dt(particle, fieldset, time): @pytest.mark.parametrize('mode', ['scipy', 'jit']) @pytest.mark.parametrize('dt', [1e-2, 1e-6]) -def test_small_dt(fieldset, mode, dt, npart=10): - pset = ParticleSet(fieldset, pclass=ptype[mode], lon=np.zeros(npart), +def test_small_dt(mode, dt, npart=10): + pset = ParticleSet(pclass=ptype[mode], lon=np.zeros(npart), lat=np.zeros(npart), time=np.arange(0, npart)*dt*10) def DoNothing(particle, fieldset, time):