Merge pull request #633 from OceanParcels/no_field_gradient

Removing Field.gradient() function

parcels/examples/example_globcurrent.py

@@ -175,19 +175,3 @@ class MyParticle(ptype[mode]):
     pset = ParticleSet(fieldset, pclass=MyParticle, lon=[25], lat=[-35], time=time)
 
     pset.execute(AdvectionRK4, runtime=delta(days=1), dt=dt)
-
-
-@pytest.mark.parametrize('deferred_load', [True, False])
-@pytest.mark.parametrize('use_xarray', [True, False])
-def test_globcurrent_deferred_fieldset_gradient(deferred_load, use_xarray):
-    fieldset = set_globcurrent_fieldset(deferred_load=deferred_load, use_xarray=use_xarray)
-    (dU_dx, dU_dy) = fieldset.U.gradient()
-    fieldset.add_field(dU_dy)
-
-    pset = ParticleSet(fieldset, pclass=JITParticle, lon=25, lat=-35)
-    pset.execute(AdvectionRK4, runtime=delta(days=1), dt=delta(days=1))
-
-    tdim = 3 if deferred_load else 366
-    assert(dU_dx.data.shape == (tdim, 41, 81))
-    assert(fieldset.dU_dy.data.shape == (tdim, 41, 81))
-    assert(dU_dx is fieldset.U.gradientx)

parcels/field.py

@@ -122,8 +122,6 @@ def __init__(self, name, data, lon=None, lat=None, depth=None, time=None, grid=N
                 self.data[self.data > self.vmax] = 0.
 
         self._scaling_factor = None
-        (self.gradientx, self.gradienty) = (None, None)  # to store if Field is a gradient() of another field
-        self.is_gradient = False
 
         # Variable names in JIT code
         self.dimensions = kwargs.pop('dimensions', None)
@@ -413,34 +411,6 @@ def cell_areas(self):
             self.calc_cell_edge_sizes()
         return self.grid.cell_edge_sizes['x'] * self.grid.cell_edge_sizes['y']
 
-    def gradient(self, update=False, tindex=None):
-        """Method to calculate horizontal gradients of Field.
-                Returns two Fields: the zonal and meridional gradients,
-                on the same Grid as the original Field, using numpy.gradient() method
-                Names of these grids are dNAME_dx and dNAME_dy, where NAME is the name
-                of the original Field"""
-        tindex = range(self.grid.tdim) if tindex is None else tindex
-        if not self.grid.cell_edge_sizes:
-            self.calc_cell_edge_sizes()
-        if self.grid.defer_load and isinstance(self.data, DeferredArray):
-            (dFdx, dFdy) = (None, None)
-        else:
-            dFdy = np.gradient(self.data[tindex, :], axis=-2) / self.grid.cell_edge_sizes['y']
-            dFdx = np.gradient(self.data[tindex, :], axis=-1) / self.grid.cell_edge_sizes['x']
-        if update:
-            if self.gradientx.data is None:
-                self.gradientx.data = np.zeros_like(self.data)
-                self.gradienty.data = np.zeros_like(self.data)
-            self.gradientx.data[tindex, :] = dFdx
-            self.gradienty.data[tindex, :] = dFdy
-        else:
-            dFdx_fld = Field('d%s_dx' % self.name, dFdx, grid=self.grid)
-            dFdy_fld = Field('d%s_dy' % self.name, dFdy, grid=self.grid)
-            dFdx_fld.is_gradient = True
-            dFdy_fld.is_gradient = True
-            (self.gradientx, self.gradienty) = (dFdx_fld, dFdy_fld)
-            return (dFdx_fld, dFdy_fld)
-
     def search_indices_vertical_z(self, z):
         grid = self.grid
         z = np.float32(z)

parcels/fieldset.py

@@ -774,7 +774,7 @@ def computeTimeChunk(self, time, dt):
 
         # load in new data
         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:
+            if type(f) in [VectorField, NestedField, SummedField] or not f.grid.defer_load or f.dataFiles is None:
                 continue
             g = f.grid
             if g.update_status == 'first_updated':  # First load of data
@@ -805,8 +805,6 @@ def computeTimeChunk(self, time, dt):
                     f.data[tind, :] = np.where(f.data[tind, :] < f.vmin, 0, f.data[tind, :])
                 if f.vmax is not None:
                     f.data[tind, :] = np.where(f.data[tind, :] > f.vmax, 0, f.data[tind, :])
-                if f.gradientx is not None:
-                    f.gradient(update=True, tindex=tind)
 
         # do user-defined computations on fieldset data
         if self.compute_on_defer:

parcels/kernels/diffusion.py

@@ -30,8 +30,14 @@ def SpatiallyVaryingBrownianMotion2D(particle, fieldset, time):
     Rx = random.uniform(-1., 1.) * math.sqrt(2*math.fabs(particle.dt)*kh_zonal/r)
 
     # Deterministic 'boost' out of areas of low diffusivity
-    dKdx = fieldset.dKh_zonal_dx[time, particle.depth, particle.lat, particle.lon]
-    dKdy = fieldset.dKh_meridional_dy[time, particle.depth, particle.lat, particle.lon]
+    dx = .01  # for spherical coords, dx is in degrees
+    Kyp1 = fieldset.Kh_meridional[time, particle.depth, particle.lat+dx, particle.lon]
+    Kym1 = fieldset.Kh_meridional[time, particle.depth, particle.lat-dx, particle.lon]
+    dKdy = (Kyp1-Kym1) / (2*dx)
+    Kxp1 = fieldset.Kh_zonal[time, particle.depth, particle.lat, particle.lon+dx]
+    Kxm1 = fieldset.Kh_zonal[time, particle.depth, particle.lat, particle.lon-dx]
+    dKdx = (Kxp1-Kxm1) / (2*dx)
+
     CorrectionX = dKdx * math.fabs(particle.dt)
     CorrectionY = dKdy * math.fabs(particle.dt)
 

tests/test_diffusion.py

@@ -1,6 +1,5 @@
 from parcels import (FieldSet, Field, RectilinearZGrid, ParticleSet, BrownianMotion2D,
-                     SpatiallyVaryingBrownianMotion2D, JITParticle, ScipyParticle,
-                     Geographic, GeographicPolar)
+                     SpatiallyVaryingBrownianMotion2D, JITParticle, ScipyParticle)
 from parcels import rng as random
 from datetime import timedelta as delta
 import numpy as np
