Skip to content

Commit 2ce77c1

Browse files
Merge pull request #2212 from OceanParcels/removing_nextloop_variables
Removing `_nextloop` particle variables
2 parents 29a4309 + bc6394c commit 2ce77c1

File tree

7 files changed

+33
-54
lines changed

7 files changed

+33
-54
lines changed

docs/examples/tutorial_kernelloop.ipynb

Lines changed: 11 additions & 25 deletions
Original file line numberDiff line numberDiff line change
@@ -26,17 +26,7 @@
2626
"\n",
2727
"When you run a Parcels simulation (i.e. a call to `pset.execute()`), the Kernel loop is the main part of the code that is executed. This part of the code loops through all particles and executes the Kernels that are defined for each particle.\n",
2828
"\n",
29-
"In order to make sure that the displacements of a particle in the different Kernels can be summed, all Kernels add to a _change_ in position (`particle_dlon`, `particle_dlat`, and `particle_ddepth`). This is important, because there are situations where movement kernels would otherwise not commute. Take the example of advecting particles by currents _and_ winds. If the particle would first be moved by the currents and then by the winds, the result could be different from first moving by the winds and then by the currents. Instead, by adding the changes in position, the ordering of the Kernels has no consequence on the particle displacement."
30-
]
31-
},
32-
{
33-
"cell_type": "markdown",
34-
"metadata": {},
35-
"source": [
36-
"<div class=\"alert alert-info\">\n",
37-
"\n",
38-
"__Note__ that the variables `particle_dlon`, `particle_dlat`, and `particle_ddepth` are defined by the wrapper function that Parcels generates for each Kernel. This is why you don't have to define these variables yourself when writing a Kernel. See [here](https://github.com/OceanParcels/parcels/blob/daa4b062ed8ae0b2be3d87367d6b45599d6f95db/parcels/kernel.py#L277-L294) for the implementation of the wrapper functions.\n",
39-
"</div>"
29+
"In order to make sure that the displacements of a particle in the different Kernels can be summed, all Kernels add to a _change_ in position (`particles.dlon`, `particles.dlat`, and `particles.ddepth`). This is important, because there are situations where movement kernels would otherwise not commute. Take the example of advecting particles by currents _and_ winds. If the particle would first be moved by the currents and then by the winds, the result could be different from first moving by the winds and then by the currents. Instead, by adding the changes in position, the ordering of the Kernels has no consequence on the particle displacement."
4030
]
4131
},
4232
{
@@ -52,29 +42,25 @@
5242
"source": [
5343
"Below is a structured overview of the Kernel loop is implemented. Note that this is for longitude only, but the same process is applied for latitude and depth.\n",
5444
"\n",
55-
"1. Define an extra variable `particle.lon_nextloop` for each particle, which is the longitude at the end of the Kernel loop. Inititalise it to `particle.lon`.\n",
56-
"\n",
57-
"2. Also define an extra variable `particle.time_nextloop` for each particle, which is the time at the end of the Kernel loop. Inititalise it to `particle.time`.\n",
45+
"1. Initialise an extra Variable `particles.lon=0` and `particles.time_nextloop = particles.time`\n",
5846
"\n",
59-
"3. Within the Kernel loop, for each particle:<br>\n",
47+
"2. Within the Kernel loop, for each particle:<br>\n",
6048
"\n",
61-
" 1. Update `particle.lon` with `particle.lon_nextloop`<br>\n",
49+
" 1. Update `particles.lon += particles.dlon`<br>\n",
6250
"\n",
63-
" 2. Update `particle.time` with `particle.time_nextloop`<br>\n",
51+
" 2. Set variable `particles.dlon = 0`<br>\n",
6452
"\n",
65-
" 3. Set local variable `particle_dlon = 0`<br>\n",
53+
" 3. Update `particles.time = particles.time_nextloop`\n",
6654
"\n",
6755
" 4. For each Kernel in the list of Kernels:\n",
6856
" \n",
6957
" 1. Execute the Kernel\n",
7058
" \n",
71-
" 2. Update `particle_dlon` by adding the change in longitude, if needed<br>\n",
59+
" 2. Update `particles.dlon` by adding the change in longitude, if needed<br>\n",
7260
"\n",
73-
" 5. Update `particle.lon_nextloop` with `particle.lon + particle_dlon`<br>\n",
74-
" \n",
75-
" 6. Update `particle.time_nextloop` with `particle.time + particle.dt`<br>\n",
61+
" 5. Update `particles.time_nextloop += particles.dt`<br>\n",
7662
"\n",
77-
" 7. If `outputdt` is a multiple of `particle.time`, write `particle.lon` and `particle.time` to zarr output file<br>\n",
63+
" 6. If `outputdt` is a multiple of `particle.time`, write `particle.lon` and `particle.time` to zarr output file<br>\n",
7864
"\n",
7965
"Besides having commutable Kernels, the main advantage of this implementation is that, when using Field Sampling with e.g. `particle.temp = fieldset.Temp[particle.time, particle.depth, particle.lat, particle.lon]`, the particle location stays the same throughout the entire Kernel loop. Additionally, this implementation ensures that the particle location is the same as the location of the sampled field in the output file."
8066
]
@@ -268,12 +254,12 @@
268254
"### 1. Avoid updating particle locations directly in Kernels\n",
269255
"It is better not to update `particle.lon` directly in a Kernel, as it can interfere with the loop above. Assigning a value to `particle.lon` in a Kernel will throw a warning. \n",
270256
"\n",
271-
"Instead, update the local variable `particle_dlon`.\n",
257+
"Instead, update the local variable `particle.dlon`.\n",
272258
"\n",
273259
"### 2. Be careful with updating particle variables that do not depend on Fields.\n",
274260
"While assigning the interpolated value of a `Field` to a Particle goes well in the loop above, this is not necessarily so for assigning other attributes. For example, a line like `particle.age += particle.dt` is executed directly so may result in the age being `dt` at `time = 0` in the output file. \n",
275261
"\n",
276-
"A workaround is to either initialise the age to `-dt`, or to increase the `age` only when `particle.time > 0` (using an `if` statement).\n",
262+
"A workaround is to either initialise the age to `-dt`, or to increase the `age` only when `particle.time > 0` (using an `np.where` statement).\n",
277263
"\n",
278264
"\n",
279265
"### 3. The last time is not written to file\n",

parcels/kernel.py

Lines changed: 7 additions & 8 deletions
Original file line numberDiff line numberDiff line change
@@ -118,21 +118,20 @@ def add_positionupdate_kernels(self):
118118
def Setcoords(particles, fieldset): # pragma: no cover
119119
import numpy as np # noqa
120120

121+
particles.lon += particles.dlon
122+
particles.lat += particles.dlat
123+
particles.depth += particles.ddepth
124+
121125
particles.dlon = 0
122126
particles.dlat = 0
123127
particles.ddepth = 0
124-
particles.lon = particles.lon_nextloop
125-
particles.lat = particles.lat_nextloop
126-
particles.depth = particles.depth_nextloop
128+
127129
particles.time = particles.time_nextloop
128130

129-
def Updatecoords(particles, fieldset): # pragma: no cover
130-
particles.lon_nextloop = particles.lon + particles.dlon
131-
particles.lat_nextloop = particles.lat + particles.dlat
132-
particles.depth_nextloop = particles.depth + particles.ddepth
131+
def UpdateTime(particles, fieldset): # pragma: no cover
133132
particles.time_nextloop = particles.time + particles.dt
134133

135-
self._pyfuncs = (Setcoords + self + Updatecoords)._pyfuncs
134+
self._pyfuncs = (Setcoords + self + UpdateTime)._pyfuncs
136135

137136
def check_fieldsets_in_kernels(self, pyfunc): # TODO v4: this can go into another method? assert_is_compatible()?
138137
"""

parcels/particle.py

Lines changed: 0 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -170,13 +170,11 @@ def get_default_particle(spatial_dtype: np.float32 | np.float64) -> ParticleClas
170170
dtype=spatial_dtype,
171171
attrs={"standard_name": "longitude", "units": "degrees_east", "axis": "X"},
172172
),
173-
Variable("lon_nextloop", dtype=spatial_dtype, to_write=False),
174173
Variable(
175174
"lat",
176175
dtype=spatial_dtype,
177176
attrs={"standard_name": "latitude", "units": "degrees_north", "axis": "Y"},
178177
),
179-
Variable("lat_nextloop", dtype=spatial_dtype, to_write=False),
180178
Variable(
181179
"depth",
182180
dtype=spatial_dtype,
@@ -185,7 +183,6 @@ def get_default_particle(spatial_dtype: np.float32 | np.float64) -> ParticleClas
185183
Variable("dlon", dtype=spatial_dtype, to_write=False),
186184
Variable("dlat", dtype=spatial_dtype, to_write=False),
187185
Variable("ddepth", dtype=spatial_dtype, to_write=False),
188-
Variable("depth_nextloop", dtype=spatial_dtype, to_write=False),
189186
Variable(
190187
"time",
191188
dtype=_SAME_AS_FIELDSET_TIME_INTERVAL.VALUE,

parcels/particlefile.py

Lines changed: 3 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -262,9 +262,9 @@ def write_latest_locations(self, pset, time):
262262
time :
263263
Time at which to write ParticleSet. Note that typically this would be pset.time_nextloop
264264
"""
265-
for var in ["lon", "lat", "depth", "time"]:
266-
pset._data[f"{var}"] = pset._data[f"{var}_nextloop"]
267-
265+
for var in ["lon", "lat", "depth"]:
266+
pset._data[f"{var}"] += pset._data[f"d{var}"]
267+
pset._data["time"] = pset._data["time_nextloop"]
268268
self.write(pset, time)
269269

270270

parcels/particleset.py

Lines changed: 0 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -130,9 +130,6 @@ def __init__(
130130
lat=lat,
131131
depth=depth,
132132
time=time,
133-
lon_nextloop=lon,
134-
lat_nextloop=lat,
135-
depth_nextloop=depth,
136133
time_nextloop=time,
137134
trajectory=trajectory_ids,
138135
),

tests/v4/test_advection.py

Lines changed: 11 additions & 11 deletions
Original file line numberDiff line numberDiff line change
@@ -102,8 +102,8 @@ def test_advection_zonal_periodic():
102102
pset = ParticleSet(fieldset, pclass=PeriodicParticle, lon=startlon, lat=[0.5, 0.5])
103103
pset.execute([AdvectionEE, periodicBC], runtime=np.timedelta64(40, "s"), dt=np.timedelta64(1, "s"))
104104
np.testing.assert_allclose(pset.total_dlon, 4, atol=1e-5)
105-
np.testing.assert_allclose(pset.lon_nextloop, startlon, atol=1e-5)
106-
np.testing.assert_allclose(pset.lat_nextloop, 0.5, atol=1e-5)
105+
np.testing.assert_allclose(pset.lon + pset.dlon, startlon, atol=1e-5)
106+
np.testing.assert_allclose(pset.lat, 0.5, atol=1e-5)
107107

108108

109109
def test_horizontal_advection_in_3D_flow(npart=10):
@@ -246,7 +246,7 @@ def test_radialrotation(npart=10):
246246
pset = parcels.ParticleSet(fieldset, lon=lon, lat=lat, time=starttime)
247247
pset.execute(parcels.AdvectionRK4, endtime=np.timedelta64(10, "m"), dt=dt)
248248

249-
theta = 2 * np.pi * (pset.time_nextloop - starttime) / np.timedelta64(24 * 3600, "s")
249+
theta = 2 * np.pi * (pset.time - starttime) / np.timedelta64(24 * 3600, "s")
250250
true_lon = (lon - 30.0) * np.cos(theta) + 30.0
251251
true_lat = -(lon - 30.0) * np.sin(theta) + 30.0
252252

@@ -300,11 +300,11 @@ def truth_moving(x_0, y_0, t):
300300
lon = x_0 + ds.u_g * t + (ds.u_0 - ds.u_g) / ds.f * np.sin(ds.f * t)
301301
return lon, lat
302302

303-
exp_lon, exp_lat = truth_moving(start_lon, start_lat, pset.time_nextloop[0])
304-
np.testing.assert_allclose(pset.lon_nextloop, exp_lon, rtol=rtol)
305-
np.testing.assert_allclose(pset.lat_nextloop, exp_lat, rtol=rtol)
303+
exp_lon, exp_lat = truth_moving(start_lon, start_lat, pset.time[0])
304+
np.testing.assert_allclose(pset.lon, exp_lon, rtol=rtol)
305+
np.testing.assert_allclose(pset.lat, exp_lat, rtol=rtol)
306306
if method == "RK4_3D":
307-
np.testing.assert_allclose(pset.depth_nextloop, exp_lat, rtol=rtol)
307+
np.testing.assert_allclose(pset.depth, exp_lat, rtol=rtol)
308308

309309

310310
@pytest.mark.parametrize(
@@ -347,9 +347,9 @@ def truth_moving(x_0, y_0, t):
347347
)
348348
return lon, lat
349349

350-
exp_lon, exp_lat = truth_moving(start_lon, start_lat, pset.time_nextloop[0])
351-
np.testing.assert_allclose(pset.lon_nextloop, exp_lon, rtol=rtol)
352-
np.testing.assert_allclose(pset.lat_nextloop, exp_lat, rtol=rtol)
350+
exp_lon, exp_lat = truth_moving(start_lon, start_lat, pset.time[0])
351+
np.testing.assert_allclose(pset.lon, exp_lon, rtol=rtol)
352+
np.testing.assert_allclose(pset.lat, exp_lat, rtol=rtol)
353353

354354

355355
@pytest.mark.parametrize(
@@ -469,7 +469,7 @@ def periodicBC(particles, fieldset): # pragma: no cover
469469

470470
pset = parcels.ParticleSet(fieldset, lon=lonp, lat=latp)
471471
pset.execute([AdvectionEE, periodicBC], runtime=runtime, dt=np.timedelta64(6, "h"))
472-
np.testing.assert_allclose(pset.lat_nextloop, latp, atol=1e-1)
472+
np.testing.assert_allclose(pset.lat, latp, atol=1e-1)
473473

474474

475475
@pytest.mark.parametrize("method", ["RK4", "RK4_3D"])

tests/v4/test_particleset.py

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -123,7 +123,7 @@ def Addlon(particles, fieldset): # pragma: no cover
123123
particles.dlon += particles.dt / np.timedelta64(1, "s")
124124

125125
pset.execute(Addlon, dt=np.timedelta64(2, "s"), runtime=np.timedelta64(8, "s"), verbose_progress=False)
126-
assert np.allclose([p.lon_nextloop for p in pset], [8 - t for t in times])
126+
assert np.allclose([p.lon + p.dlon for p in pset], [8 - t for t in times])
127127

128128

129129
def test_pset_add_explicit(fieldset):

0 commit comments

Comments
 (0)