Merge pull request #435 from OceanParcels/changing_FieldList_to_SummedFields

Changing FieldLists to SummedFields

Author: erikvansebille

parcels/codegenerator.py

@@ -1,5 +1,4 @@
-from parcels.field import Field, VectorField
-from parcels.fieldset import FieldList, VectorFieldList
+from parcels.field import Field, VectorField, SummedField, SummedVectorField
 from parcels.loggers import logger
 import ast
 import cgen as c
@@ -20,12 +19,12 @@ def __getattr__(self, attr):
         if isinstance(getattr(self.obj, attr), Field):
             return FieldNode(getattr(self.obj, attr),
                              ccode="%s->%s" % (self.ccode, attr))
-        elif isinstance(getattr(self.obj, attr), VectorFieldList):
-            return VectorFieldListNode(getattr(self.obj, attr),
-                                       ccode="%s->%s" % (self.ccode, attr))
-        elif isinstance(getattr(self.obj, attr), FieldList) or isinstance(getattr(self.obj, attr), list):
-            return FieldListNode(getattr(self.obj, attr),
-                                 ccode="%s->%s" % (self.ccode, attr))
+        elif isinstance(getattr(self.obj, attr), SummedVectorField):
+            return SummedVectorFieldNode(getattr(self.obj, attr),
+                                         ccode="%s->%s" % (self.ccode, attr))
+        elif isinstance(getattr(self.obj, attr), SummedField) or isinstance(getattr(self.obj, attr), list):
+            return SummedFieldNode(getattr(self.obj, attr),
+                                   ccode="%s->%s" % (self.ccode, attr))
         elif isinstance(getattr(self.obj, attr), VectorField):
             return VectorFieldNode(getattr(self.obj, attr),
                                    ccode="%s->%s" % (self.ccode, attr))
@@ -60,24 +59,24 @@ def __init__(self, field, args, var, var2, var3):
         self.var3 = var3  # third variable for UVW interpolation
 
 
-class FieldListNode(IntrinsicNode):
+class SummedFieldNode(IntrinsicNode):
     def __getitem__(self, attr):
-        return FieldListEvalNode(self.obj, attr)
+        return SummedFieldEvalNode(self.obj, attr)
 
 
-class FieldListEvalNode(IntrinsicNode):
+class SummedFieldEvalNode(IntrinsicNode):
     def __init__(self, fields, args, var):
         self.fields = fields
         self.args = args
         self.var = var  # the variable in which the interpolated field is written
 
 
-class VectorFieldListNode(IntrinsicNode):
+class SummedVectorFieldNode(IntrinsicNode):
     def __getitem__(self, attr):
-        return VectorFieldListEvalNode(self.obj, attr)
+        return SummedVectorFieldEvalNode(self.obj, attr)
 
 
-class VectorFieldListEvalNode(IntrinsicNode):
+class SummedVectorFieldEvalNode(IntrinsicNode):
     def __init__(self, field, args, var, var2, var3):
         self.field = field
         self.args = args
@@ -210,18 +209,18 @@ def visit_Subscript(self, node):
 
         # If we encounter field evaluation we replace it with a
         # temporary variable and put the evaluation call on the stack.
-        if isinstance(node.value, FieldListNode):
+        if isinstance(node.value, SummedFieldNode):
             tmp = [self.get_tmp() for _ in node.value.obj]
             # Insert placeholder node for field eval ...
-            self.stmt_stack += [FieldListEvalNode(node.value, node.slice, tmp)]
+            self.stmt_stack += [SummedFieldEvalNode(node.value, node.slice, tmp)]
             # .. and return the name of the temporary that will be populated
             return ast.Name(id='+'.join(tmp))
-        elif isinstance(node.value, VectorFieldListNode):
+        elif isinstance(node.value, SummedVectorFieldNode):
             tmp = [self.get_tmp() for _ in node.value.obj.U]
             tmp2 = [self.get_tmp() for _ in node.value.obj.U]
             tmp3 = [self.get_tmp() if node.value.obj.W else None for _ in node.value.obj.U]
             # Insert placeholder node for field eval ...
-            self.stmt_stack += [VectorFieldListEvalNode(node.value, node.slice, tmp, tmp2, tmp3)]
+            self.stmt_stack += [SummedVectorFieldEvalNode(node.value, node.slice, tmp, tmp2, tmp3)]
             # .. and return the name of the temporary that will be populated
             if all(tmp3):
                 return ast.Tuple([ast.Name(id='+'.join(tmp)), ast.Name(id='+'.join(tmp2)), ast.Name(id='+'.join(tmp3))], ast.Load())
