Repository: OceanParcels/parcels (python)

Issue:
FieldSet.add_field() should throw error if Field name already exist
Currently, `FieldSet.add_field()` doen's't throw an error when a user tries to add a field with a name that already exists. This can be very confusing, so Parcels should check for this

Pull request description:
Adding check for duplicate Field name in FieldSet.add_field()
This fixes #562

Patch:
diff --git a/parcels/codegenerator.py b/parcels/codegenerator.py
index 0082fc92..5dbb21fc 100644
--- a/parcels/codegenerator.py
+++ b/parcels/codegenerator.py
@@ -1,4 +1,4 @@
-from parcels.field import Field, VectorField, SummedField, SummedVectorField, NestedField
+from parcels.field import Field, VectorField, SummedField, NestedField
 from parcels.tools.loggers import logger
 import ast
 import cgen as c
@@ -27,12 +27,13 @@ class FieldSetNode(IntrinsicNode):
             else:
                 return NestedFieldNode(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))
+            if isinstance(getattr(self.obj, attr)[0], VectorField):
+                return SummedVectorFieldNode(getattr(self.obj, attr),
+                                             ccode="%s->%s" % (self.ccode, attr))
+            else:
+                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))
@@ -85,8 +86,8 @@ class SummedVectorFieldNode(IntrinsicNode):
 
 
 class SummedVectorFieldEvalNode(IntrinsicNode):
-    def __init__(self, field, args, var, var2, var3):
-        self.field = field
+    def __init__(self, fields, args, var, var2, var3):
+        self.fields = fields
         self.args = args
         self.var = var  # the variable in which the interpolated field is written
         self.var2 = var2  # second variable for UV interpolation
@@ -257,9 +258,9 @@ class IntrinsicTransformer(ast.NodeTransformer):
             # .. and return the name of the temporary that will be populated
             return ast.Name(id='+'.join(tmp))
         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]
+            tmp = [self.get_tmp() for _ in range(len(node.value.obj))]
+            tmp2 = [self.get_tmp() for _ in range(len(node.value.obj))]
+            tmp3 = [self.get_tmp() if list.__getitem__(node.value.obj, 0).vector_type == '3D' else None for _ in range(len(node.value.obj))]
             # Insert placeholder node for field eval ...
             self.stmt_stack += [SummedVectorFieldEvalNode(node.value, node.slice, tmp, tmp2, tmp3)]
             # .. and return the name of the temporary that will be populated
@@ -276,7 +277,7 @@ class IntrinsicTransformer(ast.NodeTransformer):
         elif isinstance(node.value, VectorFieldNode):
             tmp = self.get_tmp()
             tmp2 = self.get_tmp()
-            tmp3 = self.get_tmp() if node.value.obj.W else None
+            tmp3 = self.get_tmp() if node.value.obj.vector_type == '3D' else None
             # Insert placeholder node for field eval ...
             self.stmt_stack += [VectorFieldEvalNode(node.value, node.slice, tmp, tmp2, tmp3)]
             # .. and return the name of the temporary that will be populated
@@ -291,7 +292,7 @@ class IntrinsicTransformer(ast.NodeTransformer):
         elif isinstance(node.value, NestedVectorFieldNode):
             tmp = self.get_tmp()
             tmp2 = self.get_tmp()
-            tmp3 = self.get_tmp() if list.__getitem__(node.value.obj, 0).W else None
+            tmp3 = self.get_tmp() if list.__getitem__(node.value.obj, 0).vector_type == '3D' else None
             self.stmt_stack += [NestedVectorFieldEvalNode(node.value, node.slice, tmp, tmp2, tmp3)]
             if tmp3:
                 return ast.Tuple([ast.Name(id=tmp), ast.Name(id=tmp2), ast.Name(id=tmp3)], ast.Load())
@@ -412,22 +413,13 @@ class KernelGenerator(ast.NodeVisitor):
         decl = c.Static(c.DeclSpecifier(c.Value("ErrorCode", node.name), spec='inline'))
         args = [c.Pointer(c.Value(self.ptype.name, "particle")),
                 c.Value("double", "time")]
