Re-ordering of Kernel actions to vectorially sum particle displacements

docs/conf.py

@@ -360,7 +360,8 @@ def linkcode_resolve(domain, info):
     'examples/tutorial_interaction': '_static/pulled_particles_twoatractors_line.gif',
     'examples/documentation_LargeRunsOutput': '_static/harddrive.png',
     'examples/tutorial_unitconverters': '_static/globe-icon.jpg',
-    'examples/documentation_geospatial': '_images/tutorial_geospatial_google_earth.png'
+    'examples/documentation_geospatial': '_images/tutorial_geospatial_google_earth.png',
+    'examples/tutorial_kernelloop': '_static/loop-icon.jpeg',
 }
 # -- Options for LaTeX output ---------------------------------------------
 

docs/documentation.rst

@@ -21,7 +21,6 @@ Parcels has several documentation and tutorial Jupyter notebooks which go throug
    examples/documentation_indexing.ipynb
    examples/tutorial_nemo_curvilinear.ipynb
    examples/tutorial_nemo_3D.ipynb
-   examples/tutorial_SummedFields.ipynb
    examples/tutorial_NestedFields.ipynb
    examples/tutorial_timevaryingdepthdimensions.ipynb
    examples/tutorial_periodic_boundaries.ipynb
@@ -47,6 +46,7 @@ Parcels has several documentation and tutorial Jupyter notebooks which go throug
    examples/tutorial_particle_field_interaction.ipynb
    examples/tutorial_interaction.ipynb
    examples/tutorial_analyticaladvection.ipynb
+   examples/tutorial_kernelloop.ipynb
 
 
 .. nbgallery::

docs/examples/example_brownian.py

@@ -50,7 +50,7 @@ def test_brownian_example(mode, mesh, npart=3000):
     ys = pset.lat * mesh_conversion(mesh)
     xs = pset.lon * mesh_conversion(mesh)  # since near equator, we do not need to care about curvature effect
 
-    tol = 200  # 200m tolerance
+    tol = 250  # 250m tolerance
     assert np.allclose(np.std(xs), expected_std_x, atol=tol)
     assert np.allclose(np.std(ys), expected_std_y, atol=tol)
     assert np.allclose(np.mean(xs), 0, atol=tol)

docs/examples/example_dask_chunk_OCMs.py

@@ -183,13 +183,13 @@
 
     def periodicBC(particle, fieldSet, time):
         if particle.lon > 15.0:
-            particle.lon -= 15.0
+            particle_dlon -= 15.0  # noqa
         if particle.lon < 0:
-            particle.lon += 15.0
+            particle_dlon += 15.0  # noqa
         if particle.lat > 60.0:
-            particle.lat -= 11.0
+            particle_dlat -= 11.0  # noqa
         if particle.lat < 49.0:
-            particle.lat += 11.0
+            particle_dlat += 11.0  # noqa
 
     pset = ParticleSet.from_list(field_set, ptype[mode], lon=lonp, lat=latp)
     pfile = ParticleFile("nemo_particles_chunk", pset, outputdt=delta(days=1))

docs/examples/example_globcurrent.py

