Merge pull request #319 from OceanParcels/auto_advanceTime2

Auto advance time

parcels/examples/example_globcurrent.py

@@ -9,24 +9,22 @@
 ptype = {'scipy': ScipyParticle, 'jit': JITParticle}
 
 
-def set_globcurrent_fieldset(filename=None, indices={}):
+def set_globcurrent_fieldset(filename=None, indices={}, full_load=False):
     if filename is None:
         filename = path.join(path.dirname(__file__), 'GlobCurrent_example_data',
                              '20*-GLOBCURRENT-L4-CUReul_hs-ALT_SUM-v02.0-fv01.0.nc')
     filenames = {'U': filename, 'V': filename}
     variables = {'U': 'eastward_eulerian_current_velocity', 'V': 'northward_eulerian_current_velocity'}
     dimensions = {'lat': 'lat', 'lon': 'lon', 'time': 'time'}
-    return FieldSet.from_netcdf(filenames, variables, dimensions, indices)
+    return FieldSet.from_netcdf(filenames, variables, dimensions, indices, full_load=full_load)
 
 
 def test_globcurrent_fieldset():
     fieldset = set_globcurrent_fieldset()
     assert(fieldset.U.lon.size == 81)
     assert(fieldset.U.lat.size == 41)
-    assert(fieldset.U.data.shape == (365, 41, 81))
     assert(fieldset.V.lon.size == 81)
     assert(fieldset.V.lat.size == 41)
-    assert(fieldset.V.data.shape == (365, 41, 81))
 
     indices = {'lon': [5], 'lat': range(20, 30)}
     fieldsetsub = set_globcurrent_fieldset(indices=indices)
@@ -46,19 +44,17 @@ def test_globcurrent_fieldset_advancetime(mode, dt, substart, subend, lonstart,
                          '20*-GLOBCURRENT-L4-CUReul_hs-ALT_SUM-v02.0-fv01.0.nc')
     files = sorted(glob(str(basepath)))
 
-    fieldsetsub = set_globcurrent_fieldset(files[substart:subend])
+    fieldsetsub = set_globcurrent_fieldset(files[0:10])
     psetsub = ParticleSet.from_list(fieldset=fieldsetsub, pclass=ptype[mode], lon=[lonstart], lat=[latstart])
 
-    fieldsetall = set_globcurrent_fieldset(files[0:10])
+    fieldsetall = set_globcurrent_fieldset(files[0:10], full_load=True)
     psetall = ParticleSet.from_list(fieldset=fieldsetall, pclass=ptype[mode], lon=[lonstart], lat=[latstart])
     if dt < 0:
         psetsub[0].time = fieldsetsub.U.time[-1]
         psetall[0].time = fieldsetall.U.time[-1]
 
     for i in irange:
         psetsub.execute(AdvectionRK4, runtime=delta(days=1), dt=dt)
-        fieldsetsub.advancetime(set_globcurrent_fieldset(files[i]))
-
         psetall.execute(AdvectionRK4, runtime=delta(days=1), dt=dt)
 
     assert abs(psetsub[0].lon - psetall[0].lon) < 1e-4

parcels/examples/example_moving_eddies.py

@@ -148,6 +148,34 @@ def test_moving_eddies_file(fieldsetfile, mode):
     assert(pset[1].lon < 2.0 and 48.8 < pset[1].lat < 48.85)
 
 
+@pytest.mark.parametrize('mode', ['scipy', 'jit'])
+def test_periodic_and_computeTimeChunk_eddies(mode):
+    filename = path.join(path.dirname(__file__), 'MovingEddies_data', 'moving_eddies')
+    fieldset = FieldSet.from_parcels(filename)
+    fieldset.add_periodic_halo(zonal=True, meridional=True)
+    pset = ParticleSet.from_list(fieldset=fieldset,
+                                 pclass=ptype[mode],
+                                 lon=[3.3, 3.3],
+                                 lat=[46.0, 47.8])
+
+    def periodicBC(particle, fieldset, time, dt):
+        if particle.lon < fieldset.halo_west:
+            particle.lon += fieldset.halo_east - fieldset.halo_west
+        elif particle.lon > fieldset.halo_east:
+            particle.lon -= fieldset.halo_east - fieldset.halo_west
+        if particle.lat < fieldset.halo_south:
+            particle.lat += fieldset.halo_north - fieldset.halo_south
+        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
+
+    kernels = pset.Kernel(AdvectionRK4)+slowlySouthWestward+periodicBC
+    pset.execute(kernels, runtime=delta(days=6), dt=delta(hours=1))
+
+
 if __name__ == "__main__":
     p = ArgumentParser(description="""
 Example of particle advection around an idealised peninsula""")

parcels/examples/example_ofam.py

@@ -7,27 +7,17 @@
 ptype = {'scipy': ScipyParticle, 'jit': JITParticle}
 
 
-def set_ofam_fieldset():
+def set_ofam_fieldset(full_load=False):
     filenames = {'U': path.join(path.dirname(__file__), 'OFAM_example_data', 'OFAM_simple_U.nc'),
                  'V': path.join(path.dirname(__file__), 'OFAM_example_data', 'OFAM_simple_V.nc')}
     variables = {'U': 'u', 'V': 'v'}
     dimensions = {'lat': 'yu_ocean', 'lon': 'xu_ocean', 'depth': 'st_ocean',
                   'time': 'Time'}
-    return FieldSet.from_netcdf(filenames, variables, dimensions, allow_time_extrapolation=True)
-
-
-def test_ofam_fieldset():
-    fieldset = set_ofam_fieldset()
-    assert(fieldset.U.lon.size == 2001)
-    assert(fieldset.U.lat.size == 601)
-    assert(fieldset.U.data.shape == (4, 601, 2001))
-    assert(fieldset.V.lon.size == 2001)
-    assert(fieldset.V.lat.size == 601)
-    assert(fieldset.V.data.shape == (4, 601, 2001))
+    return FieldSet.from_netcdf(filenames, variables, dimensions, allow_time_extrapolation=True, full_load=full_load)
 
 
 def test_ofam_fieldset_fillvalues():
-    fieldset = set_ofam_fieldset()
+    fieldset = set_ofam_fieldset(full_load=True)
     # V.data[0, 0, 150] is a landpoint, that makes NetCDF4 generate a masked array, instead of an ndarray
     assert(fieldset.V.data[0, 0, 150] == 0)