Merge pull request #503 from OceanParcels/new_kernel_order

New kernel order

parcels/codegenerator.py

@@ -909,10 +909,12 @@ def generate(self, funcname, field_args, const_args, kernel_ast, c_include):
         particle_backup = c.Statement("%s particle_backup" % self.ptype.name)
         sign_end_part = c.Assign("sign_end_part", "endtime - particles[p].time > 0 ? 1 : -1")
         dt_pos = c.Assign("__dt", "fmin(fabs(particles[p].dt), fabs(endtime - particles[p].time))")
+        pdt_eq_dt_pos = c.Assign("particles[p].dt", "__dt * sign_dt")
         dt_0_break = c.If("particles[p].dt == 0", c.Statement("break"))
         notstarted_continue = c.If("(sign_end_part != sign_dt) && (particles[p].dt != 0)",
                                    c.Statement("continue"))
         body = [c.Statement("set_particle_backup(&particle_backup, &(particles[p]))")]
+        body += [pdt_eq_dt_pos]
         body += [c.Assign("res", "%s(&(particles[p]), %s)" % (funcname, fargs_str))]
         body += [c.Assign("particles[p].state", "res")]  # Store return code on particle
         body += [c.If("res == SUCCESS", c.Block([c.Statement("particles[p].time += sign_dt * __dt"),

parcels/examples/example_moving_eddies.py

@@ -186,7 +186,7 @@ def test_periodic_and_computeTimeChunk_eddies(mode):
                                  lon=[3.3, 3.3],
                                  lat=[46.0, 47.8])
 
-    def periodicBC(particle, fieldset, time, dt):
+    def periodicBC(particle, fieldset, time):
         if particle.lon < fieldset.halo_west:
             particle.lon += fieldset.halo_east - fieldset.halo_west
         elif particle.lon > fieldset.halo_east:
@@ -196,9 +196,9 @@ def periodicBC(particle, fieldset, time, dt):
         elif particle.lat > fieldset.halo_north:
             particle.lat -= fieldset.halo_north - fieldset.halo_south
 
-    def slowlySouthWestward(particle, fieldset, time, dt):
-        particle.lon = particle.lon - 5 * dt / 1e5
-        particle.lat -= 3 * dt / 1e5
+    def slowlySouthWestward(particle, fieldset, time):
+        particle.lon = particle.lon - 5 * particle.dt / 1e5
+        particle.lat -= 3 * particle.dt / 1e5
 
     kernels = pset.Kernel(AdvectionRK4)+slowlySouthWestward+periodicBC
     pset.execute(kernels, runtime=delta(days=6), dt=delta(hours=1))

parcels/examples/example_nemo_curvilinear.py

@@ -27,7 +27,7 @@ def run_nemo_curvilinear(mode, outfile):
     lonp = 30 * np.ones(npart)
     latp = [i for i in np.linspace(-70, 88, npart)]
 
-    def periodicBC(particle, pieldSet, time, dt):
+    def periodicBC(particle, fieldSet, time):
         if particle.lon > 180:
             particle.lon -= 360
 

parcels/examples/example_peninsula.py

@@ -67,8 +67,8 @@ def peninsula_fieldset(xdim, ydim, mesh='flat'):
     return FieldSet.from_data(data, dimensions, mesh=mesh)
 
 
-def UpdateP(particle, fieldset, time, dt):
-    particle.p = fieldset.P[time, particle.lon, particle.lat, particle.depth]
+def UpdateP(particle, fieldset, time):
+    particle.p = fieldset.P[time, particle.depth, particle.lat, particle.lon]
 
 
 def pensinsula_example(fieldset, npart, mode='jit', degree=1,
@@ -124,7 +124,7 @@ def test_peninsula_fieldset(mode, mesh):
     err_adv = np.array([abs(p.p_start - p.p) for p in pset])
     assert(err_adv <= 1.e-3).all()
     # Test Field sampling accuracy by comparing kernel against Field sampling
-    err_smpl = np.array([abs(p.p - pset.fieldset.P[0., p.lon, p.lat, p.depth]) for p in pset])
+    err_smpl = np.array([abs(p.p - pset.fieldset.P[0., p.depth, p.lat, p.lon]) for p in pset])
     assert(err_smpl <= 1.e-3).all()
 
 
@@ -147,7 +147,7 @@ def test_peninsula_file(mode, mesh):
     err_adv = np.array([abs(p.p_start - p.p) for p in pset])
     assert(err_adv <= 1.e-3).all()
     # Test Field sampling accuracy by comparing kernel against Field sampling
-    err_smpl = np.array([abs(p.p - pset.fieldset.P[0., p.lon, p.lat, p.depth]) for p in pset])
+    err_smpl = np.array([abs(p.p - pset.fieldset.P[0., p.depth, p.lat, p.lon]) for p in pset])
     assert(err_smpl <= 1.e-3).all()
 
 

parcels/examples/example_recursive_errorhandling.py

@@ -35,13 +35,13 @@ def test_recursive_errorhandling(mode, xdim=2, ydim=2):
     pset = ParticleSet.from_line(fieldset=fieldset, pclass=ptype[mode],
                                  start=(0.5, 0.5), finish=(0.5, 0.5), size=10)
 
-    def TestLon(particle, fieldset, time, dt):
+    def TestLon(particle, fieldset, time):
         """Kernel to check whether a longitude is larger than fieldset.minlon.
         If not, the Kernel throws an error"""
         if particle.lon <= fieldset.minlon:
             return ErrorCode.Error
 
-    def Error_RandomiseLon(particle, fieldset, time, dt):
+    def Error_RandomiseLon(particle, fieldset, time):
         """Error handling kernel that draws a new longitude.
         Note that this new longitude can be smaller than fieldset.minlon"""
         particle.lon = random.uniform(0., 1.)

parcels/examples/example_stommel.py

@@ -52,8 +52,8 @@ def stommel_fieldset(xdim=200, ydim=200):
     return FieldSet.from_data(data, dimensions, mesh='flat', transpose=True)
 
 
-def UpdateP(particle, fieldset, time, dt):
-    particle.p = fieldset.P[time, particle.lon, particle.lat, particle.depth]
+def UpdateP(particle, fieldset, time):
+    particle.p = fieldset.P[time, particle.depth, particle.lat, particle.lon]
 
 
 def stommel_example(npart=1, mode='jit', verbose=False, method=AdvectionRK4):

parcels/examples/parcels_tutorial.ipynb

@@ -1295,10 +1295,10 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "def WestVel(particle, fieldset, time, dt):\n",
+    "def WestVel(particle, fieldset, time):\n",
     "    if time > 86400:\n",
     "        uvel = -2.\n",
-    "        particle.lon += uvel * dt"
+    "        particle.lon += uvel * particle.dt"
    ]
   },
   {
@@ -1632,8 +1632,8 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "def SampleP(particle, fieldset, time, dt):  # Custom function that samples fieldset.P at particle location\n",
-    "    particle.p = fieldset.P[time, particle.lon, particle.lat, particle.depth]\n",
+    "def SampleP(particle, fieldset, time):  # Custom function that samples fieldset.P at particle location\n",
+    "    particle.p = fieldset.P[time, particle.depth, particle.lat, particle.lon]\n",
     "\n",
     "k_sample = pset.Kernel(SampleP)    # Casting the SampleP function to a kernel."
    ]
@@ -1740,7 +1740,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "def TotalDistance(particle, fieldset, time, dt):\n",
+    "def TotalDistance(particle, fieldset, time):\n",
     "    # Calculate the distance in latitudinal direction (using 1.11e2 kilometer per degree latitude)\n",
     "    lat_dist = (particle.lat - particle.prev_lat) * 1.11e2\n",
     "    # Calculate the distance in longitudinal direction, using cosine(latitude) - spherical earth\n",

parcels/examples/tutorial_Agulhasparticles.ipynb

@@ -94,7 +94,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "def DeleteParticle(particle, fieldset, time, dt):\n",
+    "def DeleteParticle(particle, fieldset, time):\n",
     "    particle.delete()"
    ]
   },