@@ -7,7 +7,6 @@
 
 from parcels import (
     AdvectionRK4,
-    ErrorCode,
     Field,
     FieldSet,
     JITParticle,
@@ -68,8 +67,8 @@
     fieldsetall = set_globcurrent_fieldset(files[0:10], deferred_load=False, use_xarray=use_xarray)
     psetall = ParticleSet.from_list(fieldset=fieldsetall, pclass=ptype[mode], lon=[lonstart], lat=[latstart])
     if dt < 0:
-        psetsub[0].time = fieldsetsub.U.grid.time[-1]
-        psetall[0].time = fieldsetall.U.grid.time[-1]
+        psetsub[0].time_nextloop = fieldsetsub.U.grid.time[-1]
+        psetall[0].time_nextloop = fieldsetall.U.grid.time[-1]
 
     psetsub.execute(AdvectionRK4, runtime=delta(days=7), dt=dt)
     psetall.execute(AdvectionRK4, runtime=delta(days=7), dt=dt)
@@ -101,7 +100,7 @@
         fieldset = set_globcurrent_fieldset(time_periodic=delta(days=365), deferred_load=deferred_load)
 
         class MyParticle(ptype[mode]):
-            sample_var = Variable('sample_var', initial=fieldset.U)
+            sample_var = Variable('sample_var', initial=0.)
 
         pset = ParticleSet(fieldset, pclass=MyParticle, lon=25, lat=-35, time=fieldset.U.grid.time[0])
 
@@ -184,20 +183,6 @@
     pset.execute(AdvectionRK4, dt=0.)
 
 
-@pytest.mark.parametrize('mode', ['scipy', 'jit'])
-@pytest.mark.parametrize('dt', [-300, 300])
-@pytest.mark.parametrize('use_xarray', [True, False])
-def test_globcurrent_variable_fromfield(mode, dt, use_xarray):
-    fieldset = set_globcurrent_fieldset(use_xarray=use_xarray)
-
-    class MyParticle(ptype[mode]):
-        sample_var = Variable('sample_var', initial=fieldset.U)
-    time = fieldset.U.grid.time[0] if dt > 0 else fieldset.U.grid.time[-1]
-    pset = ParticleSet(fieldset, pclass=MyParticle, lon=[25], lat=[-35], time=time)
-
-    pset.execute(AdvectionRK4, runtime=delta(days=1), dt=dt)
-
-
 @pytest.mark.parametrize('mode', ['scipy', 'jit'])
 @pytest.mark.parametrize('dt', [-300, 300])
 @pytest.mark.parametrize('with_starttime', [True, False])
@@ -235,10 +220,7 @@
 
     def DeleteP0(particle, fieldset, time):
         if particle.id == 0:
-            return ErrorCode.ErrorOutOfBounds  # we want to pass through recov loop
-
-    def DeleteParticle(particle, fieldset, time):
-        particle.delete()
+            particle.delete()
 
     pset0 = ParticleSet(fieldset, pclass=ptype[mode],
                         lon=[25, 25],
@@ -247,8 +229,7 @@
 
     pset0.execute(pset0.Kernel(DeleteP0)+AdvectionRK4,
                   runtime=delta(days=rundays),
-                  dt=delta(minutes=5),
-                  recovery={ErrorCode.ErrorOutOfBounds: DeleteParticle})
+                  dt=delta(minutes=5))
 
     pset1 = ParticleSet(fieldset, pclass=ptype[mode],
                         lon=[25, 25],
@@ -273,7 +254,7 @@
     pset = ParticleSet(fieldset, pclass=ptype[mode], lon=25, lat=-35)
     pfile = pset.ParticleFile(filename, outputdt=delta(hours=6))
     pset.execute(AdvectionRK4, runtime=delta(days=1), dt=dt, output_file=pfile)
-    pfile.close()
+    pfile.write_latest_locations(pset, max(pset.time_nextloop))
 
     restarttime = np.nanmax if dt > 0 else np.nanmin
     pset_new = ParticleSet.from_particlefile(fieldset, pclass=ptype[mode], filename=filename, restarttime=restarttime)

docs/examples/example_mitgcm.py

@@ -35,9 +35,9 @@
 
     def periodicBC(particle, fieldset, time):
         if particle.lon < 0:
-            particle.lon += fieldset.domain_width
+            particle_dlon += fieldset.domain_width  # noqa
         elif particle.lon > fieldset.domain_width:
-            particle.lon -= fieldset.domain_width
+            particle_dlon -= fieldset.domain_width  # noqa
 
     # Release particles 5 cells away from the Eastern boundary
     pset = ParticleSet.from_line(
@@ -54,7 +54,6 @@
     pset.execute(
         kernels, runtime=delta(days=5), dt=delta(minutes=30), output_file=pfile
     )
-    pfile.close()
 
 
 def test_mitgcm_output_compare():

docs/examples/example_moving_eddies.py

@@ -166,9 +166,12 @@
     pset.execute(method, endtime=0, dt=-dt,
                  output_file=pset.ParticleFile(name=outfile, outputdt=outputdt))
 
+    # Also include last timestep
+    for var in ['lon', 'lat', 'depth', 'time']:
+        pset.collection.setallvardata(f"{var}", pset.collection.getvardata(f"{var}_nextloop"))
+
     assert np.allclose(pset.lon, lons)
     assert np.allclose(pset.lat, lats)
-    return pset
 
 
 @pytest.mark.parametrize('mode', ['scipy', 'jit'])
@@ -177,6 +180,9 @@
     fieldset = moving_eddies_fieldset(mesh=mesh)
     outfile = tmpdir.join("EddyParticle")
     pset = moving_eddies_example(fieldset, outfile, 2, mode=mode)
+    # Also include last timestep
+    for var in ['lon', 'lat', 'depth', 'time']:
+        pset.collection.setallvardata(f"{var}", pset.collection.getvardata(f"{var}_nextloop"))
     if mesh == 'flat':
         assert (pset[0].lon < 2.2e5 and 1.1e5 < pset[0].lat < 1.2e5)
         assert (pset[1].lon < 2.2e5 and 3.7e5 < pset[1].lat < 3.8e5)
@@ -200,6 +206,9 @@
     fieldset = FieldSet.from_parcels(fieldsetfile(mesh, tmpdir), extra_fields={'P': 'P'})
     outfile = tmpdir.join("EddyParticle")
     pset = moving_eddies_example(fieldset, outfile, 2, mode=mode)
+    # Also include last timestep
+    for var in ['lon', 'lat', 'depth', 'time']:
+        pset.collection.setallvardata(f"{var}", pset.collection.getvardata(f"{var}_nextloop"))
     if mesh == 'flat':
         assert (pset[0].lon < 2.2e5 and 1.1e5 < pset[0].lat < 1.2e5)
         assert (pset[1].lon < 2.2e5 and 3.7e5 < pset[1].lat < 3.8e5)
@@ -226,17 +235,17 @@
 
     def periodicBC(particle, fieldset, time):
         if particle.lon < fieldset.halo_west:
-            particle.lon += fieldset.halo_east - fieldset.halo_west
+            particle_dlon += fieldset.halo_east - fieldset.halo_west  # noqa
         elif particle.lon > fieldset.halo_east:
-            particle.lon -= fieldset.halo_east - fieldset.halo_west
+            particle_dlon -= fieldset.halo_east - fieldset.halo_west  # noqa
         if particle.lat < fieldset.halo_south:
-            particle.lat += fieldset.halo_north - fieldset.halo_south
+            particle_dlat += fieldset.halo_north - fieldset.halo_south  # noqa
         elif particle.lat > fieldset.halo_north:
-            particle.lat -= fieldset.halo_north - fieldset.halo_south
+            particle_dlat -= fieldset.halo_north - fieldset.halo_south  # noqa
 
     def slowlySouthWestward(particle, fieldset, time):
-        particle.lon = particle.lon - 5 * particle.dt / 1e5
-        particle.lat -= 3 * particle.dt / 1e5
+        particle_dlon -= 5 * particle.dt / 1e5  # noqa
+        particle_dlat -= 3 * particle.dt / 1e5  # noqa
 
     kernels = pset.Kernel(AdvectionRK4)+slowlySouthWestward+periodicBC
     pset.execute(kernels, runtime=delta(days=6), dt=delta(hours=1))

docs/examples/example_nemo_curvilinear.py

@@ -49,7 +49,7 @@
 
     def periodicBC(particle, fieldSet, time):
         if particle.lon > 180:
-            particle.lon -= 360
+            particle_dlon -= 360  # noqa
 
     pset = ParticleSet.from_list(field_set, ptype[mode], lon=lonp, lat=latp)
     pfile = ParticleFile(outfile, pset, outputdt=delta(days=1))

docs/examples/example_peninsula.py

@@ -98,6 +98,8 @@
 
 
 def UpdateP(particle, fieldset, time):
+    if time == 0:
+        particle.p_start = fieldset.P[time, particle.depth, particle.lat, particle.lon]
     particle.p = fieldset.P[time, particle.depth, particle.lat, particle.lon]
 
 
@@ -132,7 +134,7 @@
         # JIT compilation requires a-priori knowledge of the particle
         # data structure, so we define additional variables here.
         p = Variable('p', dtype=np.float32, initial=0.)
-        p_start = Variable('p_start', dtype=np.float32, initial=fieldset.P)
+        p_start = Variable('p_start', dtype=np.float32, initial=0)
 
     # Initialise particles
     if fieldset.U.grid.mesh == 'flat':

docs/examples/example_stommel.py

@@ -22,7 +22,7 @@
     timer,
 )
 
-pset_modes = ['soa', 'aos']
+pset_modes = ['soa']
 ptype = {'scipy': ScipyParticle, 'jit': JITParticle}
 method = {'RK4': AdvectionRK4, 'EE': AdvectionEE, 'RK45': AdvectionRK45}
 pset_type = {'soa': {'pset': ParticleSetSOA, 'pfile': ParticleFileSOA, 'kernel': KernelSOA},
@@ -77,6 +77,8 @@
 
 
 def UpdateP(particle, fieldset, time):
+    if time == 0:
+        particle.p_start = fieldset.P[time, particle.depth, particle.lat, particle.lon]
     particle.p = fieldset.P[time, particle.depth, particle.lat, particle.lon]
 
 
@@ -100,9 +102,15 @@
     timer.psetinit = timer.Timer('Pset_init', parent=timer.pset)
     ParticleClass = JITParticle if mode == 'jit' else ScipyParticle
 
+    # Execute for 600 days, with 1-hour timesteps and 5-day output
+    runtime = delta(days=600)
+    dt = delta(hours=1)
+    outputdt = delta(days=5)
+
     class MyParticle(ParticleClass):
         p = Variable('p', dtype=np.float32, initial=0.)
-        p_start = Variable('p_start', dtype=np.float32, initial=fieldset.P)
+        p_start = Variable('p_start', dtype=np.float32, initial=0.)
+        next_dt = Variable('next_dt', dtype=np.float64, initial=dt.total_seconds())
         age = Variable('age', dtype=np.float32, initial=0.)
 
     pset = pset_type[pset_mode]['pset'].from_line(fieldset, size=npart, pclass=MyParticle, repeatdt=repeatdt,
@@ -111,10 +119,6 @@
     if verbose:
         print(f"Initial particle positions:\n{pset}")
 
-    # Execute for 30 days, with 1hour timesteps and 12-hourly output
-    runtime = delta(days=600)
-    dt = delta(hours=1)
-    outputdt = delta(days=5)
     maxage = runtime.total_seconds() if maxage is None else maxage
     fieldset.add_constant('maxage', maxage)
     print("Stommel: Advecting %d particles for %s" % (npart, runtime))