[fluidnumericsJoe]

Description

This PR adds an example of limiting particle timestep by local cfl in schism tutorial

Checklist

• This PR targets the correct branch (main for normal development, v3-support for v3 support)

[erikvansebille]

Thanks for adding this extra context, @fluidnumericsJoe! Below a few small comments (I can't put comments directly on the notebook in GitHub)

1. Is "altitude" the right word in the comment of cell [5]? That to me suggests (3D) height?
2. Avoid variable names like h, e and p in that same cell - make them more explicit so users understand what's going on?
3. CFL number hasn't been introduced in cell [5]. Could you move that whole cell down to the section of cell [9]? That will also make it easier to refer to this part of the tutorial directly, because all information is in one subsection
4. Make explicit that the math symbol $C$ is the same as CFL_target
5. Should we add a link to the CFL-section on the tutorial_dt_integrators too?

[fluidnumericsJoe]

1. "Altitude": reworded. The text now describes the cell width as the shortest distance from a corner of the triangle to the opposite edge (twice the area divided by the longest edge), and notes that this is one of several valid measures.

2. Variable names: p, e, h are now corners (points), edges, edge_lengths, cell_area and cell_width, both in that cell and in the kernel and release cells.

3. Moving the cell: the cell-width code now lives in the "Limit the timestep by the local CFL number" section, together with the kernel. It opens by introducing the CFL number, $C = |u|,\Delta t/\Delta x$. To make the move possible, cell_width is no longer part of the velocity UxDataset. Instead it's built as its own small FieldSet and added with fieldset = fieldset + .... this doesn't seem to impact runtime for this example

4. $C$ vs CFL_target: now stated explicitly: "with $C$ = CFL_target and $\Delta x$ = cell_width". There's also a sentence on how fieldset.add_context makes CFL_target available in the kernel.

5. Link to tutorial_dt_integrators: added. I linked to "A priori estimation of dt", which makes the same $U,\Delta t < \Delta x$ argument for a structured grid.

[erikvansebille]

Thanks for making these updates, however I don't see the link back from the dt_integration notebook to this notebook?

[fluidnumericsJoe]

however I don't see the link back from the dt_integration notebook to this notebook?

It's in the section titled "Limit the timestep by the local CFL number" under the text "See a priori estimation of dt for the same reasoning applied to a structured grid."

It's also in the last cell titled "Adapting this tutorial to your own SCHISM setup" under the text " choosing a timestep and integration method"

[erikvansebille]

Sorry, I meant the other way around. So a link from the https://docs.parcels-code.org/en/latest/user_guide/examples/tutorial_dt_integrators.html notebook to this new notebook

[erikvansebille]

Sorry, I meant the other way around. So a link from the https://docs.parcels-code.org/en/latest/user_guide/examples/tutorial_dt_integrators.html notebook to this new notebook

I have now added such a link in e88f57d - as I was working on merging your code with the new main

[erikvansebille]

I now realise that the CFLLimitedDt Kernel in this PR expects that particles.dt is reset at the end of each loop - but this is not true anymore since #3019. I think that the implementation in the notebook now only decreases particles.dt, but never increases them - leading to a timeout in runtime

def CFLLimitedDt(particles, fieldset):
    """Shrink each particle's dt so it moves at most CFL_target cell widths per step."""
    cell_width = fieldset.cell_width[particles]
    u, v = fieldset.UV[particles]
    dt_cfl = (
        fieldset.CFL_target * cell_width / np.maximum(np.sqrt(u * u + v * v), 1e-12)
    )
    particles.dt = np.sign(particles.dt) * np.minimum(np.abs(particles.dt), dt_cfl)

Should the np.minimum() be removed in that last line? So that it is always simply dt_cfl? If I do that, I get an OutsideTimeInterval though...

@fluidnumericsJoe, any thoughts?

[erikvansebille]

I got the "old" behaviour back by using the particles.next_dt variable, which is also used for RK45

diff --git a/docs/user_guide/convert_schism.ipynb b/docs/user_guide/convert_schism.ipynb
index 5692fb20..661a7c7a 100644
--- a/docs/user_guide/convert_schism.ipynb
+++ b/docs/user_guide/convert_schism.ipynb
@@ -271,7 +271,10 @@
     "SchismParticle = parcels.Particle.add_variable(\n",
-    "    parcels.Variable(\"out_of_bounds\", dtype=np.int32, initial=0)\n",
+    "    [\n",
+    "        parcels.Variable(\"out_of_bounds\", dtype=np.int32, initial=0),\n",
+    "        parcels.Variable(\"next_dt\", dtype=np.float32, initial=0),\n",
+    "    ]\n",
     ")\n",
@@ -363,6 +366,7 @@
     "    dt_cfl = (\n",
     "        fieldset.CFL_target * cell_width / np.maximum(np.sqrt(u * u + v * v), 1e-12)\n",
     "    )\n",
+    "    particles.next_dt = particles.dt.copy()\n",
     "    particles.dt = np.sign(particles.dt) * np.minimum(np.abs(particles.dt), dt_cfl)\n",
diff --git a/src/parcels/_core/kernel.py b/src/parcels/_core/kernel.py
index 5d216749..ebd7c6e2 100644
--- a/src/parcels/_core/kernel.py
+++ b/src/parcels/_core/kernel.py
@@ -127,8 +127,8 @@ class Kernel:
-        if hasattr(self.fieldset, "RK45_tol"):
-            # Update dt in case it's increased in RK45 kernel
+        if "next_dt" in particles._data:
+            # Update dt in case it's set by the kernel
             particles.dt = particles.next_dt

This is of course not ideal, as it requires using the next_dt (and thus would not allow users to change the particles.dt in the kernel), but at least is shows that indeed the change in #3019 is now breaking this kernel