Support for creating a ParticleSet with empty FieldSet

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]

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:

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:

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,19 +150,19 @@ 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
 
     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):