Skip to content
Merged
Show file tree
Hide file tree
Changes from 3 commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
104 changes: 104 additions & 0 deletions docs/user_guide/examples/tutorial_statuscodes.md
Original file line number Diff line number Diff line change
@@ -0,0 +1,104 @@
---
file_format: mystnb
kernelspec:
name: python3
---

# Working with Status Codes

In order to capture errors in the [Kernel loop](explanation_kernelloop.md), Parcels uses a Status Code system. There are several Status Codes, listed below.

```{code-cell}
import parcels

for statuscode, val in parcels.StatusCode.__dict__.items():
if statuscode.startswith("__"):
continue
print(f"{statuscode} = {val}")
```

Once an error is thrown (for example, a Field Interpolation error), then the `particles.state` is updated to the corresponding status code. This gives you the flexibility to write a Kernel that checks for a status code and does something with it.

For example, you can write a Kernel that checks for `particles.state == parcels.StatusCode.ErrorOutOfBounds` and deletes the particle, and then append this custom Kernel to the Kernel list in `pset.execute()`.

```
def DeleteOutOfBounds(particles, fieldset):
out_of_bounds = particles.state == parcels.StatusCode.ErrorOutOfBounds
particles[out_of_bounds].state = parcels.StatusCode.Delete


def DeleteAnyError(particles, fieldset):
any_error = particles.state >= 50 # This captures all Errors
particles[any_error].state = parcels.StatusCode.Delete
```

But of course, you can also write code for more sophisticated behaviour than just deleting the particle. It's up to you! Note that if you don't delete the particle, you will have to update the `particles.state = parcels.StatusCode.Evaluate` yourself. For example:

```
def Move1DegreeWest(particles, fieldset):
out_of_bounds = particles.state == parcels.StatusCode.ErrorOutOfBounds
particles[out_of_bounds].dlon -= 1.0
particles[out_of_bounds].state = parcels.StatusCode.Evaluate
```

Or, if you want to make sure that particles don't escape through the water surface
Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
Or, if you want to make sure that particles don't escape through the water surface
When `particles.dz` is too large, they will be displaced through the surface, throwing an `ErrorThroughSurface`. We can adapt the particles.dz of these particles to move no further than the surface:


```{code-cell}
def KeepInOcean(particles, fieldset):
# find particles that move through the surface
through_surface = particles.state == parcels.StatusCode.ErrorThroughSurface

# move particles to surface
Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
# move particles to surface
# move particles from current z exactly to the surface

particles[through_surface].dz = fieldset.W.grid.depth[0] - particles[through_surface].z
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This line may be a bit magical to new users. Does this need a bit of information why this works?

Or even replace with

Suggested change
particles[through_surface].dz = fieldset.W.grid.depth[0] - particles[through_surface].z
particles[through_surface].z = fieldset.W.grid.depth[0]
particles[through_surface].dz = 0

even though we advice not to change particles.z directly?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I have expanded both the description before the code-cell and the comment slightly. I personally like adapting dz because I might forget the second line in your suggestion otherwise and would worry others do the same, but it does make things cleaner.


# change state from error to evaluate
particles[through_surface].state = parcels.StatusCode.Evaluate
```

Kernel functions such as the ones above can then be added to the list of kernels in `pset.execute()`.

Let's add the `KeepInOcean` Kernel to an particle simulation where particles move through the surface:

```{code-cell}
import numpy as np
from parcels._datasets.structured.generated import simple_UV_dataset

ds = simple_UV_dataset(dims=(1, 2, 5, 4), mesh="flat").isel(time=0)

dx, dy = 1.0 / len(ds.XG), 1.0 / len(ds.YG)

# Add W velocity that pushes through surface
ds["W"] = ds["U"] - 0.1 # 0.1 m/s towards the surface

grid = parcels.XGrid.from_dataset(ds, mesh="flat")
U = parcels.Field("U", ds["U"], grid, interp_method=parcels.interpolators.XLinear)
V = parcels.Field("V", ds["V"], grid, interp_method=parcels.interpolators.XLinear)
W = parcels.Field("W", ds["W"], grid, interp_method=parcels.interpolators.XLinear)
UVW = parcels.VectorField("UVW", U, V, W)
fieldset = parcels.FieldSet([U, V, W, UVW])
```

If we advect particles with the `AdvectionRK2_3D` kernel, Parcels will raise a `FieldOutOfBoundSurfaceError`:

```{code-cell}
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Will it be a problem for the GitHub Action that this cell fails?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The raises-exception tag in the following line is recognized by ReadTheDocs, and all checks pass, so I think it’s fine

:tags: [raises-exception]
pset = parcels.ParticleSet(fieldset, parcels.Particle, z=[0.5], lat=[2], lon=[1.5])
kernels = [parcels.kernels.AdvectionRK2_3D]
pset.execute(kernels, runtime=np.timedelta64(1, "m"), dt=np.timedelta64(1, "s"), verbose_progress=False)
```

When we add the `KeepInOcean` Kernel, particles will stay at the surface:

```{code-cell}
pset = parcels.ParticleSet(fieldset, parcels.Particle, z=[0.5], lat=[2], lon=[1.5])

kernels = [parcels.kernels.AdvectionRK2_3D, KeepInOcean]

pset.execute(kernels,runtime=np.timedelta64(20, "s"), dt=np.timedelta64(1, "s"), verbose_progress=False)

print(f"particle z at end of run = {pset.z}")
```

```{note}
Kernels that control what to do with `particles.state` should typically be added at the _end_ of the Kernel list, because otherwise later Kernels may overwrite the `particles.state` or the `particles.dlon` variables (see [Kernel loop explanation](explanation_kernelloop.md)).
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The link doesn't work; tries to go to #explanation_kernelloop on the page itself?

```
1 change: 1 addition & 0 deletions docs/user_guide/index.md
Original file line number Diff line number Diff line change
Expand Up @@ -45,6 +45,7 @@ examples/tutorial_delaystart.ipynb
:titlesonly:

examples/tutorial_sampling.ipynb
examples/tutorial_statuscodes.ipynb
examples/tutorial_gsw_density.ipynb
examples/tutorial_Argofloats.ipynb
```
Expand Down