test_fd_derivation.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "from pystencils.session import *\n",
    "from pystencils.fd.derivation import *"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 2D standard stencils\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "stencil = [(-1, 0), (1, 0), (0, -1), (0, 1), (0, 0)]\n",
    "standard_2d_00 = FiniteDifferenceStencilDerivation((0,0), stencil)\n",
    "f = ps.fields(\"f: [2D]\")\n",
    "standard_2d_00_res = standard_2d_00.get_stencil()\n",
    "res = standard_2d_00_res.apply(f.center)\n",
    "expected = f[-1, 0] - 2 * f[0, 0] + f[1, 0]\n",
    "assert res == expected"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Finite difference stencil of accuracy 2, isotropic error: False"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "assert standard_2d_00_res.accuracy == 2\n",
    "assert not standard_2d_00_res.is_isotropic\n",
    "standard_2d_00_res"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}0 & 0 & 0\\\\1 & -2 & 1\\\\0 & 0 & 0\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "[[0, 0, 0], [1, -2, 1], [0, 0, 0]]"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "standard_2d_00.get_stencil().as_array()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 2D isotropic stencils\n",
    "\n",
    "## second x-derivative"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Finite difference stencil of accuracy 2, isotropic error: True"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "stencil = [(i, j) for i in (-1, 0, 1) for j in (-1, 0, 1)]\n",
    "isotropic_2d_00 = FiniteDifferenceStencilDerivation((0,0), stencil)\n",
    "isotropic_2d_00_res = isotropic_2d_00.get_stencil(isotropic=True)\n",
    "assert isotropic_2d_00_res.is_isotropic\n",
    "assert isotropic_2d_00_res.accuracy == 2\n",
    "isotropic_2d_00_res"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}\\frac{1}{12} & - \\frac{1}{6} & \\frac{1}{12}\\\\\\frac{5}{6} & - \\frac{5}{3} & \\frac{5}{6}\\\\\\frac{1}{12} & - \\frac{1}{6} & \\frac{1}{12}\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "[[1/12, -1/6, 1/12], [5/6, -5/3, 5/6], [1/12, -1/6, 1/12]]"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "isotropic_2d_00_res.as_array()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAI8AAACYCAYAAADdsLqwAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAATu0lEQVR4nO2deXAcVX6Av9cjzaHDkqyWLMuysMFHwBsOh8O7wGKSAB4XXtaLWUPWWQjBfwAudiFsMBsvGw6zVMAmWyFkKY6FBGIoF4YE22OKdbKAKY4ymDXJrnFsgW1ZsuWRNNY1Ovvlj/HI0mg0M32Mprs1X5WrPJrp1+/N75vu169f/56QUpInjxGUXFcgj3MpGO8NEQpPB0oAoaM8CfTIoHrUbMWyjQiFa4Fi9LUvGRI4IYNqxHSlsogIhQuBOqBQ56YSOC6DaseYMhNPWyIUPh9YD8wxVk0AjgAPyqD6vokysoIIhS8GHgZmWVisBnwM/I0Mqq0WlmsJIhS+E7gFmGKwiCHgA+BeGVRPDpc7Uh4RCvuAXSZ2MpJ+YLGdvkwRCpcA7wNFWdrFLhlU/zpLZRtChMJLgSctKu4dGVTXxF8k9nm+iTXiAHiBxRaVZRWXkz1xAL4pQmGrvj+rWGJhWVeIUNgff5EoT5WFOwKotrg8s1jdvkQ8wNQs70MvVrbZC5THXyTKY/XVl9nOqNVMxNWl3a5grY6BMuY/GdHZrnDn4nqW18/lwO+8FlfK3ri17SbapU8ef7HGw68d5eKrO3Vt5wbc2nYT7dInT6EXpk4b0rsTV+DWtptol93Oz3kcRF6ePIbJy5PHMOPe2xqXtctncGifn6YGL9esinDtrWPuebgWt7bdYLv0y/PYG7a/6Zk13Np2g+3Kn7byGEb/kceN3LNkJvs+DSR9b+75UX75zpEJrpEjSC/P0qp5ukrcfmK/0crkjI070suR6ntwWpstiml6eV7bf4C1y+toavDx+FuHmHNeP3t3+Xn2gWoKCiUV0wa5/7lmCt0zYp+U7Sf2s3tngFc3VqJpgmW3tXPl9V25rpYhLIpp+j5PsuHrmlkDPL71CE++fYRp9QO892aJ6QbZnb6oYMvTU3l0SyMbQ0ccKw5YFtP08iQbvq6uG8JfFJtF5vFIhEv73ccP13OytRKAvbsCeP0a626o42crawk3eXJcO+NYFFNzUW88WMied4u5bFkX/b0+Wo/V0N/rM1WmXejv9THQ76f7ZDlSQluLh+OHvTyyuZElq07y4no111XUhZTQ01lM67EatKHxp2mMjGkajF9tdUUUnrijhp883YzXB9FuD9GuKfR2l1JQ2E9JRSs+n3MPSV2RcgCkFES7iykp05i/MIrXBxdd1cPmp+w26Ss5fb2F9PVU0BWZipSn4lF9POlnE2OaBmPBHRyA9bdO5wf3tjLrnAEA/EU9eAoGkVIw0O8j0lLD8w/+QghRYWgfuUTTBNHu2HRSKRW6IhUsuCRK40EvUoMvP/NRUz+Q41qmRQixlI/fvo7OdhVN8yCBotKTKMrYh/WSxTQNmR15EoevFQ8c3Otn04ZKNm2oJHhzhKtu6qSkvJWO1mqkFEipUODtBGz/JY+hp6OUkV9vf28RU6cpLFrSxd3BmSgC7nnqWM7qlzlHEcSOnhCbU1ha3gZkHtMUZCZPsuHrpTePvf9RPKWDjtbT85YHeiuklM67KvEUDuL19tLfFxs49Pl7EIrGijURVqyJ5LZyuvgaTTvdsS/0RSnwDgKZxzQF1vZJFEVSXNZGobePyppGhob8Qojk51c7EyjuoXpmbODQF+imqq4Rj0fLca10IYQoAyIAqLWHUTyDTJlq6WNQifKYP8WUq61Mqz9EoKSHy6/7BVBtI4Em4hSa89P0KHEuvXYT/qJeamc34C+KWlB8f/w/ifL83oLCT3PupduARdhHIGvbN5YI0JTlfaRklDhQREHhFxYWfxwYPnqNkkcG1X3Enqi0gr3AR1LKj7GPQJ8Dn2Sx/BdkUM3ZPOdEcaSUUeBVwKp5R8/JoDp8KZHsWXUfsBK4DGNPj3YDHwKvyqA63FkWQlwCfAS0SCmnGSjXEkQoHOB0+1IPwb/7xkpKy5tZeOV7KT4lif0iQzKovm1ZRXUyjjix90Lh2cCNwDnoT3SgAceAbTKo7hy1z4nMz2MXgTJFCCGBHVLKYK7rkopU4mSTCR0BttkpzBXkShzIwUzCvEDWkUtxIEfTUPMCmSfX4kAO5zDnBTKOHcSBHE+AzwukH7uIAzZ4eiIvUObYSRywgTyQFygT7CYO2EQeyAuUCjuKAzaSB/ICJcOu4oDN5IG8QCOxszhgQ3kgLxDYXxywqTwwuQVygjhgY3lgcgrkFHHA5vLA5BLISeKAA+SBySGQ08QBh8gD7hbIieKAg+QBdwrkVHHAYfKAuwRysjjgQHnAHQI5XRxwqDzgbIHcIA44WB5wpkBuEQccLg84SyA3iQMukAecIZDbxAGXyAP2FsiN4oCL5AF7CuRWccBl8oC9BHKzOOBCecAeArldHHCpPJBbgSaDOJAm0YEIhc8BbgG+QWxZ5PEYAL4E/k0G1U+trKBZEpMriFC4HFgNfItYlozx08rueffPKJ7SyrwLPk+xCw1oAULAv7O0ago2FkeEwiuAZcB0xj94SGLZTj4AnpVBNZK0rPHkEaHwWcBmoFhH3fqBH8qgukfHNllnWCDF08LWYx8CZ2e0YeOBebG0cjMyW1KoK/Iy35/70KlXdhTnduDHOjfbByyXQXVMWr1Up63voU8ciB2dVurcJusMn8LmnldNU8O1WdmJNqTQ/PVDeArAnuIIYJWBTf8IuDDZG6nkmWtgRwD6VlSZIKSUH7Pqvr9H0zw0NZxlaeHakELTV3MIlMAT2+bZTZxTVABGs9YnjWkqeYxmh7fvGl5/8qeHqJpx2FKB4uIA1J75f8xfaNfls83EJWk2MX0FHm0o4O5rzqDurFhGzHUvNTlurXFfoJeqGYc5cbSepoazqD3zYEbbJWt7uSpHiZMss7qdMRlP/TaefWGUBzflNOOnaYwKNLLtiUccp4kTx0Q89Y/z7N8T4MdXz+SZv1ORjsprPZq4QHpOYfG2/+qnVRxtcL44YCqe+uSpqh3i+d0NPLnjCCfDHv77dWcv0qZHoHjbN2w7StuxCj77rfPFMRlPffJ4/ZKikthCXpde20XDF85fWytTgbx+iT8gaD40hwsWQ8uRNkeLA6bjqU+ersjpz3/xYYDaM3OeKt8SMhGoo61guI/z9R/amXFmf9LPOQmT8dTXYf7stwE2bVDxBTSq6wa47cGwru3tTKpOtDak8P5/nMlbz4OvKMq0Og+rH7Iqq3ruMBlPffJ8+7vdfPu73Sk/091RJETVE0CdlPJGXeXnmkSBILZWVdNXczjvcgje7Lg+jhCiFniF2Qt28E87FZQRq/dkEs8UWDOgF1u/soSuk1P5+g9zgG8Dzhr/iTNSIIgt1AZO7hyfAVxCe8vFNH8dwFfURWl5O75Ar9mCrZEn3FxLX0+spx673FOAwlPp9+3Dn98I3//R2L8/cQfsH+de7lnnwn3PQFNDZrdrHrjxoDh2yEQls4SUsX+93aX0dpdSN2e/2SL1ybO0avz7Vq/8b4SezjIURQPiq8t9y3jVssDi65ei1v5wzN8fe3P0aykVWpvrhl8LRWP1JeNfXGw5dHj4/wsX38P2l+w0iLoQ+AcggBCgKEMUl7UDqeO5/URaufTJs/3EfnbvDPDqxko0TbDstvZRi9OXqSdQPPF1KUullB/qKj/LiFB4PpD6cD1y5Big0NfLQJ+fDduH+Je1A3gKJIpHsvbZZqpmjD01r3lij9z2YoPVdTeKEKIVGMTn30bl9PPwBaKIU1OY0sUzDfrk6YsKtjw9lUe3NCZdOllRJNNnHZVSfkeIDFZ1txuJtxyaGuaiKEOxy/ihev72Gaibc5CtL0xh26/LuGVdW45rnBYp5X4hRAUv7lFJXEstXTzToC/Ae3cF8Po11t1Qx89W1hJu8oz3USkddu8i1b0qX6CXafWHgdg4ULRL4YyzHTPOM24sdMQzGfrkaWvxcPywl0c2N7Jk1UleXG90foi9yOQmpy/QSyTczCN/5WH7S9XMX2j6aiXnmIynPnlKyjTmL4zi9cFFV/XQeCDVvGZnoOfu+DcWdfLkjsNctxr+9dHZE1XFrGEynvrkWXBJlMaDXqQGX37mo6be2bcn9IjT3xvrZfoCvVTWtuD1C8tnJE40JuOpr8NcXqWxaEkXdwdnogi456ljura3E3rn4+zb7eOFh6tQFCj0ady18SiaNkPXfCC7YTKe+gcJV6yJsGJNRPd2dsLIRK5zL+vlH98+MupvfVH9E8rshol4Ou9y2ixWzgA0MqHMRaSSx+iltn3vabW3+CyfOjpWILu238zQSdI2pZLH6BB7s8HtsooQoozHb/8VYP1NzrhAA/0evjfrE8vKtZYI0GNw26QupJJnK7HHTvXynwa2ySrDz47//hOonf1JVu6O+wK9FPo209ud8+wcyZBBdZDYI9F6iRB77HgM48ojg+onwP1ApreIm4D1Mqju0Fu7bDIq6UB/bxGK5xZgN+YO44lEgbeonX0rNknvMg4PAa8Tew49HRqx7+kWGVSTPsSYMtHB8IdC4VLSJDqQQdV2M+tSZasQobAfKCJVooOV88IUenfy8v+keoRaAzpkUB3uFyQmVzDeguwgQmEPUMrp2Q+JSKBHBtWUo+gZyeNErEhzcmo+0g4pZdDAtrYWyApcealuh/w4dkgwlW1cJ48dxInjdoFcJY+dxInjZoFcI48dxYnjVoFcIY+dxYnjRoEcL48TxInjNoEcLY+TxInjJoEcK48TxYnjFoEcKY+TxYnjBoEcJ48bxInjdIEcJY+bxInjZIEcI48bxYnjVIEcIY+bxYnjRIFsL89kECeO0wSytTyTSZw4ThLItvJMRnHiOEUgW8ozmcWJ4wSBbCdPXpzT2F0gW8mTF2csdhbINvLkxRkfuwpkC3ny4qTHjgLlXJ68OJljN4FyKk9eHP3YSaCcyZMXxzh2ESgn8uTFMY8dBJpwefLiWEeuBZpQefLiWE8uBUqaVk6Ewn8MXE7sYfjxEwGMRRLLAfMh8KkMqsMPwttJHBEKX0BsaYMSUrXvhrugrPJMEQqvTVGcBhwHfiOD6lFLK5ohUsqPhRCLgI+EEMcTn40XoXAAuAJYwDgrFadAI5Zz6TcyqI7KvTQm0YEIhe8FVuvcQTJel0H1p2A7cdYBf5nRhxsPzMMX6KZqRiZSDAJ3yaC600z9zJAsuYIIhcuBl4HMFl4Zn37gdhlUd8X/MOq0JULh6cBtJncS53oRCi+wmTizyVQc/RQQy2eUM8Y5hf0F5sWBWIqd+0b+IbHPcyH6TlOpiZy4ApuIc4qLs1z+zFM/wJyRRKCLLCx+ngiFy+IvEuUpMl18d0cp7S1VDA0pvPL4P8fLtYE4AIEJ2Id/AvaRkpEC8cHWmxga9BBums5An97+TjKGHbFmsbY4UsLJcDVSKnR3VMR3ZhNxMmegv5CezlIA+qLFdLRVUFrRPrzUUGqsO3KbYLgTPTjwEc1fx9L8CkVSWWM28fpw+/Rdqne2K9y5uJ7l9XM58LuxaeaiXcVIKZAytoPKmk+BPnN1zQHRrhI6208v4tHRVkV7iy9l2+1JAx7P6SUBertLGRoaHfN0MU2BPnn8xRoPv3aUi6/uTPp+Z3slUp4us/XYQiCnfQBDFE3pYOQSlwWF/ZRW9KVsuz1ZhqaNPrt0nywb9TpdTFOgT55CL0ydljxJdX+fl4F+P0JIhJAUT4nwndU/l1LmZOzDFAUFQ3j9sZzFQkhKytpStt2+/JqzL/ovfIFuhJBIKeiKTGXk8IyJdlnX5xkaKMBTMEBJeRvFUzpQFElFdatl5U80JeXtwysbF01x0tFmGCmlFKFwGDjK4EABnZEKertL0IY8eApM/xCskydQ0kOg5CvLyss1/qIehKLhC3Q7dEns0RQUDlJRdQKqTlhWpFUFuQ4hoHa2M1eymSD0y7N2+QwO7fPT1ODlmlURrr3Vdsm7s4Zb226wXfrleewN53WArcKtbTfYrpzPYc7jXPJ9HoB7lsxk36fJb13MPT/KL985kvS9SU56eZZWzdNV4vYT+41WJmds3JFejlTfg9PabFFM08vz2v4DrF1eR1ODj8ffOsSc8/rZu8vPsw9UU1AoqZg2yP3PNVPolBF7g2w/sZ/dOwO8urESTRMsu62dK6/vynW1DGFRTNP3eZINX9fMGuDxrUd48u0jTKsf4L03S0w3yO70RQVbnp7Ko1sa2Rg64lhxwLKYpj/yJBu+rq47/drjkYhJ0O/euyuA16+x7oY6fAGNHz15HLXWabcrYlgUU3NRbzxYyJ53i7lsmXN/hZnS1uLh+GEvj2xuZMmqk7y4Xk2/kQPREVPj8nRFFJ64o4afPN2M12e4GMdQUqYxf2EUrw8uuqqHxgPu6+TpjKkxeQYHYP2t0/nBva3MOmcg/QYuYMElURoPepEafPmZj5p6d7XbQEwzG+dJHL5WPHBwr59NGyrZtKGS4M0RrrrJkXeeM6a8SmPRki7uDs5EEXDPU2Zn5OUWC2KamTzJhq+X3uyO+zp6WLEmwoo1kVxXwxIsiGniacv5Uw9SMxHtmzTfYaI8Vp967HZ0moirQrudvq2sjxxZXqI8HwJWjl28b2FZVrCL7B4Z9smgGs5i+UawMgafy6A6/AMcJY8MqhHg55gXSAM2yqB6yGQ5liKD6nHgUaz9gcSJAOuyUK5ZXsMagVqJuTHMmGfVYfj55kXAFIwlOvjIhr/AYUQoPJVY+1InOsiMeKKDj2VQte1jRiIUno+xRAdDwDFi7Rt1CZ9Unjx5MmES3JTKky3y8uQxzP8Dw8XXERVoKYcAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 144x144 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(figsize=(2,2))\n",
    "isotropic_2d_00_res.visualize()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "expected_result = sp.Array([[1, -2, 1], [10, -20, 10], [1, -2, 1]]) / 12\n",
    "assert expected_result == isotropic_2d_00_res.as_array()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Isotropic laplacian"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}\\frac{1}{6} & \\frac{2}{3} & \\frac{1}{6}\\\\\\frac{2}{3} & - \\frac{10}{3} & \\frac{2}{3}\\\\\\frac{1}{6} & \\frac{2}{3} & \\frac{1}{6}\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "[[1/6, 2/3, 1/6], [2/3, -10/3, 2/3], [1/6, 2/3, 1/6]]"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "isotropic_2d_11 = FiniteDifferenceStencilDerivation((1,1), stencil)\n",
    "isotropic_2d_11_res = isotropic_2d_11.get_stencil(isotropic=True)\n",
    "iso_laplacian = isotropic_2d_00_res.as_array() + isotropic_2d_11_res.as_array()\n",
    "iso_laplacian"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "expected_result = sp.Array([[1, 4, 1], [4, -20, 4], [1, 4, 1]]) / 6\n",
    "assert iso_laplacian == expected_result"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# stencils for staggered fields\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "half = sp.Rational(1, 2)\n",
    "\n",
    "fd_points_ex = (\n",
    "    (half, 0),\n",
    "    (-half, 0),\n",
    "    (half, 1),\n",
    "    (half, -1),\n",
    "    (-half, 1),\n",
    "    (-half, -1)\n",
    ")\n",
    "assert set(fd_points_ex) == set(FiniteDifferenceStaggeredStencilDerivation(\"E\", 2).stencil)\n",
    "\n",
    "fd_points_ey = (\n",
    "    (0, half),\n",
    "    (0, -half),\n",
    "    (-1,-half),\n",
    "    (-1, half),\n",
    "    (1, -half),\n",
    "    (1, half)\n",
    ")\n",
    "assert set(fd_points_ey) == set(FiniteDifferenceStaggeredStencilDerivation(\"N\",2).stencil)\n",
    "\n",
    "fd_points_c = (\n",
    "    (half, half),\n",
    "    (-half, -half),\n",
    "    (half, -half),\n",
    "    (-half, half)\n",
    ")\n",
    "assert set(fd_points_c) ==  set(FiniteDifferenceStaggeredStencilDerivation(\"NE\",2).stencil)\n",
    "\n",
    "assert len(FiniteDifferenceStaggeredStencilDerivation(\"E\",3).points) == 10\n",
    "assert len(FiniteDifferenceStaggeredStencilDerivation(\"NE\",3).points) == 12\n",
    "assert len(FiniteDifferenceStaggeredStencilDerivation(\"TNE\",3).points) == 8"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "c = ps.fields(\"c: [2D]\")\n",
    "c3 = ps.fields(\"c3: [3D]\")\n",
    "\n",
    "assert FiniteDifferenceStaggeredStencilDerivation(\"E\", 2, (0,)).apply(c.center) == c[1, 0] - c[0, 0]\n",
    "assert FiniteDifferenceStaggeredStencilDerivation(\"W\", 2, (0,)).apply(c.center) == c[0, 0] - c[-1, 0]\n",
    "assert FiniteDifferenceStaggeredStencilDerivation(\"N\", 2, (1,)).apply(c.center) == c[0, 1] - c[0, 0]\n",
    "assert FiniteDifferenceStaggeredStencilDerivation(\"S\", 2, (1,)).apply(c.center) == c[0, 0] - c[0, -1]\n",
    "\n",
    "assert FiniteDifferenceStaggeredStencilDerivation(\"E\", 3, (0,)).apply(c3.center) == c3[1, 0, 0] - c3[0, 0, 0]\n",
    "assert FiniteDifferenceStaggeredStencilDerivation(\"W\", 3, (0,)).apply(c3.center) == c3[0, 0, 0] - c3[-1, 0, 0]\n",
    "assert FiniteDifferenceStaggeredStencilDerivation(\"N\", 3, (1,)).apply(c3.center) == c3[0, 1, 0] - c3[0, 0, 0]\n",
    "assert FiniteDifferenceStaggeredStencilDerivation(\"S\", 3, (1,)).apply(c3.center) == c3[0, 0, 0] - c3[0, -1, 0]\n",
    "assert FiniteDifferenceStaggeredStencilDerivation(\"T\", 3, (2,)).apply(c3.center) == c3[0, 0, 1] - c3[0, 0, 0]\n",
    "assert FiniteDifferenceStaggeredStencilDerivation(\"B\", 3, (2,)).apply(c3.center) == c3[0, 0, 0] - c3[0, 0, -1]\n",
    "\n",
    "assert FiniteDifferenceStaggeredStencilDerivation(\"S\", 2, (0,)).apply(c.center) == \\\n",
    "       (c[1, 0] + c[1, -1] - c[-1, 0] - c[-1, -1])/4\n",
    "\n",
    "assert FiniteDifferenceStaggeredStencilDerivation(\"NE\", 2, (0,)).apply(c.center) + \\\n",
    "       FiniteDifferenceStaggeredStencilDerivation(\"NE\", 2, (1,)).apply(c.center) == c[1, 1] - c[0, 0]\n",
    "assert FiniteDifferenceStaggeredStencilDerivation(\"NE\", 3, (0,)).apply(c3.center) + \\\n",
    "       FiniteDifferenceStaggeredStencilDerivation(\"NE\", 3, (1,)).apply(c3.center) == c3[1, 1, 0] - c3[0, 0, 0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAPsAAAD2CAYAAAAQ7PF4AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAAO2ElEQVR4nO3de2xU14HH8e8Zv8Bg/Bq6NsZgAqYFqoglEdk2AUpK0wyhiJdIuquqqA27TdCKJis1StA26i7ZPNRtts32sdImbMWqG8EGVSnK9JEmStRo1bTk0WxLaMxrndAAY+MBY3tsz5z9Y8AYPJ4ZwzzunfP7SBGemXt9jnT99T2+84ix1iIipS9Q7AmISGEodhFHlBd7AjKWCUcMsAJYBtQCJk9DWeA08HMbCr6ZpzFkFBOOVAG3ATcBk3P4rWPAG0DYhoLnU46tv9m9x4QjDwN/WeBhd9pQcHeBx3SKCUfKgV3A0jwOcxD4q1TBaxnvMSYcaQY+X4Sh/9aEI2VFGNclt5Df0AEWALenekCxe89i8rdsT6cWmFuEcV3y5wUa54ZUdyp276lydGwXTCrQOCmPo2L3o463K/na2plsnjeX1dPnZ9z+0IEq7lk2m3Uz27ln2WwOHVDUXrH3qTruXT6btTPaefTuprTbPvuteu766Fw2ts3jsa1NDA5MaAWo2P2ovMJyy+fOse2Jkxm3HYzBzi0trFh/lj0dHazcFGXnlhYGYwWYqGTU2DzM5u1drNhwNu12//NCNT/+twZ27u3kmQNHONVZwdPfaJzIUIrdj9oWDrF2a5Q5izIX+8ZL1cTjcOdXz1A5ybJ5ew/Wwm9frC7ATCWTWzf18qmNvdTUx9Nu9+KztazcFGXe9YPUNia46/4uXtlXO5GhFHupO3awitb5McyoQ906P8axg1rK+8n7HZVc9/FLv9zbF8c4211Gz+msG1bspa7/fIDqmsRl91XXJOjv1bH3k4G+AFOmXTr7T61Lfn3+nGIvKT/dXcP6We2sn9XOA+taJrTv5Cljw+7vDTB5amKcPcSLJlUnOH/u0usgzkeTx3RKTdbHUS+X9YPbv3CO279w7qr2bVsQ4ydP12MTjCzlO/9YxZovncnhDCXfZs4b5Mj/VgHJn4P33prEtIY4ddOzjl1ndj+yCYj1G4YGk0+9xPrNuE/DLLm1j0AA9nynjsEBw96n6gC4cVVfoaYraQwPJY9fIg6JePLr4aGx2336zigv763l8DuVRLsCPPtkAys2RCcylM7sfnTiWDlbb7pu5Pb6We00Ng+z+3dHAHhgXQsLl/bzxYe6qayCHbs+4Nv3NfGjb05nxpxBduz6gEpdn/OEHz7SyHPfvfQU2mv7p7FxWxert0S5d/kcvvfqUZrbhvnkHX0cf7ebHZtaGRwwLL2tly8/3DWRofRGGI8x4cg64PEiDb/JhoLvFGnskmfCkQeBLQUYar8NBf/uyju1jBdxhGKX0bTMKw0pj6Ni957+Io49UMSxXVCoi6Ipf4YUu/ccANK/dHI8XX9qItrVcJXjRoAjV7mvZOf1Yo6j2D3GhoIR4PsT3nF4qIz+89PojTYw8YuuCeAxGwrqhTb59Trw82KNoavxHmXCkUXAcrL9DLqXn1vJ/x1aiTFxrr9lL4uX/SGLYRIkP4PulzYUPH5NE5asmHAkwHtvraJ98QKSn0GXqw8qGSD5GXSv2VBwOOXYit3/jDEB4ENg+oW7fmWtXVbEKckVjDHlwDpgB7AIaLLWdhdyDnpRTWlYxeWfgnKjMabNWnusSPORC4wxQeArwH1ABVADvFno0EF/s5eK+0n+EF0UAP6mSHMRwBiz2BjzX0An8BDQQPIYnQMeKcqctIz3N2NMM3CUsZ871gNMt9am/PtNcu+Kpfp8ksfkyk/s7SK5hC/4cdGZ3f/uJvWLKMqANQWei+vWAHtIfkJwNWND7weeLNYvYJ3ZfSzFhbkr6UJdARljKoDnST6LkupjvwaAVmttpKATu0Bndn9bBUxL8/gnjDFtBZqL86y1Q8Ba4FXGvhoxDuwrVuigq/F+1wP8ftTtJRf+fePCv8NAijdHS75Ya4eMMW+R/L+yDHDpWZIY8ESx5gVaxpcUY4wF3rXWLij2XFxljPkH4O+Bx0j+7X5xSf+mtXZJml3zTst4kRwZFfo/Wmsf5NKS3lKkp9tG0zJeJAeuCP3rMLKkXwvcAewv5vxAsYtcs1ShX3Thot2PizGvK2kZL3IN0oXuNYpd5Cr5KXRQ7CJXxW+hg2IXmTA/hg6KXWRC/Bo6KHaRrPk5dFDsIlnxe+ig2EUyKoXQQbGLpFUqoYNiFxlXKYUOil0kpVILHRS7yBilGDoodpHLlGrooNhFRpRy6KDYRYDSDx0Uu4gToYNiF8e5EjoodnGYS6GDYhdHuRY6KHZxkIuhg2IXx7gaOih2cYjLoYNiF0e4HjoodnGAQk9S7FLSFPolil1KlkK/nGKXkqTQx1LsUnIUemqKXUqKQh+fYpeSodDTU+xSEhR6ZopdfE+hZ0exi68p9OwpdvEthT4xil18SaFPnGIX31HoV0exi68o9Kun2MU3FPq1UeziCwr92il28TyFnhuKXTxNoeeOYhfPUui5pdjFkxR67il28RyFnh+KXTxFoeePYhfPUOj5pdjFExR6/il2KTqFXhiKXYpKoReOYpeiUeiFpdilKBR64Sl2KTiFXhyKXQpKoRePYpeCUejFpdilIBR68Sl2yTuF7g2KXfJKoXuHYpe8UejeotglLxS69yh2yTmF7k2KXXJKoXuXYpecUejeptglJxS69yl2uWYK3R8Uu1wThe4fil2umkL3F8UuV0Wh+49ilwlT6P6k2GVCFLp/KXbJmkL3N8UuWVHo/qfYJSOFXhoUu6Sl0EuHYpdxKfTSotglJYVeehS7jKHQS5Nil8so9NKl2GWEQi9til0Ahe4CxS4K3RGK3XEK3R2K3WEK3S2K3VEK3T3lxZ6AZMeEI2WASbtRoAyMwYQj6Y/rmqavo9CLwoQjBijL4xBxGwralGNbm/J+8QATjkwGHgA+CzRk3OH9jvmUVwzSNPvYuNv0RBrp7WnEmOO0zP0W8K82FEzkaMoyDhOOfBb4a2Ah+V1RW6AD2GVDwedGP6BlvLc9CnyebELPxsXQp9Z10TI3Bmy78J/kkQlHbgD+Bfg4+W/OAO3AP5lw5DOjH1DsHmXCkRrgMxk3zNbo0OuCXaMe2ZCzMWQ8aylOa+tH31Ds3tVErq6pjB86wAwTjujnIL9mFWnc1tE3dJC9KzfHJn3ouR1LxpPPC3JZj6uD7Dd7n6rj3uWzWTujnUfvbkq77e7HWvjKzY1sX2X5wYMVDA6kv5ovOWOMaTHGtGW1ccfblXxt7Uw2z5vL6unzM25/6EAV9yybzbqZ7dyzbDaHDlRlM4xi95vG5mE2b+9ixYazabd7+b+b2P/MFB56podnDhzmVGcFT3+jsUCzFNgHdBhjfkXne62ke9arvMJyy+fOse2Jkxm/62AMdm5pYcX6s+zp6GDlpig7t7QwGMu4q2L3m1s39fKpjb3U1MfH3aYn0sgr+6axfH0/1998itrGBHfd38Ur+2oLOFPXlZNcRt9Mx9s3c+LoPHpOBxkeHHsdpm3hEGu3RpmzKHOxb7xUTTwOd371DJWTLJu392At/PbF6ky7KvZSk4iX0dvTyImjceYv7hm5v31xjLPdZfSc1jEvtES8ApsI0Hu2npOdczjZ2Urfualpz/bjOXawitb5Mcyow9g6P8axgxmX8noFXalJJMoIBOLE+suwNsiHx5NL9+Gh5OOd781ioO/yfbb+xTsmPlzgiZa8j429yxosMBSbzJlTVZw5DVNru6lt7M76u/afD1Bdc/mLoKprEvT3ZvwlrthLyZRpZzh/tp5EooyqydB3tmIk8t5o8t+KqsqR+y6yNsUPpuTNr38G//l4Ms72xXU88Xz2sU+eMjbs/t4Ak6dmfBWkYi8l9R85Tf1HTgPQ2t7MyfeHmDkvAsBvflHNtIZmPnbD4TH77f9wkQ0FdWrPIWPMAWBJigcSfCJkWXXXGabURikvH//aSyptC2L85Ol6bIKRpXznH6tY86UzmXbV329+MzwEsX5DIg6JePLrK8/UAJ++M8rLe2s5/E4l0a4Azz7ZwIoN0cJPWDCBOMZYqib30dB0guY5h6lt7B4J3SaSx3FoMPnUaKzfjPs06ZJb+wgEYM936hgcMOx9qg6AG1f1pdx+FJ3Z/eaHjzTy3HcvPYX22v5pbNzWxeotUe5dPofvvXqU5rZhPnlHH8ff7WbHplYGBwxLb+vlyw+P96IayY8YEGXGnBP82eyacc/iJ46Vs/Wm60Zur5/VTmPzMLt/dwSAB9a1sHBpP198qJvKKtix6wO+fV8TP/rmdGbMGWTHrg+ozPxUu9715lEmHPko8HyBhtMyPseMMUuBWuCXvHD6P4CbijCNwzYUXH3xhpbxInlgrX3dWvsLa61n3j6s2L2rkD8kWt7lV7GCv2xcxe5dEQoTYbcNBSd2RVgm6pQXxlXsHmVDwTPAbwow1M8KMIbrflqkcS87tord2+4Hfk1+zvDDwAvA43n43jKKDQVfAv4ZSP/mpdw5D/w7sGf0nboa7wMmHKkHguTul3MC+JMNBXtz9P0kCxc+CLQVqMzjMMPA+zYUHPOmGsUu4ggt40UcodhFHKHYRRyh2EUcodhFHKHYRRyh2EUcodhFHKHYRRyh2EUcodhFHKHYRRyh2EUcodhFHKHYRRyh2EUcodhFHKHYRRyh2EUcodhFHKHYRRyh2EUcodhFHKHYRRyh2EUcodhFHKHYRRyh2EUcodhFHKHYRRyh2EUcodhFHKHYRRyh2EUcodhFHKHYRRyh2EUcodhFHKHYRRyh2EUcodhFHKHYRRyh2EUcodhFHKHYRRyh2EUcodhFHKHYRRyh2EUcodhFHKHYRRyh2EUcodhFHKHYRRyh2EUcodhFHKHYRRyh2EUcodhFHKHYRRyh2EUcodhFHKHYRRyh2EUcodhFHKHYRRyh2EUcodhFHKHYRRyh2EUcodhFHKHYRRyh2EUcodhFHKHYRRyh2EUcodhFHKHYRRyh2EUcodhFHKHYRRyh2EUcodhFHKHYRRyh2EUcodhFHKHYRRyh2EUcodhFHKHYRRyh2EUcodhFHKHYRRyh2EUcodhFHKHYRRyh2EUcodhFHKHYRRyh2EUcodhFHKHYRRyh2EUcodhFHPH/Ibksjy4emPAAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "d = FiniteDifferenceStaggeredStencilDerivation(\"NE\", 2, (0, 1))\n",
    "assert d.apply(c.center) == c[0,0] + c[1,1] - c[1,0] - c[0,1]\n",
    "d.visualize()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "v3 = ps.fields(\"v(3): [3D]\")\n",
    "for i in range(*v3.index_shape):\n",
    "    assert FiniteDifferenceStaggeredStencilDerivation(\"E\", 3, (0,)).apply(v3.center_vector[i]) == \\\n",
    "           v3[1,0,0](i) - v3[0,0,0](i)"
   ]
  }
 ],
 "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.9.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}