@@ -621,7 +620,7 @@ def visit_FieldNode(self, node):
         """Record intrinsic fields used in kernel"""
         self.field_args[node.obj.name] = node.obj
 
-    def visit_FieldListNode(self, node):
+    def visit_SummedFieldNode(self, node):
         """Record intrinsic fields used in kernel"""
         for fld in node.obj:
             self.field_args[fld.name] = fld
@@ -630,7 +629,7 @@ def visit_VectorFieldNode(self, node):
         """Record intrinsic fields used in kernel"""
         self.vector_field_args[node.obj.name] = node.obj
 
-    def visit_VectorFieldListNode(self, node):
+    def visit_SummedVectorFieldNode(self, node):
         """Record intrinsic fields used in kernel"""
         for fld in node.obj.U:
             self.field_args[fld.name] = fld
@@ -672,7 +671,7 @@ def visit_VectorFieldEvalNode(self, node):
         node.ccode = c.Block([c.Assign("err", ccode_eval),
                               conv_stat, c.Statement("CHECKERROR(err)")])
 
-    def visit_FieldListEvalNode(self, node):
+    def visit_SummedFieldEvalNode(self, node):
         self.visit(node.fields)
         self.visit(node.args)
         cstat = []
@@ -683,7 +682,7 @@ def visit_FieldListEvalNode(self, node):
             cstat += [c.Assign("err", ccode_eval), conv_stat, c.Statement("CHECKERROR(err)")]
         node.ccode = c.Block(cstat)
 
-    def visit_VectorFieldListEvalNode(self, node):
+    def visit_SummedVectorFieldEvalNode(self, node):
         self.visit(node.field)
         self.visit(node.args)
         cstat = []

…rcels/examples/tutorial_FieldLists.ipynb …els/examples/tutorial_SummedFields.ipynbparcels/examples/tutorial_FieldLists.ipynb renamed to parcels/examples/tutorial_SummedFields.ipynb parcels/examples/tutorial_FieldLists.ipynb renamed to parcels/examples/tutorial_SummedFields.ipynb

@@ -4,7 +4,7 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "## Tutorial on how to combine different Fields for advection into a `FieldList` object"
+    "## Tutorial on how to combine different Fields for advection into a `SummedField` object"
    ]
   },
   {
@@ -13,9 +13,9 @@
    "source": [
     "In some oceanographic applications, you may want to advect particles using a combination of different velocity data sets. For example, particles at the surface are transported by a combination of geostrophic, Ekman and Stokes flow. And often, these flows are not even on the same grid.\n",
     "\n",
-    "One option would be to write a  `Kernel` that computes the movement of particles due to each of these flows. However, in Parcels it is possible to directly combine different flows (without interpolation) and feed them into the built-in `AdvectionRK4` kernel. For that, we use so-called `FieldList` objects.\n",
+    "One option would be to write a  `Kernel` that computes the movement of particles due to each of these flows. However, in Parcels it is possible to directly combine different flows (without interpolation) and feed them into the built-in `AdvectionRK4` kernel. For that, we use so-called `SummedField` objects.\n",
     "\n",
-    "This tutorial shows how to use these `FieldLists` with a very idealised example. We start by importing the relevant modules."
+    "This tutorial shows how to use these `SummedField` with a very idealised example. We start by importing the relevant modules."
    ]
   },
   {
@@ -123,7 +123,7 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "Now comes the trick of the `FieldLists`. We can simply define a new `FieldSet` with a list of different `Fields`, as in `U=[U, Ustokes]`."
+    "Now comes the trick of the `SummedFields`. We can simply define a new `FieldSet` with a summation of different `Fields`, as in `U=U+Ustokes`."
    ]
   },
   {
@@ -132,7 +132,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "fieldset = FieldSet(U=[U, Ustokes], V=[V, Vstokes])"
+    "fieldset = FieldSet(U=U+Ustokes, V=V+Vstokes)"
    ]
   },
   {
@@ -177,9 +177,9 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "What happens under the hood is that each `Field` in the `FieldList` is interpolated separately to the particle location, and that the different velocities are added in each step of the RK4 advection. So `FieldLists` are effortless to users\n",
+    "What happens under the hood is that each `Field` in the `SummedField` is interpolated separately to the particle location, and that the different velocities are added in each step of the RK4 advection. So `SummedFields` are effortless to users\n",
     "\n",
-    "Note that `FieldLists` work for any type of `Field`, not only for velocities. Any call to a `Field` interpolation (`fieldset.fld[time, lon, lat, depth]`) will return the sum of all `Fields` in the list if `fld` is a `FieldList`."
+    "Note that `SummedFields` work for any type of `Field`, not only for velocities. Any call to a `Field` interpolation (`fieldset.fld[time, lon, lat, depth]`) will return the sum of all `Fields` in the list if `fld` is a `SummedField`."
    ]
   },
   {
@@ -206,7 +206,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython2",
-   "version": "2.7.15"
+   "version": "2.7.14"
   }
  },
  "nbformat": 4,