-        for field_name, field in self.field_args.items():
-            args += [c.Pointer(c.Value("CField", "%s" % field_name))]
-        for field_name, field in self.vector_field_args.items():
-            fieldset = field.fieldset
+        for field in self.field_args.values():
+            args += [c.Pointer(c.Value("CField", "%s" % field.ccode_name))]
+        for field in self.vector_field_args.values():
             Wname = field.W.ccode_name if field.W else 'not_defined'
-            for f in [field.U.ccode_name, field.V.ccode_name, Wname]:
-                try:
-                    # Next line will break for example if field.U was created but not added to the fieldset
-                    getattr(fieldset, f)
-                    if f not in self.field_args:
-                        args += [c.Pointer(c.Value("CField", "%s" % f))]
-                except:
-                    if f != Wname:
-                        raise RuntimeError("Field %s needed by a VectorField but it does not exist" % f)
-                    else:
-                        pass
+            for fname in [field.U.ccode_name, field.V.ccode_name, Wname]:
+                if fname not in self.field_args and fname != 'not_defined':
+                    args += [c.Pointer(c.Value("CField", "%s" % fname))]
         for const, _ in self.const_args.items():
             args += [c.Value("float", const)]
 
@@ -476,7 +468,7 @@ class KernelGenerator(ast.NodeVisitor):
                     pointer_args = True
                     continue
                 if isinstance(a, FieldNode) or isinstance(a, VectorFieldNode):
-                    a.ccode = a.obj.name
+                    a.ccode = a.obj.ccode_name
                 elif isinstance(a, ParticleNode):
                     continue
                 elif pointer_args:
@@ -696,13 +688,8 @@ class KernelGenerator(ast.NodeVisitor):
 
     def visit_SummedVectorFieldNode(self, node):
         """Record intrinsic fields used in kernel"""
-        for fld in node.obj.U:
-            self.field_args[fld.ccode_name] = fld
-        for fld in node.obj.V:
-            self.field_args[fld.ccode_name] = fld
-        if hasattr(node.obj, 'W') and node.obj.W:
-            for fld in node.obj.W:
-                self.field_args[fld.ccode_name] = fld
+        for fld in node.obj:
+            self.vector_field_args[fld.ccode_name] = fld
 
     def visit_NestedVectorFieldNode(self, node):
         """Record intrinsic fields used in kernel"""
@@ -734,7 +721,7 @@ class KernelGenerator(ast.NodeVisitor):
                           c.Statement("%s *= %s" % (node.var2, ccode_conv2))]
         else:
             statements = []
-        if node.var3:
+        if node.field.obj.vector_type == '3D':
             ccode_conv3 = node.field.obj.W.ccode_convert(*node.args.ccode)
             statements.append(c.Statement("%s *= %s" % (node.var3, ccode_conv3)))
         conv_stat = c.Block(statements)
@@ -753,28 +740,25 @@ class KernelGenerator(ast.NodeVisitor):
         node.ccode = c.Block(cstat)
 
     def visit_SummedVectorFieldEvalNode(self, node):
-        self.visit(node.field)
+        self.visit(node.fields)
         self.visit(node.args)
         cstat = []
-        if node.field.obj.W:
-            Wlist = node.field.obj.W
-        else:
-            Wlist = [None] * len(node.field.obj.U)
-        for U, V, W, var, var2, var3 in zip(node.field.obj.U, node.field.obj.V, Wlist, node.var, node.var2, node.var3):
-            vfld = VectorField(node.field.obj.name, U, V, W)
-            ccode_eval = vfld.ccode_eval(var, var2, var3, U, V, W, *node.args.ccode)
-            if U.interp_method != 'cgrid_velocity':
-                ccode_conv1 = U.ccode_convert(*node.args.ccode)
-                ccode_conv2 = V.ccode_convert(*node.args.ccode)
+        for fld, var, var2, var3 in zip(node.fields.obj, node.var, node.var2, node.var3):
+            ccode_eval = fld.ccode_eval(var, var2, var3,
+                                        fld.U, fld.V, fld.W,
+                                        *node.args.ccode)
+            if fld.U.interp_method != 'cgrid_velocity':
+                ccode_conv1 = fld.U.ccode_convert(*node.args.ccode)
+                ccode_conv2 = fld.V.ccode_convert(*node.args.ccode)
                 statements = [c.Statement("%s *= %s" % (var, ccode_conv1)),
                               c.Statement("%s *= %s" % (var2, ccode_conv2))]
             else:
                 statements = []
