Skip to content

Commit 00c4447

Browse files
reint-fischerreint-fischer
reint-fischer
authored and
reint-fischer
committed
add review suggestions
1 parent 89073c7 commit 00c4447

File tree

1 file changed

+118
-58
lines changed

1 file changed

+118
-58
lines changed

‎docs/user_guide/examples/tutorial_dt_integrators.ipynb‎

Lines changed: 118 additions & 58 deletions
Original file line numberDiff line numberDiff line change
@@ -29,7 +29,7 @@
2929
"id": "2",
3030
"metadata": {},
3131
"source": [
32-
"### A priori estimation\n",
32+
"### A priori estimation of `dt`\n",
3333
"\n",
3434
"In this example, we will estimate the appropriate timestep to compute advection of a virtual particle through an ocean current field. \n",
3535
"\n",
@@ -43,7 +43,11 @@
4343
"\n",
4444
"where $\\mathbf{v}(\\mathbf{x},t) = (u(\\mathbf{x},t), v(\\mathbf{x},t))$ describes the ocean velocity field at position $\\mathbf{x}$ at time $t$.\n",
4545
"\n",
46-
"To estimate the timescale that we want to resolve, we can think about the scales at which advection varies. Here we use the daily velocity fields at 1/12th degree horizontal resolution from the Copernicus Marine Service. This means that the velocity will vary in time at scales >= 24 hours, and in space at scales >= 1/12th degree."
46+
"To estimate the timescale that we want to resolve, we can think about the scales at which advection varies. Here we use the daily velocity fields at 1/12th degree horizontal resolution from the Copernicus Marine Service. This means that the velocity will vary in time at scales >= 24 hours, and in space at scales >= 1/12th degree.\n",
47+
"\n",
48+
"```{note}\n",
49+
"Our displacement occurs in units of longitude and latitde, but our velocity field is in m/s. Read [this guide](./tutorial_unitconverters.ipynb) to understand how Parcels converts these units under the hood.\n",
50+
"```"
4751
]
4852
},
4953
{
@@ -55,6 +59,7 @@
5559
"source": [
5660
"import matplotlib.pyplot as plt\n",
5761
"import numpy as np\n",
62+
"import pandas as pd\n",
5863
"import xarray as xr\n",
5964
"\n",
6065
"import parcels\n",
@@ -143,10 +148,23 @@
143148
"id": "7",
144149
"metadata": {},
145150
"source": [
151+
"$$\n",
152+
"\\begin{aligned}\n",
153+
"\\text{d}t < \\frac{1}{12 U_{max}}\n",
154+
"\\end{aligned}\n",
155+
"$$\n",
156+
"\n",
146157
"Using `U_max_surface_deg`, we find a second estimated limit of an appropriate `dt`:\n",
147158
"\n",
148159
"$$\n",
149160
"\\begin{aligned}\n",
161+
"\\text{d}t < \\frac{1}{12 * 1.71e-5} = 4.9e3 \\text{ seconds}\n",
162+
"\\end{aligned}\n",
163+
"$$\n",
164+
"\n",
165+
"\n",
166+
"$$\n",
167+
"\\begin{aligned}\n",
150168
"\\text{d}t < 2 \\text{ hours}\n",
151169
"\\end{aligned}\n",
152170
"$$"
@@ -184,18 +202,18 @@
184202
"id": "10",
185203
"metadata": {},
186204
"source": [
187-
"We simulate particles using a range of timesteps that differ by a factor 2-10, starting at (dt < 24 hours). We also keep track of the time it takes to run each simulation:"
205+
"We simulate particles using a range of timesteps that differ by a factor 2-10, starting at (dt < 24 hours). We also keep track of the time it takes to run each simulation:\n",
206+
"\n",
207+
"```{warning}\n",
208+
"`dt` must be chosen such that $k`dt`$ fits exactly in the FieldSet dt (24 hours), where $k$ is an integer. \n",
209+
"```"
188210
]
189211
},
190212
{
191213
"cell_type": "code",
192214
"execution_count": null,
193215
"id": "11",
194-
"metadata": {
195-
"tags": [
196-
"hide-output"
197-
]
198-
},
216+
"metadata": {},
199217
"outputs": [],
200218
"source": [
201219
"import time\n",
@@ -235,10 +253,15 @@
235253
" # time and run simulation\n",
236254
" start_time = time.time()\n",
237255
" pset.execute(\n",
238-
" parcels.kernels.AdvectionRK2, runtime=runtime, dt=dt, output_file=pfile\n",
256+
" parcels.kernels.AdvectionRK2,\n",
257+
" runtime=runtime,\n",
258+
" dt=dt,\n",
259+
" output_file=pfile,\n",
260+
" verbose_progress=False,\n",
239261
" )\n",
240262
" sim_duration_i = time.time() - start_time\n",
241-
" sim_duration[i] = sim_duration_i"
263+
" sim_duration[i] = sim_duration_i\n",
264+
" print(f\"Simulation duration = {np.round(sim_duration_i, 2)} seconds\")"
242265
]
243266
},
244267
{
@@ -443,32 +466,40 @@
443466
"ax[1].set_xlabel(\"dt (minutes)\")\n",
444467
"ax[1].grid()\n",
445468
"ax[1].legend()\n",
446-
"plt.show()\n",
447-
"print(f\"dt choices = {(dt_choices / np.timedelta64(1, 'm')).astype(int)} minutes\")\n",
448-
"print(f\"precision = {np.round(np.mean(dist_end, axis=1), 2)} km\")\n",
449-
"print(f\"sim duration = {np.round(sim_duration, 2)} s\")"
469+
"plt.show()"
450470
]
451471
},
452472
{
453473
"cell_type": "markdown",
454474
"id": "20",
455475
"metadata": {},
456476
"source": [
457-
"We can see that in our simulation advecting particles for 7 days, the effect of `dt` on the precision of our simulation is approximately linear. The precision for a simulation with a timestep of 20 minutes is order of magnitude ~100 m. The effect on the time it takes to run a simulation is not linear in our case however; it increases sharply as we decrease our timestep. This may be optimized using more efficient chunking.\n",
458-
"\n",
459-
"| `dt` | Mean separation distance after 7 days (km) | Simulation duration (s) | \n",
460-
"| ---------- | ------------------------------------------ | ------------------------ |\n",
461-
"| 12 hours | 19.52 | 0.11 |\n",
462-
"| 6 hours | 9.64 | 0.37 | \n",
463-
"| 1 hour | 1.4 | 0.94 | \n",
464-
"| 20 minutes | 0.38 | 3.19 | \n",
465-
"| 5 minutes | x | 12.34 | \n"
477+
"We can see that in our simulation advecting particles for 7 days, the effect of `dt` on the precision of our simulation is approximately linear. The precision for a simulation with a timestep of 20 minutes is order of magnitude ~100 m. The effect on the time it takes to run a simulation is not linear in our case however; it increases sharply as we decrease our timestep. This may be optimized using more efficient chunking."
466478
]
467479
},
468480
{
469-
"cell_type": "markdown",
481+
"cell_type": "code",
482+
"execution_count": null,
470483
"id": "21",
471484
"metadata": {},
485+
"outputs": [],
486+
"source": [
487+
"df = pd.DataFrame(\n",
488+
" {\n",
489+
" \"dt\": dt_choices,\n",
490+
" \"Mean separation distance (km)\": np.append(\n",
491+
" np.round(np.mean(dist_end, axis=1), 2), [None]\n",
492+
" ),\n",
493+
" \"Simulation duration (s)\": np.round(sim_duration, 2),\n",
494+
" }\n",
495+
")\n",
496+
"df"
497+
]
498+
},
499+
{
500+
"cell_type": "markdown",
501+
"id": "22",
502+
"metadata": {},
472503
"source": [
473504
"```{note}\n",
474505
"The desired precision is not always best measured by the separation distance of individual trajectories. Depending on the application of your Parcels simulation and the process you are computing, other metrics may be more suitable.\n",
@@ -477,7 +508,7 @@
477508
},
478509
{
479510
"cell_type": "markdown",
480-
"id": "22",
511+
"id": "23",
481512
"metadata": {},
482513
"source": [
483514
"## Integration schemes\n",
@@ -488,7 +519,7 @@
488519
{
489520
"cell_type": "code",
490521
"execution_count": null,
491-
"id": "23",
522+
"id": "24",
492523
"metadata": {},
493524
"outputs": [],
494525
"source": [
@@ -501,7 +532,7 @@
501532
{
502533
"cell_type": "code",
503534
"execution_count": null,
504-
"id": "24",
535+
"id": "25",
505536
"metadata": {},
506537
"outputs": [],
507538
"source": [
@@ -512,7 +543,7 @@
512543
{
513544
"cell_type": "code",
514545
"execution_count": null,
515-
"id": "25",
546+
"id": "26",
516547
"metadata": {},
517548
"outputs": [],
518549
"source": [
@@ -522,7 +553,7 @@
522553
},
523554
{
524555
"cell_type": "markdown",
525-
"id": "26",
556+
"id": "27",
526557
"metadata": {},
527558
"source": [
528559
"The higher-order methods use weighted intermediate steps in time and space to obtain a more accurate estimate of `dlat` and `dlon` for a given timestep.\n",
@@ -533,7 +564,7 @@
533564
{
534565
"cell_type": "code",
535566
"execution_count": null,
536-
"id": "27",
567+
"id": "28",
537568
"metadata": {},
538569
"outputs": [],
539570
"source": [
@@ -547,7 +578,7 @@
547578
{
548579
"cell_type": "code",
549580
"execution_count": null,
550-
"id": "28",
581+
"id": "29",
551582
"metadata": {
552583
"tags": [
553584
"hide-output"
@@ -582,48 +613,59 @@
582613
" f\"Begin simulation for (scheme: {advection_scheme.__name__}, dt={int(dt / np.timedelta64(1, 's'))} s)\"\n",
583614
" )\n",
584615
" start_time = time.time()\n",
585-
" pset.execute(advection_scheme, runtime=runtime, dt=dt, output_file=pfile)\n",
616+
" pset.execute(\n",
617+
" advection_scheme,\n",
618+
" runtime=runtime,\n",
619+
" dt=dt,\n",
620+
" output_file=pfile,\n",
621+
" verbose_progress=False,\n",
622+
" )\n",
586623
" sim_duration_ij = time.time() - start_time\n",
587-
" sim_duration[i, j] = sim_duration_ij"
624+
" sim_duration[i, j] = sim_duration_ij\n",
625+
" print(f\"Simulation duration = {np.round(sim_duration_ij, 2)} seconds\")"
588626
]
589627
},
590628
{
591629
"cell_type": "code",
592630
"execution_count": null,
593-
"id": "29",
631+
"id": "30",
594632
"metadata": {},
595633
"outputs": [],
596634
"source": [
597635
"scheme_colours = np.linspace(0, 1, len(advection_schemes), endpoint=True)\n",
598636
"# Now let's compare different advection schemes with the same timestep\n",
599-
"fig, axs = plt.subplots(nrows=1, ncols=len(dt_choices), figsize=(20, 5))\n",
637+
"fig, axs = plt.subplots(nrows=2, ncols=3, figsize=(15, 10))\n",
600638
"for i, dt in enumerate(dt_choices):\n",
601-
" axs[i].set_title(f\"{str(dt)}\")\n",
602-
" axs[i].set_xlabel(\"Longitude\")\n",
639+
" m = i // 3\n",
640+
" n = i % 3\n",
641+
" axs[m, n].set_title(f\"dt = {str(dt)}\")\n",
642+
" axs[m, n].set_xlabel(\"Longitude\")\n",
603643
" for j, advection_scheme in enumerate(advection_schemes):\n",
604644
" ds = xr.open_zarr(\n",
605645
" f\"output/{advection_scheme.__name__}_dt_{int(dt / np.timedelta64(1, 's'))}s.zarr\"\n",
606646
" )\n",
607647
" labels = [f\"{advection_scheme.__name__}\"] + [None] * (ds.lon.shape[0] - 1)\n",
608-
" axs[i].plot(\n",
648+
" axs[m, n].plot(\n",
609649
" ds.lon.T,\n",
610650
" ds.lat.T,\n",
611651
" alpha=0.75,\n",
612652
" color=plt.cm.viridis(scheme_colours[j]),\n",
613653
" label=labels,\n",
614654
" )\n",
615-
" axs[i].scatter(\n",
655+
" axs[m, n].scatter(\n",
616656
" ds.lon[:, 0], ds.lat[:, 0], c=\"r\", marker=\"s\", label=\"starting locations\"\n",
617657
" )\n",
618-
" axs[i].grid()\n",
619-
" axs[0].legend()\n",
620-
" axs[0].set_ylabel(\"Latitude\")\n",
658+
" axs[m, n].grid()\n",
659+
"axs[-1, -1].axis(\"off\")\n",
660+
"axs[0, 0].legend()\n",
661+
"axs[0, 0].set_ylabel(\"Latitude\")\n",
662+
"axs[1, 0].set_ylabel(\"Latitude\")\n",
621663
"plt.show()"
622664
]
623665
},
624666
{
625667
"cell_type": "markdown",
626-
"id": "30",
668+
"id": "31",
627669
"metadata": {},
628670
"source": [
629671
"Clearly, for longer timesteps, the RK2 and RK4 schemes perform better. However, if the timestep is appropriate, as we have determined in the previous section, then the Explicit Euler scheme does not perform notably different."
@@ -632,19 +674,22 @@
632674
{
633675
"cell_type": "code",
634676
"execution_count": null,
635-
"id": "31",
677+
"id": "32",
636678
"metadata": {},
637679
"outputs": [],
638680
"source": [
639681
"dist_end = np.zeros((len(advection_schemes) - 1, len(dt_choices), npart))\n",
640682
"\n",
641-
"fig, axs = plt.subplots(nrows=1, ncols=len(dt_choices), figsize=(20, 5))\n",
683+
"fig, axs = plt.subplots(nrows=2, ncols=3, figsize=(15, 10))\n",
684+
"\n",
642685
"for i, dt in enumerate(dt_choices):\n",
643-
" axs[i].set_title(f\"dt = {str(dt)}\")\n",
644-
" axs[i].set_xlabel(\"Time\")\n",
645-
" axs[i].tick_params(\"x\", rotation=45)\n",
646-
" axs[i].set_yscale(\"log\")\n",
647-
" axs[i].set_ylim(1e-4, 1e1)\n",
686+
" m = i // 3\n",
687+
" n = i % 3\n",
688+
" axs[m, n].set_title(f\"dt = {str(dt)}\")\n",
689+
" axs[m, n].set_xlabel(\"Time\")\n",
690+
" axs[m, n].tick_params(\"x\", rotation=45)\n",
691+
" axs[m, n].set_yscale(\"log\")\n",
692+
" axs[m, n].set_ylim(1e-4, 1e1)\n",
648693
" ds_RK4 = xr.open_zarr(\n",
649694
" f\"output/AdvectionRK4_dt_{int(dt / np.timedelta64(1, 's'))}s.zarr\"\n",
650695
" )\n",
@@ -667,23 +712,26 @@
667712
" lat_valid = ds.lat.where(~np.isnan(ds.lat).compute(), drop=True).values\n",
668713
" dist = dist_km(lon_valid, lon_valid_RK4, lat_valid, lat_valid_RK4)\n",
669714
" time_valid = ds.time.where(~np.isnan(ds.time).compute(), drop=True).values\n",
670-
" axs[i].plot(\n",
715+
" axs[m, n].plot(\n",
671716
" time_valid.T,\n",
672717
" dist.T,\n",
673718
" alpha=0.75,\n",
674719
" color=plt.cm.viridis(scheme_colours[j]),\n",
675720
" label=labels,\n",
676721
" )\n",
677722
" dist_end[j, i] = dist[:, -1]\n",
678-
" axs[i].grid()\n",
679-
" axs[0].legend()\n",
680-
" axs[0].set_ylabel(\"Distance (km)\")\n",
723+
" axs[m, n].grid()\n",
724+
"axs[-1, -1].axis(\"off\")\n",
725+
"axs[0, 0].legend()\n",
726+
"axs[0, 0].set_ylabel(\"Latitude\")\n",
727+
"axs[1, 0].set_ylabel(\"Latitude\")\n",
728+
"plt.tight_layout()\n",
681729
"plt.show()"
682730
]
683731
},
684732
{
685733
"cell_type": "markdown",
686-
"id": "32",
734+
"id": "33",
687735
"metadata": {},
688736
"source": [
689737
"By quantifying the precision of the integration methods, we can see that for a given timestep the Runge-Kutta methods perform orders of magnitude better than the Explicit Euler method. In this example, the error associated with the selected integration methods is smaller than that of the range of timesteps."
@@ -692,7 +740,7 @@
692740
{
693741
"cell_type": "code",
694742
"execution_count": null,
695-
"id": "33",
743+
"id": "34",
696744
"metadata": {},
697745
"outputs": [],
698746
"source": [
@@ -734,15 +782,15 @@
734782
},
735783
{
736784
"cell_type": "markdown",
737-
"id": "34",
785+
"id": "35",
738786
"metadata": {},
739787
"source": [
740788
"In this last figure, we see that the improvement of RK2 by orders of magnitude with respect to EE comes at a small computational cost. Since the RK2 and RK4 methods are practically indistinguishable, using the RK2 method with a timestep of 1 hour or 20 minutes, depending on the application, would be an appropriate choice for this simulation."
741789
]
742790
},
743791
{
744792
"cell_type": "markdown",
745-
"id": "35",
793+
"id": "36",
746794
"metadata": {},
747795
"source": [
748796
"### Flow conditions\n",
@@ -762,6 +810,18 @@
762810
"Check out [this notebook](https://github.com/Parcels-code/10year-anniversary-session2/blob/main/solutions/lorenz_and_lotka_volterra_solutions.ipynb) to learn how to model the Lorenz equations with Parcels!\n",
763811
"```"
764812
]
813+
},
814+
{
815+
"cell_type": "markdown",
816+
"id": "37",
817+
"metadata": {},
818+
"source": [
819+
"## Summary\n",
820+
"If you want to test the accuracy and efficiency of your own simulation, here is a brief list of things to consider:\n",
821+
"- **Use the temporal and spatial resolution of the input to estimate the timescales you need to resolve.** `dt` must fit into the FieldSet dt by a factor of an integer.\n",
822+
"- **Run the simulation for the full runtime.** Many inaccuracies grow over time so we must \n",
823+
"- **Consider which step introduces the largest error.** If we compute advection and diffusion, increasing the accuracy of one process much more than the other will not improve the results."
824+
]
765825
}
766826
],
767827
"metadata": {

0 commit comments

Comments
 (0)