parcels/field.py

@@ -11,7 +11,8 @@
                    CurvilinearSGrid, CGrid, GridCode)
 
 
-__all__ = ['Field', 'VectorField', 'Geographic', 'GeographicPolar', 'GeographicSquare', 'GeographicPolarSquare']
+__all__ = ['Field', 'VectorField', 'SummedField', 'SummedVectorField',
+           'Geographic', 'GeographicPolar', 'GeographicSquare', 'GeographicPolarSquare']
 
 
 class FieldSamplingError(RuntimeError):
@@ -164,7 +165,7 @@ class Field(object):
            2. flat: No conversion, lat/lon are assumed to be in m.
     :param grid: :class:`parcels.grid.Grid` object containing all the lon, lat depth, time
            mesh and time_origin information. Can be constructed from any of the Grid objects
-    :param fieldtype: Type of Field to be used for UnitConverter when using FieldLists
+    :param fieldtype: Type of Field to be used for UnitConverter when using SummedFields
            (either 'U', 'V', 'Kh_zonal', 'Kh_Meridional' or None)
     :param transpose: Transpose data to required (lon, lat) layout
     :param vmin: Minimum allowed value on the field. Data below this value are set to zero
@@ -976,6 +977,42 @@ def computeTimeChunk(self, data, tindex):
 
         return data
 
+    def __add__(self, field):
+        if isinstance(self, Field) and isinstance(field, Field):
+            return SummedField([self, field])
+        elif isinstance(field, SummedField):
+            field.insert(0, self)
+            return field
+
+
+class SummedField(list):
+    """Class SummedField is a list of Fields over which Field interpolation
+    is summed. This can e.g. be used when combining multiple flow fields,
+    where the total flow is the sum of all the individual flows.
+    Note that the individual Fields can be on different Grids.
+    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):
+        tmp = 0
+        for fld in self:
+            tmp += fld.eval(time, x, y, z, applyConversion)
+        return tmp
+
+    def __getitem__(self, key):
+        if isinstance(key, int):
+            return list.__getitem__(self, key)
+        else:
+            return self.eval(*key)
+
+    def __add__(self, field):
+        if isinstance(field, Field):
+            self.append(field)
+        elif isinstance(field, SummedField):
+            for fld in field:
+                self.append(fld)
+        return self
+
 
 class VectorField(object):
     """Class VectorField stores 2 or 3 fields which defines together a vector field.
@@ -1150,6 +1187,49 @@ def ccode_eval(self, varU, varV, varW, U, V, W, t, x, y, z):
                    % (varU, varV, U.interp_method.upper())
 
 
+class SummedVectorField(list):
+    """Class SummedVectorField stores 2 or 3 SummedFields which defines together a vector field.
+    This enables to interpolate them as one single vector SummedField in the kernels.
+
+    :param name: Name of the vector field
+    :param U: SummedField defining the zonal component
+    :param V: SummedField defining the meridional component
+    :param W: SummedField defining the vertical component (default: None)
+    """
+
+    def __init__(self, name, U, V, W=None):
+        self.name = name
+        self.U = U
+        self.V = V
+        self.W = W
+
+    def eval(self, time, x, y, z):
+        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)
+                zonal += tmp1
+                meridional += tmp2
+                vertical += tmp3
+            return (zonal, meridional, vertical)
+        else:
+            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)
+                zonal += tmp1
+                meridional += tmp2
+            return (zonal, meridional)
+
+    def __getitem__(self, key):
+        if isinstance(key, int):
+            return list.__getitem__(self, key)
+        else:
+            return self.eval(*key)
+
+
 class NetcdfFileBuffer(object):
     """ Class that encapsulates and manages deferred access to file data. """