-            if var3:
-                ccode_conv3 = W.ccode_convert(*node.args.ccode)
+            if fld.vector_type == '3D':
+                ccode_conv3 = fld.W.ccode_convert(*node.args.ccode)
                 statements.append(c.Statement("%s *= %s" % (var3, ccode_conv3)))
-            conv_stat = c.Block(statements)
-            cstat += [c.Assign("err", ccode_eval), conv_stat, c.Statement("CHECKERROR(err)")]
+            cstat += [c.Assign("err", ccode_eval), c.Block(statements)]
+        cstat += [c.Statement("CHECKERROR(err)")]
         node.ccode = c.Block(cstat)
 
     def visit_NestedFieldEvalNode(self, node):
@@ -806,7 +790,7 @@ class KernelGenerator(ast.NodeVisitor):
                               c.Statement("%s *= %s" % (node.var2, ccode_conv2))]
             else:
                 statements = []
-            if node.var3:
+            if fld.vector_type == '3D':
                 ccode_conv3 = fld.W.ccode_convert(*node.args.ccode)
                 statements.append(c.Statement("%s *= %s" % (node.var3, ccode_conv3)))
             cstat += [c.Assign("err", ccode_eval),
diff --git a/parcels/field.py b/parcels/field.py
index 0db69104..cc473732 100644
--- a/parcels/field.py
+++ b/parcels/field.py
@@ -13,7 +13,7 @@ import math
 from .grid import Grid, CGrid, GridCode
 
 
-__all__ = ['Field', 'VectorField', 'SummedField', 'SummedVectorField', 'NestedField']
+__all__ = ['Field', 'VectorField', 'SummedField', 'NestedField']
 
 
 class Field(object):
@@ -939,41 +939,13 @@ class Field(object):
 
     def __add__(self, field):
         if isinstance(self, Field) and isinstance(field, Field):
-            return SummedField([self, field])
+            return SummedField('_SummedField', [self, field])
         elif isinstance(field, SummedField):
+            assert isinstance(self, type(field[0])), 'Fields in a SummedField should be either all scalars or all vectors'
             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, z, y, x, applyConversion=True):
-        tmp = 0
-        for fld in self:
-            tmp += fld.eval(time, z, y, x, 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.
     This enables to interpolate them as one single vector field in the kernels.
@@ -988,12 +960,13 @@ class VectorField(object):
         self.U = U
         self.V = V
         self.W = W
+        self.vector_type = '3D' if W else '2D'
         if self.U.interp_method == 'cgrid_velocity':
             assert self.V.interp_method == 'cgrid_velocity', (
                 'Interpolation methods of U and V are not the same.')
             assert self.U.grid is self.V.grid, (
                 'Grids of U and V are not the same.')
-            if W:
+            if self.vector_type == '3D':
                 assert self.W.interp_method == 'cgrid_velocity', (
                     'Interpolation methods of U and W are not the same.')
                 assert self.W.grid is self.U.grid, (
@@ -1195,7 +1168,7 @@ class VectorField(object):
             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:
+            if self.vector_type == '3D':
                 w = self.W.eval(time, z, y, x, False)
                 w = self.W.units.to_target(w, x, y, z)
                 return (u, v, w)
@@ -1208,7 +1181,7 @@ class VectorField(object):
             if ti < grid.tdim-1 and time > grid.time[ti]:
                 t0 = grid.time[ti]
                 t1 = grid.time[ti + 1]
-                if self.W:
+                if self.vector_type == '3D':
                     (u0, v0, w0) = self.spatial_c_grid_interpolation3D(ti, z, y, x, time)
                     (u1, v1, w1) = self.spatial_c_grid_interpolation3D(ti + 1, z, y, x, time)
                     w = w0 + (w1 - w0) * ((time - t0) / (t1 - t0))
@@ -1217,7 +1190,7 @@ class VectorField(object):
                     (u1, v1) = self.spatial_c_grid_interpolation2D(ti + 1, z, y, x, time)
                 u = u0 + (u1 - u0) * ((time - t0) / (t1 - t0))
                 v = v0 + (v1 - v0) * ((time - t0) / (t1 - t0))
-                if self.W:
+                if self.vector_type == '3D':
                     return (u, v, w)
                 else:
                     return (u, v)
@@ -1225,7 +1198,7 @@ class VectorField(object):
                 # Skip temporal interpolation if time is outside
                 # of the defined time range or if we have hit an
                 # excat value in the time array.
-                if self.W:
+                if self.vector_type == '3D':
                     return self.spatial_c_grid_interpolation3D(ti, z, y, x, grid.time[ti])
                 else:
                     return self.spatial_c_grid_interpolation2D(ti, z, y, x, grid.time[ti])
@@ -1235,7 +1208,7 @@ class VectorField(object):
 
     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:
+        if self.vector_type == '3D':
             return "temporal_interpolationUVW(%s, %s, %s, %s, %s, %s, %s, " \
                    % (x, y, z, t, U.ccode_name, V.ccode_name, W.ccode_name) + \
                    "particle->cxi, particle->cyi, particle->czi, particle->cti, &%s, &%s, &%s, %s)" \
@@ -1253,47 +1226,60 @@ class DeferredArray():
         raise RuntimeError('Field is in deferred_load mode, so can''t be accessed. Use .computeTimeChunk() method to force loading of  data')
 
 
-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.
+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]).
+    SummedField is composed of either Fields or VectorFields.
 
