diff --git a/parcels/examples/example_globcurrent.py b/parcels/examples/example_globcurrent.py index c4a3e7314..b2ad2f6ef 100644 --- a/parcels/examples/example_globcurrent.py +++ b/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) diff --git a/parcels/field.py b/parcels/field.py index ec6d9a422..e389c94d8 100644 --- a/parcels/field.py +++ b/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) diff --git a/parcels/fieldset.py b/parcels/fieldset.py index c5b084048..e36a5f812 100644 --- a/parcels/fieldset.py +++ b/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: diff --git a/parcels/kernels/diffusion.py b/parcels/kernels/diffusion.py index b97ded9bb..2bcec7f99 100644 --- a/parcels/kernels/diffusion.py +++ b/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) diff --git a/tests/test_diffusion.py b/tests/test_diffusion.py index 3d7e999bd..99d7a03f1 100644 --- a/tests/test_diffusion.py +++ b/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) diff --git a/tests/test_fieldset.py b/tests/test_fieldset.py index 89e9ed606..925907690 100644 --- a/tests/test_fieldset.py +++ b/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))])