Merge branch 'master' into fix_minor_spelling_grammar_tutorials

Author: michaeldenes

.github/workflows/pypi-release.yaml .github/workflows/pypi-release.yml.github/workflows/pypi-release.yaml renamed to .github/workflows/pypi-release.yml

@@ -94,3 +94,17 @@ jobs:
           user: __token__
           password: ${{ secrets.PARCELS_PYPI_PROD_TOKEN }}
           verbose: true
+
+  test-pypi-release:
+    needs: upload-to-pypi
+    runs-on: ubuntu-latest
+    steps:
+      - uses: conda-incubator/setup-miniconda@v3
+        with:
+          activate-environment: parcels
+          python-version: 3.8
+          channels: conda-forge
+      - run: conda install -c conda-forge c-compiler pip
+      - run: pip install parcels --no-cache
+      - run: curl https://raw.githubusercontent.com/OceanParcels/parcels/master/docs/examples/example_peninsula.py > example_peninsula.py
+      - run: python example_peninsula.py

docs/examples/example_globcurrent.py

@@ -173,14 +173,6 @@ def test_globcurrent_time_extrapolation_error(mode, use_xarray):
     pset.execute(AdvectionRK4, runtime=delta(days=1), dt=delta(minutes=5))
 
 
-@pytest.mark.parametrize('mode', ['scipy', 'jit'])
-@pytest.mark.parametrize('use_xarray', [True, False])
-def test_globcurrent_dt0(mode, use_xarray):
-    fieldset = set_globcurrent_fieldset(use_xarray=use_xarray)
-    pset = ParticleSet(fieldset, pclass=ptype[mode], lon=[25], lat=[-35])
-    pset.execute(AdvectionRK4, dt=0.)
-
-
 @pytest.mark.parametrize('mode', ['scipy', 'jit'])
 @pytest.mark.parametrize('dt', [-300, 300])
 @pytest.mark.parametrize('with_starttime', [True, False])
@@ -260,3 +252,22 @@ def test_globcurrent_pset_fromfile(mode, dt, pid_offset, tmpdir):
 
     for var in ['lon', 'lat', 'depth', 'time', 'id']:
         assert np.allclose([getattr(p, var) for p in pset], [getattr(p, var) for p in pset_new])
+
+
+@pytest.mark.parametrize('mode', ['scipy', 'jit'])
+def test_error_outputdt_not_multiple_dt(mode, tmpdir):
+    # Test that outputdt is a multiple of dt
+    fieldset = set_globcurrent_fieldset()
+
+    filepath = tmpdir.join("pfile_error_outputdt_not_multiple_dt.zarr")
+
+    dt = 81.2584344538292  # number for which output writing fails
+
+    pset = ParticleSet(fieldset, pclass=ptype[mode], lon=[0], lat=[0])
+    ofile = pset.ParticleFile(name=filepath, outputdt=delta(days=1))
+
+    def DoNothing(particle, fieldset, time):
+        pass
+
+    with pytest.raises(ValueError):
+        pset.execute(DoNothing, runtime=delta(days=10), dt=dt, output_file=ofile)

docs/examples/tutorial_jit_vs_scipy.ipynb

