<commit_before><commit_msg>Add a tutorial for NestedField<commit_after>{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Tutorial on how to combine different Fields into a `NestedField` object"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In some applications, you may have access to different fields that cover each only part of the region of interest. Then, you would like to combine them all together. You may also have a field covering the entire region and another one only covering part of it, but with a higher resolution. The set of those fields form what we call nested fields\n",
    "\n",
    "It is possible to combine all those fields with kernels, either with different if/else statements depending on particle position, or using recovery kernels (if only two levels of nested fields).\n",
    "\n",
    "In Parcels, it is possible to combine all those fields into one `NestedField` object, who will try to interpolate succesively the different fields. The principle is the following:\n",
    "\n",
    "Each one it its turn, a field is interpolated: if the interpolation succeeds or if an error other than `ErrorOutOfBounds` is thrown, the function is stopped. Otherwise, next field is interpolated. NestedField returns an `ErrorOutOfBounds` only if last field is as well out of boundaries.\n",
    "\n",
    "This means that the order of the fields in the `NestedField` matters.\n",
    "\n",
    "This tutorial shows how to use these `NestedField` with a very idealised example."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "from parcels import Field, NestedField, FieldSet, ParticleSet, JITParticle, plotTrajectoriesFile, AdvectionRK4\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's first define a zonal and meridional velocity field defined on a high resolution (dx = 100m) 2kmx2km grid with a flat mesh. The zonal velocity is uniform and 1 m/s, and the meridional velocity is equal to .5 * cos(lon / 200 * pi / 2)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "dim = 21\n",
    "lon = np.linspace(0., 2e3, dim, dtype=np.float32)\n",
    "lat = np.linspace(0., 2e3, dim, dtype=np.float32)\n",
    "lon_g, lat_g = np.meshgrid(lon, lat)\n",
    "V1_data = .5*np.cos(lon_g / 200 * np.pi/2)\n",
    "U1 = Field('U1', np.ones((dim, dim), dtype=np.float32),\n",
    "          lon=lon, lat=lat)\n",
    "V1 = Field('V1', V1_data, grid=U1.grid)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's now define the same velocity field on a low resolution (dx = 2km) 20kmx20km grid."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "dim = 11\n",
    "lon = np.linspace(0., 20e3, dim, dtype=np.float32)\n",
    "lat = np.linspace(0., 20e3, dim, dtype=np.float32)\n",
    "lon_g, lat_g = np.meshgrid(lon, lat)\n",
    "V2_data = .5*np.cos(lon_g / 200 * np.pi/2)\n",
    "U2 = Field('U2', np.ones((dim, dim), dtype=np.float32),\n",
    "          lon=lon, lat=lat)\n",
    "V2 = Field('V2', V2_data, grid=U2.grid)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We now combine those fields into a `NestedField` and create the fieldset"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "U = NestedField('U', [U1, U2])\n",
    "V = NestedField('V', [V1, V2])\n",
    "fieldset = FieldSet(U, V)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "INFO: Compiled JITParticleAdvectionRK4 ==> /var/folders/h0/01fvrmn11qb62yjw7v1kn62r0000gq/T/parcels-503/3b9fbe7fb6d38c4c87ac6f6017636211.so\n"
     ]
    },
    {
     "data": {
      "image/png": "<KEY>",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10cd8d850>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "pset = ParticleSet(fieldset, pclass=JITParticle, lon=[0], lat=[500])\n",
    "output_file = pset.ParticleFile(name='NestedFieldParticle.nc', outputdt=50)\n",
    "pset.execute(AdvectionRK4, runtime=10000, dt=10, output_file=output_file)\n",
    "plotTrajectoriesFile('NestedFieldParticle.nc');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As we observe, there is a change of dynamic at lon=2000, which corresponds to the change of grid.\n",
    "\n",
    "The analytical solution to the problem:\n",
    "\\begin{align}\n",
    "dx/dt &= 1\\\\\n",
    "dy/dt &= .5 \\cos(x \\pi/400)\\\\\n",
    "\\text{with  } x(0) &= 0, y(0) = 500\n",
    "\\end{align}\n",
    "\n",
    "is\n",
    "\\begin{align}\n",
    "x(t) &= t\\\\\n",
    "y(t) &= 500 + 200/\\pi \\sin(t \\pi / 400)\n",
    "\\end{align}\n",
    "which is captured by the HR field but not the LR one."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Keep track of the field interpolated\n",
    "\n",
    "For different reasons, you may want to keep track of the field you have interpolated. You can do that easily by creating another field that share the grid with original fields.\n",
    "Watch out that this operation has a cost of a full interpolation operation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "F1 = Field('F1', np.ones((U1.grid.ydim, U1.grid.xdim), dtype=np.float32), grid=U1.grid)\n",
    "F2 = Field('F2', 2*np.ones((U2.grid.ydim, U2.grid.xdim), dtype=np.float32), grid=U2.grid)\n",
    "F = NestedField('F', [F1, F2])\n",
    "fieldset.add_field(F)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "INFO: Compiled SampleParticleSampleNestedFieldIndex ==> /var/folders/h0/01fvrmn11qb62yjw7v1kn62r0000gq/T/parcels-503/f11a0f2e8057f47494303cfd2166c89c.so\n",
      "WARNING: dt or runtime are zero, or endtime is equal to Particle.time. The kernels will be executed once, without incrementing time\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Particle (0, 500) interpolates Field #1\n",
      "Particle (10000, 500) interpolates Field #2\n"
     ]
    }
   ],
   "source": [
    "from parcels import Variable\n",
    "\n",
    "def SampleNestedFieldIndex(particle, fieldset, time, dt):\n",
    "    particle.f = fieldset.F[time, particle.lon, particle.lat, particle.depth]\n",
    "\n",
    "class SampleParticle(JITParticle):\n",
    "    f = Variable('f', dtype=np.int32)\n",
    "    \n",
    "pset = ParticleSet(fieldset, pclass= SampleParticle, lon=[0], lat=[500])\n",
    "pset.execute(SampleNestedFieldIndex, runtime=0, dt=0)\n",
    "print('Particle (%g, %g) interpolates Field #%d' % (pset[0].lon, pset[0].lat, pset[0].f))\n",
    "\n",
    "pset[0].lon = 10000\n",
    "pset.execute(SampleNestedFieldIndex, runtime=0, dt=0)\n",
    "print('Particle (%g, %g) interpolates Field #%d' % (pset[0].lon, pset[0].lat, pset[0].f))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.15"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
{'fraction_non_alphanumeric': 0.18148330058939097, 'fraction_numerical': 0.033030451866404714, 'mean_word_length': 2.978993649242794, 'pattern_counts': {'":': 99, '<': 5, '<?xml version=': 0, '>': 7, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '6071952', 'n_tokens_mistral': 3190, 'n_tokens_neox': 2999, 'n_words': 759}