@@ -71,14 +70,6 @@ def test_fieldKh_SpatiallyVaryingBrownianMotion(mesh, mode, xdim=200, ydim=100):
     fieldset.add_field(Field('Kh_zonal', Kh, grid=grid))
     fieldset.add_field(Field('Kh_meridional', Kh, grid=grid))
 
-    dKh_zonal_dx, _ = fieldset.Kh_zonal.gradient()
-    _, dKh_meridional_dy = fieldset.Kh_meridional.gradient()
-    fieldset.add_field(dKh_zonal_dx)
-    fieldset.add_field(dKh_meridional_dy)
-    if mesh == 'spherical':
-        fieldset.dKh_zonal_dx.units = GeographicPolar()
-        fieldset.dKh_meridional_dy.units = Geographic()
-
     npart = 100
     runtime = delta(days=1)
 

tests/test_fieldset.py

@@ -5,7 +5,6 @@
 import datetime
 import numpy as np
 import xarray as xr
-import math
 import pytest
 from os import path
 import cftime
@@ -289,26 +288,6 @@ def test_fieldset_cellareas(mesh):
             assert np.allclose(cell_areas[y, :], cell_areas[y, 0], rtol=1e-3)
 
 
-@pytest.mark.parametrize('mesh', ['flat', 'spherical'])
-def test_fieldset_gradient(mesh):
-    data, dimensions = generate_fieldset(5, 3)
-    fieldset = FieldSet.from_data(data, dimensions, mesh=mesh)
-
-    # Calculate field gradients for testing against numpy gradients.
-    dFdx, dFdy = fieldset.V.gradient()
-
-    # Create numpy fields.
-    conv_factor = 6.371e6 * np.pi / 180. if mesh == 'spherical' else 1.
-    np_dFdx = np.gradient(fieldset.V.data[0, :, :], (np.diff(fieldset.V.lon) * conv_factor)[0], axis=1)
-    np_dFdy = np.gradient(fieldset.V.data[0, :, :], (np.diff(fieldset.V.lat) * conv_factor)[0], axis=0)
-    if mesh == 'spherical':
-        for y in range(np_dFdx.shape[0]):
-            np_dFdx[:, y] /= math.cos(fieldset.V.grid.lat[y] * math.pi / 180.)
-
-    assert np.allclose(dFdx.data, np_dFdx, rtol=5e-2)  # Field gradient dx.
-    assert np.allclose(dFdy.data, np_dFdy, rtol=5e-2)  # Field gradient dy.
-
-
 def addConst(particle, fieldset, time):
     particle.lon = particle.lon + fieldset.movewest + fieldset.moveeast
 
@@ -461,8 +440,6 @@ def test_fieldset_defer_loading_function(zdim, scale_fac, tmpdir, filename='test
     # testing for scaling factors
     fieldset.U.set_scaling_factor(scale_fac)
 
-    dFdx, dFdy = fieldset.V.gradient()
-
     dz = np.gradient(fieldset.U.depth)
     DZ = np.moveaxis(np.tile(dz, (fieldset.U.grid.ydim, fieldset.U.grid.xdim, 1)), [0, 1, 2], [1, 2, 0])
 
@@ -475,7 +452,6 @@ def compute(fieldset):
     fieldset.compute_on_defer = compute
     fieldset.computeTimeChunk(1, 1)
     assert np.allclose(fieldset.U.data, scale_fac*(zdim-1.)/zdim)
-    assert np.allclose(dFdx.data, 0)
 
     pset = ParticleSet(fieldset, JITParticle, 0, 0)
 
@@ -484,7 +460,6 @@ def DoNothing(particle, fieldset, time):
 
     pset.execute(DoNothing, dt=3600)
     assert np.allclose(fieldset.U.data, scale_fac*(zdim-1.)/zdim)
-    assert np.allclose(dFdx.data, 0)
 
 
 @pytest.mark.parametrize('maxlatind', [3, pytest.param(2, marks=pytest.mark.xfail(strict=True))])