-    :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)
+    :param name: Name of the SummedField
+    :param F: List of fields. F can be a scalar Field, a VectorField, or the zonal component (U) of the VectorField
+    :param V: List of fields defining the meridional component of a VectorField, if F is the zonal component. (default: None)
+    :param W: List of fields defining the vertical component of a VectorField, if F and V are the zonal and meridional components (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, 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, z, y, x)
-                zonal += tmp1
-                meridional += tmp2
-                vertical += tmp3
-            return (zonal, meridional, vertical)
+    def __init__(self, name, F, V=None, W=None):
+        if V is None:
+            if isinstance(F[0], VectorField):
+                vector_type = F[0].vector_type
+            for Fi in F:
+                assert isinstance(Fi, Field) or (isinstance(Fi, VectorField) and Fi.vector_type == vector_type), 'Components of a SummedField must be Field or VectorField'
+                self.append(Fi)
+        elif W is None:
+            for (i, Fi, Vi) in zip(range(len(F)), F, V):
+                assert isinstance(Fi, Field) and isinstance(Vi, Field), \
+                    'F, and V components of a SummedField must be Field'
+                self.append(VectorField(name+'_%d' % i, Fi, Vi))
         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, z, y, x)
-                zonal += tmp1
-                meridional += tmp2
-            return (zonal, meridional)
+            for (i, Fi, Vi, Wi) in zip(range(len(F)), F, V, W):
+                assert isinstance(Fi, Field) and isinstance(Vi, Field) and isinstance(Wi, Field), \
+                    'F, V and W components of a SummedField must be Field'
+                self.append(VectorField(name+'_%d' % i, Fi, Vi, Wi))
+        self.name = name
 
     def __getitem__(self, key):
         if isinstance(key, int):
             return list.__getitem__(self, key)
         else:
-            return self.eval(*key)
+            vals = []
+            val = None
+            for iField in range(len(self)):
+                    val = list.__getitem__(self, iField).eval(*key)
+                    vals.append(val)
+            return tuple(np.sum(vals, 0)) if isinstance(val, tuple) else np.sum(vals)
+
+    def __add__(self, field):
+        if isinstance(field, Field):
+            assert isinstance(self[0], type(field)), 'Fields in a SummedField should be either all scalars or all vectors'
+            self.append(field)
+        elif isinstance(field, SummedField):
+            assert isinstance(self[0], type(field[0])), 'Fields in a SummedField should be either all scalars or all vectors'
+            for fld in field:
+                self.append(fld)
+        return self
 
 
 class NestedField(list):
@@ -1304,22 +1290,29 @@ class NestedField(list):
     NestedField returns an `ErrorOutOfBounds` only if last field is as well out of boundaries.
     NestedField is composed of either Fields or VectorFields.
 