@@ -5,7 +5,7 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "# JIT Particles and Scipy particles\n"
+    "# JIT-vs-Scipy Particles"
    ]
   },
   {

docs/examples/tutorial_parcels_structure.ipynb

@@ -302,7 +302,7 @@
     "    ```python\n",
     "    value = fieldset.U[time, particle.depth, particle.lat, particle.lon]\n",
     "    ```\n",
-    "    or simple\n",
+    "    or simply\n",
     "    ```python\n",
     "    value = fieldset.U[particle]\n",
     "    ```\n",
@@ -321,7 +321,14 @@
     "  \n",
     "  - It is advised _not_ to update the particle location (`particle.lon`, `particle.lat`, `particle.depth`, and/or `particle.time`) directly, as that can negatively interfere with the way that particle movements by different kernels are vectorially added. Use `particle_dlon`, `particle_dlat`, `particle_ddepth`, and/or `particle_dtime` instead. See also the [kernel loop tutorial](https://docs.oceanparcels.org/en/latest/examples/tutorial_kernelloop.html).\n",
     "  \n",
-    "  - Note that one has to be careful with writing kernels for vector fields on Curvilinear grids. While Parcels automatically rotates the U and V field when necessary, this is not the case for for example wind data. In that case, a custom rotation function will have to be written.\n"
+    "  - Note that one has to be careful with writing kernels for vector fields on Curvilinear grids. While Parcels automatically rotates the U and V field when necessary, this is not the case for for example wind data. In that case, a custom rotation function will have to be written.\n",
+    "\n",
+    "  <div class=\"alert alert-info\">\n",
+    "  A note on Field interpoaltion notation\n",
+    "\n",
+    "  Note that for the interpolation of a `Field`, the second option (`value = fieldset.U[particle]`) is not only a short-hand notation for the (`value = fieldset.U[time, particle.depth, particle.lat, particle.lon]`); it is actually a _faster_ way to interpolate the field at the particle location in Scipy mode, as described in [this section of the JIT-vs-Scipy tutorial](https://docs.oceanparcels.org/en/latest/examples/tutorial_jit_vs_scipy.html#Further-digging-into-Scipy-mode:-adding-particle-keyword-to-Field-sampling).\n",
+    "  \n",
+    "  </div>"
    ]
   },
   {

parcels/compilation/codegenerator.py

@@ -346,7 +346,7 @@ def visit_Call(self, node):
             convert = True
             if "applyConversion" in node.keywords:
                 k = node.keywords["applyConversion"]
-                if isinstance(k, ast.NameConstant):
+                if isinstance(k, ast.Constant):
                     convert = k.value
 
             # convert args to Index(Tuple(*args))
@@ -364,7 +364,7 @@ def visit_Call(self, node):
             convert = True
             if "applyConversion" in node.keywords:
                 k = node.keywords["applyConversion"]
-                if isinstance(k, ast.NameConstant):
+                if isinstance(k, ast.Constant):
                     convert = k.value
 
             # convert args to Index(Tuple(*args))
@@ -463,8 +463,7 @@ def _check_FieldSamplingArguments(ccode):
     def visit_FunctionDef(self, node):
         # Generate "ccode" attribute by traversing the Python AST
         for stmt in node.body:
-            if not (hasattr(stmt, 'value') and type(stmt.value) is ast.Str):  # ignore docstrings
-                self.visit(stmt)
+            self.visit(stmt)
 
         # Create function declaration and argument list
         decl = c.Static(c.DeclSpecifier(c.Value("StatusCode", node.name), spec='inline'))
@@ -492,7 +491,7 @@ def visit_FunctionDef(self, node):
             body += [c.Statement(f"particles->{coord}[pnum] = particles->{coord}_nextloop[pnum]")]
         body += [c.Statement("particles->time[pnum] = particles->time_nextloop[pnum]")]
 
-        body += [stmt.ccode for stmt in node.body if not (hasattr(stmt, 'value') and type(stmt.value) is ast.Str)]
+        body += [stmt.ccode for stmt in node.body]
 
         for coord in ['lon', 'lat', 'depth']:
             body += [c.Statement(f"particles->{coord}_nextloop[pnum] = particles->{coord}[pnum] + particle_d{coord}")]
@@ -878,14 +877,14 @@ def visit_Print(self, node):
         node.ccode = c.Statement(f'printf("{stat}\\n", {vars})')
 
     def visit_Constant(self, node):
-        if node.s == 'parcels_customed_Cfunc_pointer_args':
-            node.ccode = node.s
-        elif isinstance(node.s, str):
+        if node.value == 'parcels_customed_Cfunc_pointer_args':
+            node.ccode = node.value
+        elif isinstance(node.value, str):
             node.ccode = ''  # skip strings from docstrings or comments
-        elif isinstance(node.s, bool):
-            node.ccode = "1" if node.s is True else "0"
+        elif isinstance(node.value, bool):
+            node.ccode = "1" if node.value is True else "0"
         else:
-            node.ccode = str(node.n)
+            node.ccode = str(node.value)
 
 
 class LoopGenerator:

parcels/fieldset.py

@@ -1,3 +1,5 @@
+import importlib.util
+import sys
 from copy import deepcopy
 from glob import glob
 from os import path
@@ -1062,6 +1064,32 @@ def from_xarray_dataset(cls, ds, variables, dimensions, mesh='spherical', allow_
         v = fields.pop('V', None)
         return cls(u, v, fields=fields)
 
+    @classmethod
+    def from_modulefile(cls, filename, modulename="create_fieldset", **kwargs):
+        """Initialises FieldSet data from a file containing a python module file with a create_fieldset() function.
+
+        Parameters
+        ----------
+        filename: path to a python file containing at least a function which returns a FieldSet object.
+        modulename: name of the function in the python file that returns a FieldSet object. Default is "create_fieldset".
+        """
+        # check if filename exists
+        if not path.exists(filename):
+            raise IOError(f"FieldSet module file {filename} does not exist")
+
+        # Importing the source file directly (following https://docs.python.org/3/library/importlib.html#importing-a-source-file-directly)
+        spec = importlib.util.spec_from_file_location(modulename, filename)
+        fieldset_module = importlib.util.module_from_spec(spec)
+        sys.modules[modulename] = fieldset_module
+        spec.loader.exec_module(fieldset_module)
+
+        if not hasattr(fieldset_module, modulename):
+            raise IOError(f"{filename} does not contain a {modulename} function")
+        fieldset = getattr(fieldset_module, modulename)(**kwargs)
+        if not isinstance(fieldset, FieldSet):
+            raise IOError(f"Module {filename}.{modulename} does not return a FieldSet object")
+        return fieldset
+
     def get_fields(self):
         """Returns a list of all the :class:`parcels.field.Field` and :class:`parcels.field.VectorField`
         objects associated with this FieldSet.

parcels/particleset.py

@@ -878,8 +878,10 @@ def execute(self, pyfunc=AdvectionRK4, pyfunc_inter=None, endtime=None, runtime=
             runtime = runtime.total_seconds()
         if isinstance(dt, delta):
             dt = dt.total_seconds()
-        if dt > 0 and dt <= 1e-6:
+        if abs(dt) <= 1e-6:
             raise ValueError('Time step dt is too small')
+        if (dt * 1e6) % 1 != 0:
+            raise ValueError('Output interval should not have finer precision than 1e-6 s')
         outputdt = output_file.outputdt if output_file else np.infty
         if isinstance(outputdt, delta):
             outputdt = outputdt.total_seconds()

pyproject.toml

@@ -19,6 +19,21 @@ classifiers = [
   "Topic :: Scientific/Engineering",
   "Intended Audience :: Science/Research",
 ]
+dependencies = [
+  "cgen",
+  "cftime",
+  "numpy",
+  "dask",
+  "cftime",
+  "psutil",
+  "netCDF4",
+  "zarr",
+  "tqdm",
+  "pymbolic",
+  "pytest",
+  "scipy",
+  "xarray",
+]
 
 [project.urls]
 homepage = "https://oceanparcels.org/"

tests/test_data/fieldset_nemo.py

@@ -0,0 +1,18 @@
+from os import path
+
+import parcels
+
+
+def create_fieldset(indices=None):
+    data_path = path.join(path.dirname(__file__))
+
+    filenames = {'U': {'lon': path.join(data_path, 'mask_nemo_cross_180lon.nc'),
+                       'lat': path.join(data_path, 'mask_nemo_cross_180lon.nc'),
+                       'data': path.join(data_path, 'Uu_eastward_nemo_cross_180lon.nc')},
+                 'V': {'lon': path.join(data_path, 'mask_nemo_cross_180lon.nc'),
+                       'lat': path.join(data_path, 'mask_nemo_cross_180lon.nc'),
+                       'data': path.join(data_path, 'Vv_eastward_nemo_cross_180lon.nc')}}
+    variables = {'U': 'U', 'V': 'V'}
+    dimensions = {'lon': 'glamf', 'lat': 'gphif', 'time': 'time_counter'}
+    indices = indices or {}
+    return parcels.FieldSet.from_nemo(filenames, variables, dimensions, indices=indices)

tests/test_data/fieldset_nemo_error.py

@@ -0,0 +1,21 @@
+from os import path
+
+import parcels
+
+
+def random_function_name():
+    data_path = path.join(path.dirname(__file__))
+
+    filenames = {'U': {'lon': path.join(data_path, 'mask_nemo_cross_180lon.nc'),
+                       'lat': path.join(data_path, 'mask_nemo_cross_180lon.nc'),
+                       'data': path.join(data_path, 'Uu_eastward_nemo_cross_180lon.nc')},
+                 'V': {'lon': path.join(data_path, 'mask_nemo_cross_180lon.nc'),
+                       'lat': path.join(data_path, 'mask_nemo_cross_180lon.nc'),
+                       'data': path.join(data_path, 'Vv_eastward_nemo_cross_180lon.nc')}}
+    variables = {'U': 'U', 'V': 'V'}
+    dimensions = {'lon': 'glamf', 'lat': 'gphif', 'time': 'time_counter'}
+    return parcels.FieldSet.from_nemo(filenames, variables, dimensions)
+
+
+def none_returning_function():
+    return None