parcels/examples/tutorial_Argofloats.ipynb

@@ -21,7 +21,7 @@
    "outputs": [],
    "source": [
     "# Define the new Kernel that mimics Argo vertical movement\n",
-    "def ArgoVerticalMovement(particle, fieldset, time, dt):\n",
+    "def ArgoVerticalMovement(particle, fieldset, time):\n",
     "    driftdepth = 1000  # maximum depth in m\n",
     "    maxdepth = 2000  # maximum depth in m\n",
     "    vertical_speed = 0.10  # sink and rise speed in m/s\n",
@@ -30,26 +30,26 @@
     "\n",
     "    if particle.cycle_phase == 0:\n",
     "        # Phase 0: Sinking with vertical_speed until depth is driftdepth\n",
-    "        particle.depth += vertical_speed * dt\n",
+    "        particle.depth += vertical_speed * particle.dt\n",
     "        if particle.depth >= driftdepth:\n",
     "            particle.cycle_phase = 1\n",
     "\n",
     "    elif particle.cycle_phase == 1:\n",
     "        # Phase 1: Drifting at depth for drifttime seconds\n",
-    "        particle.drift_age += dt\n",
+    "        particle.drift_age += particle.dt\n",
     "        if particle.drift_age >= drifttime:\n",
     "            particle.drift_age = 0  # reset drift_age for next cycle\n",
     "            particle.cycle_phase = 2\n",
     "\n",
     "    elif particle.cycle_phase == 2:\n",
     "        # Phase 2: Sinking further to maxdepth\n",
-    "        particle.depth += vertical_speed * dt\n",
+    "        particle.depth += vertical_speed * particle.dt\n",
     "        if particle.depth >= maxdepth:\n",
     "            particle.cycle_phase = 3\n",
     "\n",
     "    elif particle.cycle_phase == 3:\n",
     "        # Phase 3: Rising with vertical_speed until at surface\n",
-    "        particle.depth -= vertical_speed * dt\n",
+    "        particle.depth -= vertical_speed * particle.dt\n",
     "        #particle.temp = fieldset.temp[time, particle.lon, particle.lat, particle.depth]  # if fieldset has temperature\n",
     "        if particle.depth <= fieldset.mindepth:\n",
     "            particle.depth = fieldset.mindepth\n",
@@ -62,7 +62,7 @@
     "            particle.cycle_phase = 0\n",
     "            particle.cycle_age = 0\n",
     "\n",
-    "    particle.cycle_age += dt  # update cycle_age"
+    "    particle.cycle_age += particle.dt  # update cycle_age"
    ]
   },
   {

parcels/examples/tutorial_NestedFields.ipynb

@@ -208,8 +208,8 @@
    "source": [
     "from parcels import Variable\n",
     "\n",
-    "def SampleNestedFieldIndex(particle, fieldset, time, dt):\n",
-    "    particle.f = fieldset.F[time, particle.lon, particle.lat, particle.depth]\n",
+    "def SampleNestedFieldIndex(particle, fieldset, time):\n",
+    "    particle.f = fieldset.F[time, particle.depth, particle.lat, particle.lon]\n",
     "\n",
     "class SampleParticle(JITParticle):\n",
     "    f = Variable('f', dtype=np.int32)\n",

parcels/examples/tutorial_nemo_curvilinear.ipynb

@@ -132,7 +132,7 @@
     "latp = [i for i in np.linspace(-70, 88, npart)]\n",
     "\n",
     "# Create a periodic boundary condition kernel\n",
-    "def periodicBC(particle, fieldset, time, dt):\n",
+    "def periodicBC(particle, fieldset, time):\n",
     "    if particle.lon > 180:\n",
     "        particle.lon -= 360\n",
     "\n",

parcels/examples/tutorial_periodic_boundaries.ipynb

@@ -95,7 +95,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "def periodicBC(particle, fieldset, time, dt):\n",
+    "def periodicBC(particle, fieldset, time):\n",
     "    if particle.lon < fieldset.halo_west:\n",
     "        particle.lon += fieldset.halo_east - fieldset.halo_west\n",
     "    elif particle.lon > fieldset.halo_east:\n",

parcels/field.py

@@ -175,6 +175,8 @@ def from_netcdf(cls, filenames, variable, dimensions, indices=None, grid=None,
                     raise NotImplementedError('longitude and latitude dimensions are currently processed together from one single file')
                 lonlat_filename = filenames['lon'][0]
                 if 'depth' in dimensions:
+                    if not isinstance(filenames['depth'], Iterable) or isinstance(filenames['depth'], str):
+                        filenames['depth'] = [filenames['depth']]
                     if len(filenames['depth']) != 1:
                         raise NotImplementedError('Vertically adaptive meshes not implemented for from_netcdf()')
                     depth_filename = filenames['depth'][0]
@@ -186,11 +188,11 @@ def from_netcdf(cls, filenames, variable, dimensions, indices=None, grid=None,
         indices = {} if indices is None else indices.copy()
         for ind in indices.values():
             assert np.min(ind) >= 0, \
-                ('Negative indices are currently not allowed in Parcels. ' +
-                 'This is related to the non-increasing dimension it could generate ' +
-                 'if the domain goes from lon[-4] to lon[6] for example. ' +
-                 'Please raise an issue on https://github.com/OceanParcels/parcels/issues ' +
-                 'if you would need such feature implemented.')
+                ('Negative indices are currently not allowed in Parcels. '
+                 + 'This is related to the non-increasing dimension it could generate '
+                 + 'if the domain goes from lon[-4] to lon[6] for example. '
+                 + 'Please raise an issue on https://github.com/OceanParcels/parcels/issues '
+                 + 'if you would need such feature implemented.')
 
         with NetcdfFileBuffer(lonlat_filename, dimensions, indices, netcdf_engine) as filebuffer:
             lon, lat = filebuffer.read_lonlat
@@ -726,7 +728,7 @@ def depth_index(self, depth, lat, lon):
         else:
             return depth_index.argmin() - 1 if depth_index.any() else 0
 
-    def eval(self, time, x, y, z, applyConversion=True):
+    def eval(self, time, z, y, x, applyConversion=True):
         """Interpolate field values in space and time.
 
         We interpolate linearly in time and apply implicit unit
@@ -752,12 +754,12 @@ def eval(self, time, x, y, z, applyConversion=True):
         else:
             return value
 
-    def ccode_eval(self, var, t, x, y, z):
+    def ccode_eval(self, var, t, z, y, x):
         # Casting interp_methd to int as easier to pass on in C-code
         return "temporal_interpolation(%s, %s, %s, %s, %s, particle->cxi, particle->cyi, particle->czi, particle->cti, &%s, %s)" \
             % (x, y, z, t, self.name, var, self.interp_method.upper())
 
-    def ccode_convert(self, _, x, y, z):
+    def ccode_convert(self, _, z, y, x):
         return self.units.ccode_to_target(x, y, z)
 
     @property
@@ -924,10 +926,10 @@ class SummedField(list):
     Also note that, since SummedFields are lists, the individual Fields can
     still be queried through their list index (e.g. SummedField[1]).
     """
-    def eval(self, time, x, y, z, applyConversion=True):
+    def eval(self, time, z, y, x, applyConversion=True):
         tmp = 0
         for fld in self:
-            tmp += fld.eval(time, x, y, z, applyConversion)
+            tmp += fld.eval(time, z, y, x, applyConversion)
         return tmp
 
     def __getitem__(self, key):
@@ -1151,18 +1153,18 @@ def spatial_c_grid_interpolation3D(self, ti, z, y, x, time):
             (u, v, w) = self.spatial_c_grid_interpolation3D_full(ti, z, y, x, time)
         else:
             (u, v) = self.spatial_c_grid_interpolation2D(ti, z, y, x, time)
-            w = self.W.eval(time, x, y, z, False)
+            w = self.W.eval(time, z, y, x, False)
             w = self.W.units.to_target(w, x, y, z)
         return (u, v, w)
 
-    def eval(self, time, x, y, z):
+    def eval(self, time, z, y, x):
         if self.U.interp_method != 'cgrid_velocity':
-            u = self.U.eval(time, x, y, z, False)
-            v = self.V.eval(time, x, y, z, False)
+            u = self.U.eval(time, z, y, x, False)
+            v = self.V.eval(time, z, y, x, False)
             u = self.U.units.to_target(u, x, y, z)
             v = self.V.units.to_target(v, x, y, z)
             if self.W is not None:
-                w = self.W.eval(time, x, y, z, False)
+                w = self.W.eval(time, z, y, x, False)
                 w = self.W.units.to_target(w, x, y, z)
                 return (u, v, w)
             else:
@@ -1199,7 +1201,7 @@ def eval(self, time, x, y, z):
     def __getitem__(self, key):
         return self.eval(*key)
 
-    def ccode_eval(self, varU, varV, varW, U, V, W, t, x, y, z):
+    def ccode_eval(self, varU, varV, varW, U, V, W, t, z, y, x):
         # Casting interp_methd to int as easier to pass on in C-code
         if varW:
             return "temporal_interpolationUVW(%s, %s, %s, %s, %s, %s, %s, " \
@@ -1235,13 +1237,13 @@ def __init__(self, name, U, V, W=None):
         self.V = V
         self.W = W
 
-    def eval(self, time, x, y, z):
+    def eval(self, time, z, y, x):
         zonal = meridional = vertical = 0
         if self.W is not None:
             for (U, V, W) in zip(self.U, self.V, self.W):
                 vfld = VectorField(self.name, U, V, W)
                 vfld.fieldset = self.fieldset
-                (tmp1, tmp2, tmp3) = vfld.eval(time, x, y, z)
+                (tmp1, tmp2, tmp3) = vfld.eval(time, z, y, x)
                 zonal += tmp1
                 meridional += tmp2
                 vertical += tmp3
@@ -1250,7 +1252,7 @@ def eval(self, time, x, y, z):
             for (U, V) in zip(self.U, self.V):
                 vfld = VectorField(self.name, U, V)
                 vfld.fieldset = self.fieldset
-                (tmp1, tmp2) = vfld.eval(time, x, y, z)
+                (tmp1, tmp2) = vfld.eval(time, z, y, x)
                 zonal += tmp1
                 meridional += tmp2
             return (zonal, meridional)

parcels/fieldset.py

@@ -434,18 +434,6 @@ def add_periodic_halo(self, zonal=False, meridional=False, halosize=5):
             if isinstance(value, Field):
                 value.add_periodic_halo(zonal, meridional, halosize)
 
-    def eval(self, x, y):
-        """Evaluate the zonal and meridional velocities (u,v) at a point (x,y)
-
-        :param x: zonal point to evaluate
-        :param y: meridional point to evaluate
-
-        :return u, v: zonal and meridional velocities at point"""
-
-        u = self.U.eval(x, y)
-        v = self.V.eval(x, y)
-        return u, v
-
     def write(self, filename):
         """Write FieldSet to NetCDF file using NEMO convention
 

parcels/kernel.py

@@ -85,6 +85,9 @@ def __init__(self, fieldset, ptype, pyfunc=None, funcname=None,
             self.pyfunc = user_ctx[self.funcname]
         else:
             self.pyfunc = pyfunc
+        assert len(inspect.getargspec(self.pyfunc).args) == 3, \
+            'Since Parcels v2.0, kernels do only take 3 arguments: particle, fieldset, time !! AND !! Argument order in field interpolation is time, depth, lat, lon.'
+
         self.name = "%s%s" % (ptype.name, self.funcname)
 
         # Generate the kernel function and add the outer loop
@@ -208,7 +211,8 @@ def execute_python(self, pset, endtime, dt):
                 for var in ptype.variables:
                     p_var_back[var.name] = getattr(p, var.name)
                 try:
-                    res = self.pyfunc(p, pset.fieldset, p.time, sign_dt * dt_pos)
+                    p.dt = sign_dt * dt_pos
+                    res = self.pyfunc(p, pset.fieldset, p.time)
                 except FieldSamplingError as fse:
                     res = ErrorCode.ErrorOutOfBounds
                     p.exception = fse
@@ -274,7 +278,7 @@ def remove_deleted(pset):
                 else:
                     recovery_kernel = recovery_map[p.state]
                     p.state = ErrorCode.Success
-                    recovery_kernel(p, self.fieldset, p.time, dt)
+                    recovery_kernel(p, self.fieldset, p.time)
 
             # Remove all particles that signalled deletion
             remove_deleted(pset)