-    :param name: Name of the Nested field
-    :param U: List of fields (order matters). U can be a scalar Field, a VectorField, or the zonal component of the VectorField
-    :param V: List of fields defining the meridional component (default: None)
-    :param W: List of fields defining the vertical component (default: None)
+    :param name: Name of the NestedField
+    :param F: List of fields (order matters). F can be a scalar Field, a VectorField, or the zonal component (U) of the VectorField
+    :param V: List of fields defining the meridional component of a VectorField, if F is the zonal component. (default: None)
+    :param W: List of fields defining the vertical component of a VectorField, if F and V are the zonal and meridional components (default: None)
     """
 
-    def __init__(self, name, U, V=None, W=None):
+    def __init__(self, name, F, V=None, W=None):
         if V is None:
-            for Ui in U:
-                self.append(Ui)
+            if isinstance(F[0], VectorField):
+                vector_type = F[0].vector_type
+            for Fi in F:
+                assert isinstance(Fi, Field) or (isinstance(Fi, VectorField) and Fi.vector_type == vector_type), 'Components of a NestedField must be Field or VectorField'
+                self.append(Fi)
         elif W is None:
-            for (i, Ui, Vi) in zip(range(len(U)), U, V):
-                self.append(VectorField(name+'_%d' % i, Ui, Vi))
+            for (i, Fi, Vi) in zip(range(len(F)), F, V):
+                assert isinstance(Fi, Field) and isinstance(Vi, Field), \
+                    'F, and V components of a NestedField must be Field'
+                self.append(VectorField(name+'_%d' % i, Fi, Vi))
         else:
-            for (i, Ui, Vi, Wi) in zip(range(len(U)), U, V, W):
-                self.append(VectorField(name+'_%d' % i, Ui, Vi, Wi))
+            for (i, Fi, Vi, Wi) in zip(range(len(F)), F, V, W):
+                assert isinstance(Fi, Field) and isinstance(Vi, Field) and isinstance(Wi, Field), \
+                    'F, V and W components of a NestedField must be Field'
+                self.append(VectorField(name+'_%d' % i, Fi, Vi, Wi))
         self.name = name
 
     def __getitem__(self, key):
diff --git a/parcels/fieldset.py b/parcels/fieldset.py
index 627f00eb..ec503515 100644
--- a/parcels/fieldset.py
+++ b/parcels/fieldset.py
@@ -1,4 +1,4 @@
-from parcels.field import Field, VectorField, SummedField, SummedVectorField, NestedField
+from parcels.field import Field, VectorField, SummedField, NestedField
 from parcels.gridset import GridSet
 from parcels.grid import Grid
 from parcels.tools.loggers import logger
@@ -103,15 +103,17 @@ class FieldSet(object):
         :param name: Name of the :class:`parcels.field.Field` object to be added
         """
         name = field.name if name is None else name
+        if hasattr(self, name):  # check if Field with same name already exists when adding new Field
+            raise RuntimeError("FieldSet already has a Field with name '%s'" % name)
         if isinstance(field, SummedField):
             setattr(self, name, field)
+            field.name = name
             for fld in field:
                 self.gridset.add_grid(fld)
                 fld.fieldset = self
         elif isinstance(field, NestedField):
             setattr(self, name, field)
             for fld in field:
-                setattr(self, fld.name, fld)
                 self.gridset.add_grid(fld)
                 fld.fieldset = self
         elif isinstance(field, list):
@@ -130,7 +132,6 @@ class FieldSet(object):
         vfield.fieldset = self
         if isinstance(vfield, NestedField):
             for f in vfield:
-                setattr(self, f.name, f)
                 f.fieldset = self
 
     def check_complete(self):
@@ -151,14 +152,14 @@ class FieldSet(object):
             g.time_origin = self.time_origin
         if not hasattr(self, 'UV'):
             if isinstance(self.U, SummedField):
-                self.add_vector_field(SummedVectorField('UV', self.U, self.V))
+                self.add_vector_field(SummedField('UV', self.U, self.V))
             elif isinstance(self.U, NestedField):
                 self.add_vector_field(NestedField('UV', self.U, self.V))
             else:
                 self.add_vector_field(VectorField('UV', self.U, self.V))
         if not hasattr(self, 'UVW') and hasattr(self, 'W'):
             if isinstance(self.U, SummedField):
