parcel sampling debug

parcel_interp_debug.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# debug issue with interpolation with parcels simulations\n",
    "\n",
    "See [GH issue](https://github.com/OceanParcels/parcels/issues/1122)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "INFO: Compiled ParcelsRandom ==> /var/folders/pb/vp119fcj1pn3dc3x2dz4dkt80000gn/T/parcels-501/libparcels_random_bfa75db7-25e5-46d7-afd4-6a657983c50f.so\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "import xarray as xr\n",
    "from parcels import FieldSet, ParticleSet, JITParticle, Variable, AdvectionRK4, ParticleFile"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "class SampleParticle(JITParticle):\n",
    "    p = Variable('p', dtype=np.float32)\n",
    "\n",
    "def SampleP(particle, fieldset, time):\n",
    "    particle.p = fieldset.P[time, particle.depth, particle.lat, particle.lon]\n",
    "    \n",
    "def SamplePoff(particle, fieldset, time):\n",
    "    \"\"\" offset sampling by dt\n",
    "    \"\"\"\n",
    "    particle.p = fieldset.P[time+particle.dt, particle.depth, particle.lat, particle.lon]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### generate input fields\n",
    "\n",
    "\n",
    "\\begin{align}\n",
    "x &= [0,1], \\\\\n",
    "y &= [0,1], \\\\\n",
    "u(x,y,t) &= 0.5, \\\\\n",
    "v(x,y,t) &= 0, \\\\\n",
    "p(x,y,t) &= 10t + x, \\\\\n",
    "\\end{align}\n",
    "\n",
    "We'll follow a parcel that is initially at $x=0$ and $y=0$ with $p(x=0,y=0,t=0)=0$.\n",
    "\n",
    "After a unit time step, the parcel is at $x=0.5$ and $y=0$ with $p(x=0.5,y=0,t=1)=10.5$\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "WARNING: Casting time data to np.float64\n",
      "WARNING: Casting field data to np.float32\n"
     ]
    }
   ],
   "source": [
    "dims = [5, 4, 2]\n",
    "dx, dy = 1./dims[0], 1./dims[1]\n",
    "dimensions = {'lon': np.linspace(0., 1., dims[0], dtype=np.float32),\n",
    "              'lat': np.linspace(0., 1., dims[1], dtype=np.float32),\n",
    "              'time': np.arange(dims[2], dtype=np.float32),\n",
    "             }\n",
    "\n",
    "p = (10*np.arange(dims[2])[None, None, :]\n",
    "     + dimensions['lon'][:, None, None]\n",
    "     + np.zeros(dims[1], dtype=np.float32)[None, :, None]\n",
    "    )\n",
    "\n",
    "data = {'U': np.ones(dims, dtype=np.float32)/2,\n",
    "        'V': np.zeros(dims, dtype=np.float32),\n",
    "        'P': p   \n",
    "       }\n",
    "fieldset = FieldSet.from_data(data, dimensions, mesh='flat', transpose=True)\n",
    "\n",
    "xv, yv = np.meshgrid(np.arange(0, 1, 0.5), np.arange(0, 1, 0.5))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "### 1: Base configuration"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "INFO: Compiled SampleParticleSampleP ==> /var/folders/pb/vp119fcj1pn3dc3x2dz4dkt80000gn/T/parcels-501/0120dab411f0fda7d0c2708c638cd1f6_0.so\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "P[0](lon=0.000000, lat=0.000000, depth=0.000000, p=0.000000, time=1.000000)\n",
       "P[1](lon=0.500000, lat=0.000000, depth=0.000000, p=0.500000, time=1.000000)\n",
       "P[2](lon=0.000000, lat=0.500000, depth=0.000000, p=0.000000, time=1.000000)\n",
       "P[3](lon=0.500000, lat=0.500000, depth=0.000000, p=0.500000, time=1.000000)"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pset = ParticleSet(fieldset, pclass=SampleParticle, lon=xv.flatten(), lat=yv.flatten())\n",
    "pfile = ParticleFile(\"tmp1.nc\", pset, outputdt=1)\n",
    "pset.execute(SampleP, endtime=1, dt=1, output_file=pfile)\n",
    "pfile.close()\n",
    "pset"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "lon=0.0\n",
      "p=0.0\n",
      "time=1.0\n"
     ]
    }
   ],
   "source": [
    "def print_nc(n):\n",
    "    ds = xr.open_dataset(f\"tmp{n}.nc\").isel(obs=1)\n",
    "    print(\"lon={}\".format(ds[\"lon\"].values[0]))\n",
    "    print(\"p={}\".format(ds[\"p\"].values[0]))\n",
    "    print(\"time={}\".format(ds[\"time\"].values[0]/np.timedelta64(1,'s')))\n",
    "\n",
    "print_nc(1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Information stamped at $t$ is:\n",
    "\n",
    "- Position: $x(t-dt)$\n",
    "- Variable $p$ evaluated at $x(t-dt)$ and $t-dt$\n",
    "\n",
    "Not so bad of a situation as only the timestamp is incorrect."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "### 2: manually sum kernels with RK4+SampleP"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "INFO: Compiled SampleParticleAdvectionRK4SampleP ==> /var/folders/pb/vp119fcj1pn3dc3x2dz4dkt80000gn/T/parcels-501/a14b6229de31167b0795e58411e31226_0.so\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "P[4](lon=0.500000, lat=0.000000, depth=0.000000, p=0.500000, time=1.000000)\n",
       "P[5](lon=1.000000, lat=0.000000, depth=0.000000, p=1.000000, time=1.000000)\n",
       "P[6](lon=0.500000, lat=0.500000, depth=0.000000, p=0.500000, time=1.000000)\n",
       "P[7](lon=1.000000, lat=0.500000, depth=0.000000, p=1.000000, time=1.000000)"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pset = ParticleSet(fieldset, pclass=SampleParticle, lon=xv.flatten(), lat=yv.flatten())\n",
    "kernels = AdvectionRK4 + pset.Kernel(SampleP)\n",
    "pfile = ParticleFile(\"tmp2.nc\", pset, outputdt=1)\n",
    "pset.execute(kernels, endtime=1, dt=1, output_file=pfile)\n",
    "pfile.close()\n",
    "pset"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "lon=0.5\n",
      "p=0.5\n",
      "time=1.0\n"
     ]
    }
   ],
   "source": [
    "print_nc(2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Information stamped at $t$ is:\n",
    "\n",
    "- Position : $x(t)$\n",
    "- Variable $p$ evaluated at $x(t)$ and $t-dt$\n",
    "\n",
    "**!Danger!: the interpolation does not use consistent position and time**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "\n",
    "### 3: manually sum kernels with SampleP+RK4\n",
    "\n",
    "(personal note: used in production run)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "INFO: Compiled SampleParticleSamplePAdvectionRK4 ==> /var/folders/pb/vp119fcj1pn3dc3x2dz4dkt80000gn/T/parcels-501/89eacaae3f5adca1917d7934eeff2fb2_0.so\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "P[8](lon=0.500000, lat=0.000000, depth=0.000000, p=0.000000, time=1.000000)\n",
       "P[9](lon=1.000000, lat=0.000000, depth=0.000000, p=0.500000, time=1.000000)\n",
       "P[10](lon=0.500000, lat=0.500000, depth=0.000000, p=0.000000, time=1.000000)\n",
       "P[11](lon=1.000000, lat=0.500000, depth=0.000000, p=0.500000, time=1.000000)"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pset = ParticleSet(fieldset, pclass=SampleParticle, lon=xv.flatten(), lat=yv.flatten())\n",
    "kernels = pset.Kernel(SampleP) + AdvectionRK4\n",
    "pfile = ParticleFile(\"tmp3.nc\", pset, outputdt=1)\n",
    "pset.execute(kernels, endtime=1, dt=1, output_file=pfile)\n",
    "pfile.close()\n",
    "pset"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "lon=0.5\n",
      "p=0.0\n",
      "time=1.0\n"
     ]
    }
   ],
   "source": [
    "print_nc(3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Information stamped at $t$ is:\n",
    "\n",
    "- Position: $x(t)$\n",
    "- Variable $p$ evaluated at $x(t-dt)$ and $t-dt$\n",
    "\n",
    "Not so bad of a situation but need to keep in mind position and $p$ refer to differnt timestamps"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "\n",
    "### 4: manually sum kernels with RK4+SamplePoff"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "INFO: Compiled SampleParticleAdvectionRK4SamplePoff ==> /var/folders/pb/vp119fcj1pn3dc3x2dz4dkt80000gn/T/parcels-501/335ced1adbde87f6ebe3b6a81a495b85_0.so\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "P[12](lon=0.500000, lat=0.000000, depth=0.000000, p=10.500000, time=1.000000)\n",
       "P[13](lon=1.000000, lat=0.000000, depth=0.000000, p=11.000000, time=1.000000)\n",
       "P[14](lon=0.500000, lat=0.500000, depth=0.000000, p=10.500000, time=1.000000)\n",
       "P[15](lon=1.000000, lat=0.500000, depth=0.000000, p=11.000000, time=1.000000)"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pset = ParticleSet(fieldset, pclass=SampleParticle, lon=xv.flatten(), lat=yv.flatten())\n",
    "kernels = AdvectionRK4 + pset.Kernel(SamplePoff)\n",
    "pfile = ParticleFile(\"tmp4.nc\", pset, outputdt=1)\n",
    "pset.execute(kernels, endtime=1, dt=1, output_file=pfile)\n",
    "pfile.close()\n",
    "pset"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "lon=0.5\n",
      "p=10.5\n",
      "time=1.0\n"
     ]
    }
   ],
   "source": [
    "print_nc(4)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Information stamped at $t$ is:\n",
    "\n",
    "- Position: $x(t)$\n",
    "- Variable $p$ evaluated at $x(t)$ and $t$\n",
    "\n",
    "**This appear to be the correct solution**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "\n",
    "### 5: manually sum kernels with SamplePoff+RK4"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "INFO: Compiled SampleParticleSamplePoffAdvectionRK4 ==> /var/folders/pb/vp119fcj1pn3dc3x2dz4dkt80000gn/T/parcels-501/2782524329af42865d9a1753a606d3fa_0.so\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "P[16](lon=0.500000, lat=0.000000, depth=0.000000, p=10.000000, time=1.000000)\n",
       "P[17](lon=1.000000, lat=0.000000, depth=0.000000, p=10.500000, time=1.000000)\n",
       "P[18](lon=0.500000, lat=0.500000, depth=0.000000, p=10.000000, time=1.000000)\n",
       "P[19](lon=1.000000, lat=0.500000, depth=0.000000, p=10.500000, time=1.000000)"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pset = ParticleSet(fieldset, pclass=SampleParticle, lon=xv.flatten(), lat=yv.flatten())\n",
    "kernels = pset.Kernel(SamplePoff) + AdvectionRK4\n",
    "pfile = ParticleFile(\"tmp5.nc\", pset, outputdt=1)\n",
    "pset.execute(kernels, endtime=1, dt=1, output_file=pfile)\n",
    "pfile.close()\n",
    "pset"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "lon=0.5\n",
      "p=10.0\n",
      "time=1.0\n"
     ]
    }
   ],
   "source": [
    "print_nc(5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Information stamped at $t$ is:\n",
    "\n",
    "- Position: $x(t)$\n",
    "- Variable $p$ evaluated at $x(t-dt)$ and $t$\n",
    "\n",
    "**!Danger!: the interpolation does not use consistent position and time**"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}