-                self.add_vector_field(SummedVectorField('UVW', self.U, self.V, self.W))
+                self.add_vector_field(SummedField('UVW', self.U, self.V, self.W))
             elif isinstance(self.U, NestedField):
                 self.add_vector_field(NestedField('UVW', self.U, self.V, self.W))
             else:
@@ -166,7 +167,7 @@ class FieldSet(object):
 
         ccode_fieldnames = []
         counter = 1
-        for fld in self.fields:
+        for fld in self.get_fields():
                 if fld.name not in ccode_fieldnames:
                     fld.ccode_name = fld.name
                 else:
@@ -496,17 +497,20 @@ class FieldSet(object):
         v = fields.pop('V', None)
         return cls(u, v, fields=fields)
 
-    @property
-    def fields(self):
+    def get_fields(self):
         """Returns a list of all the :class:`parcels.field.Field` and :class:`parcels.field.VectorField`
         objects associated with this FieldSet"""
         fields = []
         for v in self.__dict__.values():
             if type(v) in [Field, VectorField]:
-                fields.append(v)
-            elif isinstance(v, SummedField):
+                if v not in fields:
+                    fields.append(v)
+            elif type(v) in [NestedField, SummedField]:
+                if v not in fields:
+                    fields.append(v)
                 for v2 in v:
-                    fields.append(v2)
+                    if v2 not in fields:
+                        fields.append(v2)
         return fields
 
     def add_constant(self, name, value):
@@ -546,7 +550,7 @@ class FieldSet(object):
         if hasattr(self, 'V'):
             self.V.write(filename, varname='vomecrty')
 
-        for v in self.fields:
+        for v in self.get_fields():
             if (v.name is not 'U') and (v.name is not 'V'):
                 v.write(filename)
 
@@ -562,7 +566,7 @@ class FieldSet(object):
         for gnew in fieldset_new.gridset.grids:
             gnew.advanced = False
 
-        for fnew in fieldset_new.fields:
+        for fnew in fieldset_new.get_fields():
             if isinstance(fnew, VectorField):
                 continue
             f = getattr(self, fnew.name)
@@ -582,16 +586,16 @@ class FieldSet(object):
 
         for g in self.gridset.grids:
             g.update_status = 'not_updated'
-        for f in self.fields:
-            if isinstance(f, VectorField) or not f.grid.defer_load:
+        for f in self.get_fields():
+            if type(f) in [VectorField, NestedField, SummedField] or not f.grid.defer_load:
                 continue
             if f.grid.update_status == 'not_updated':
                 nextTime_loc = f.grid.computeTimeChunk(f, time, signdt)
             nextTime = min(nextTime, nextTime_loc) if signdt >= 0 else max(nextTime, nextTime_loc)
 
         # load in new data
-        for f in self.fields:
-            if isinstance(f, VectorField) or not f.grid.defer_load or f.is_gradient or f.dataFiles is None:
+        for f in self.get_fields():
+            if type(f) in [VectorField, NestedField, SummedField] or not f.grid.defer_load or f.is_gradient or f.dataFiles is None:
                 continue
             g = f.grid
             if g.update_status == 'first_updated':  # First load of data
diff --git a/parcels/include/index_search.h b/parcels/include/index_search.h
index d08c0ec3..999f2abd 100644
--- a/parcels/include/index_search.h
+++ b/parcels/include/index_search.h
@@ -254,7 +254,7 @@ static inline ErrorCode search_indices_curvilinear(type_coord x, type_coord y, t
 
   *xsi = *eta = -1;
   int maxIterSearch = 1e6, it = 0;
-  double tol = 1e-12;
+  double tol = 1e-10;
   while ( (*xsi < -tol) || (*xsi > 1+tol) || (*eta < -tol) || (*eta > 1+tol) ){
     double xgrid_loc[4] = {xgrid[*yi][*xi], xgrid[*yi][*xi+1], xgrid[*yi+1][*xi+1], xgrid[*yi+1][*xi]};
     if (sphere_mesh){ //we are on the sphere
diff --git a/parcels/kernel.py b/parcels/kernel.py
index 416fb4ff..1f4729e8 100644
--- a/parcels/kernel.py
+++ b/parcels/kernel.py
@@ -114,15 +114,12 @@ class Kernel(object):
             self.field_args = kernelgen.field_args
             self.vector_field_args = kernelgen.vector_field_args
             fieldset = self.fieldset
-            for fname in self.vector_field_args:
-                f = getattr(fieldset, fname)
-                Wname = f.W.name if f.W else 'not_defined'
-                for sF in [f.U.name, f.V.name, Wname]:
-                    if sF not in self.field_args:
-                        try:
-                            self.field_args[sF] = getattr(fieldset, sF)
-                        except:
-                            continue
+            for f in self.vector_field_args.values():
+                Wname = f.W.ccode_name if f.W else 'not_defined'
+                for sF_name, sF_component in zip([f.U.ccode_name, f.V.ccode_name, Wname], ['U', 'V', 'W']):
+                    if sF_name not in self.field_args:
+                        if sF_name != 'not_defined':
+                            self.field_args[sF_name] = getattr(f, sF_component)
             self.const_args = kernelgen.const_args
             loopgen = LoopGenerator(fieldset, ptype)
             if path.isfile(c_include):


Test patch:
diff --git a/tests/test_fieldset.py b/tests/test_fieldset.py
index a595d9ad..6fa5a81d 100644
--- a/tests/test_fieldset.py
+++ b/tests/test_fieldset.py
@@ -158,6 +158,25 @@ def test_add_field(xdim, ydim, tmpdir, filename='test_add'):
     fieldset.write(filepath)
 
 
+@pytest.mark.parametrize('dupobject', ['same', 'new'])
+def test_add_duplicate_field(dupobject):
+    data, dimensions = generate_fieldset(100, 100)
+    fieldset = FieldSet.from_data(data, dimensions)
+    field = Field('newfld', fieldset.U.data, lon=fieldset.U.lon, lat=fieldset.U.lat)
+    fieldset.add_field(field)
+    error_thrown = False
+    try:
+        if dupobject == 'same':
+            fieldset.add_field(field)
+        elif dupobject == 'new':
+            field2 = Field('newfld', np.ones((2, 2)), lon=np.array([0, 1]), lat=np.array([0, 2]))
+            fieldset.add_field(field2)
+    except RuntimeError:
+        error_thrown = True
+
+    assert error_thrown
+
+
 @pytest.mark.parametrize('mesh', ['flat', 'spherical'])
 def test_fieldset_celledgesizes(mesh):
     data, dimensions = generate_fieldset(10, 7)
diff --git a/tests/test_fieldset_sampling.py b/tests/test_fieldset_sampling.py
index 051c4bb3..2fa2bfc3 100644
--- a/tests/test_fieldset_sampling.py
+++ b/tests/test_fieldset_sampling.py
@@ -429,7 +429,6 @@ def test_summedfields(mode, with_W, k_sample_p, mesh):
     fieldsetS = FieldSet(U1+U2, V1+V2)
 
     conv = 1852*60 if mesh == 'spherical' else 1.
-    assert np.allclose(fieldsetS.U.eval(0, 0, 0, 0)*conv, 0.3)
     assert np.allclose(fieldsetS.U[0, 0, 0, 0]*conv, 0.3)
 
     P1 = Field('P', 30*np.ones((zdim*gf, ydim*gf, xdim*gf), dtype=np.float32), grid=U1.grid)
@@ -452,6 +451,7 @@ def test_summedfields(mode, with_W, k_sample_p, mesh):
     assert np.isclose(pset[0].p, 60)
     assert np.isclose(pset[0].lon*conv, 0.6, atol=1e-3)
     assert np.isclose(pset[0].lat, 0.9)
+    assert np.allclose(fieldsetS.UV[0][0, 0, 0, 0], [.2/conv, 0])
 
 
 @pytest.mark.parametrize('mode', ['jit', 'scipy'])
@@ -502,3 +502,4 @@ def test_nestedfields(mode, k_sample_p):
     pset.execute(AdvectionRK4+pset.Kernel(k_sample_p), runtime=1, dt=1, recovery={ErrorCode.ErrorOutOfBounds: Recover})
     assert np.isclose(pset[0].lat, -1)
     assert np.isclose(pset[0].p, 999)
+    assert np.allclose(fieldset.UV[0][0, 0, 0, 0], [.1, .2])