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": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAGkAAABLCAYAAABz5qkHAAAACXBIWXMAAA7EAAAOxAGVKw4bAAAFRElEQVR4Ae1dUW7UMBDdRXyj0kqI7+0NlvYElBuUI5TegIoj0BtQToDaI7QngO4NyjdCoq24wPJe1rOysk5iB3vtJmPJtWM74/F7njixve50MpnM4W/hXe5quVy+d2VoWjwEptPpHaTNXBKB//S5lXGBOAvb7qd9ofFkCHyG5J2a9ENcHzNtCi+WtA/WlBSiUoCDdZGgy7oldaqGG8n4H/hH+H34bxCyQJjdDV03WtISfgbAJ00e+Ry3jiUfcZpn9SyVtFzhEHUj1vCAdDl5hkinQy/9gEIk8UoKI05r4vWlpOUIx6CbF0kAn2941w4SviNtDqDqg56jaLKkwevmS9JRA8S0Jrqm/FVu2r9NdQ9Gt06SLCu5b8F61pKXLGssunWSBIR3DcrSM23Q24izy6WKj0I3H5IE4LZxZ08KZQoHrZsPSW3WIj2Z30453Ch06yTJvGqTACHEJkN6cJaZirHo1kmSYYSv30KITZIQ53o9t8uljA9eN1+S+MF64ECaU0MLq0c7iiRPGoVuc8DoMy3EKaD6tNAD0uYgqXE6aRt50GFwuhFreMC3nNhLFUhrdW+Q+wnfJjNTilPpbyFk0XrXdjIHrZsuVWynEwXXAmNYL1X4jknBlegN8RBQkuJhmUySkpQM2niClaR4WCaTpCQlgzaeYCUpHpbJJClJyaCNJ1hJiodlMklKUjJo4wlWkuJhmUxSyNxdMiVSCzZ7Ibixk262CiZnhcw7GnWag2BL4pwS/INpeLPkQnKMnl9ByKnx76AaJ4Vvkde00yi69qirN25eliQNheZcgWVPdC0ARm9YJIG0oBNbFsg6Q5u44ZNrUS/tvJjxmLh5rSehYdWaERrBxnH9aUfSSg6hJ9e87uo6Io0Eda6j1e/re426gnBD+fV6UvDjDjc/NVdZv+nVLt2Lfyp4Pe5cLXsqaej5XBB0OT5B6LJsollV7fd3DJa0gQSsii8MHFvPQeLjRoHCEkZJEjj4As+fmp4VxodTneIfd2YsuXFq35x4AgKcey8gjwRdI/+0+fayckjSC/jfJixLO2hjHkdN40qQviDoI294IgS9hqrkpfoR2V+Er+AZDtaBIL7S8nfBawtC2oy+0Eb/gl7kxe+XfoU2wlstEME3uUObIHMziSve9RmTiv+usFE3lsIP1yvEZf6ObdiF58bOc7t8wnhv3LxJQgM54LJhR6YhN0j7gfgtGsozIEp1JIiPtGo8qimZ/BspBm7eJDkeFbX2lnkJvaO8dPRtXQzcxvqd1BfzLPcpSVlgD6tUSQrDK0tpJSkL7GGVKklheGUprSRlgT2sUiUpDK8spZWkLLCHVaokheGVpbSSlAX2sEq9p4Uo1kxQ6smRYRhHwY3T+J1bm1BGT47scRRCX9xwX9iWLliQnhwJ1EJdLNx8x6TBn84YSoBn+Si4+ZIka0h13WQ7VFN+vXyK66a6B6NbJ0kwWVlRvG9BeNaSlyxrLLp1kgSEuRpLJz1zdbX620acXS5VfBS6+ZAkAItFybUd7tkXGeKD1s2HpDZrkZ6sJ0du9sxouHWShDV6ecwJIbY60oOTb+iwK5X4WHTrJMkAMvjTGYX4yGEU3HxJ4raoA0cD9ORIByhWUjTc5hDqMy00uNMZ8bisfr2YMgS2vXDDfetpoZAJVu5f05MjAUKg+2/c9OTIQMS3VRwf6rSk6p9c+Y5J29JN63EgoCQ5QCktSUkqjRGHPkqSA5TSkpSk0hhx6GOTdIc3imXN82NMXWIEgPkG9qhyjT2/kzjv1vRT+UVi/VT8CgH+AlHmQTcw+QfKv9HtYFDs7AAAAABJRU5ErkJggg==\n",
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}0 & 0 & 0\\\\1 & -2 & 1\\\\0 & 0 & 0\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "⎡0  0   0⎤\n",
       "⎢        ⎥\n",
       "⎢1  -2  1⎥\n",
       "⎢        ⎥\n",
       "⎣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": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAIIAAABNCAYAAABqvbycAAAACXBIWXMAAA7EAAAOxAGVKw4bAAAH5ElEQVR4Ae2dT47kNBTGq0esUYsFQiyrN6x7NCeguUHPcAJmbtAtViwQGvXcoAchViygkThA96xYDswNZtiwQEjQrVkhsSh+X8oOSSpV5YrtxEnZ0qvYjmO/973n/3FltlgsZr40m82OoWvffFJ8PlXZQvP1Dhl2dgcHBzKAJyaDeeeMEnwwVdli8iVlLtbQlUst5dkT6LVL2rGlSVW2XfmSfqBWPUsn1RbhuUnMpXRvSl/2jB2BCwQ4bAjxgPCp4qqGcIFlZMU3kJpKEN2qotcc3YyMoDCEe7U7ObC3CGRD2FvV1wXPhlDHY29D2RD2VvV1wUMZQnM0Wi9l3KFUZQvKl5chMOqcQ5foWVMT+a+gs3Hrfcl9qrLF4usAsbWg9Ct0lKePSyPYl1+MSlNHLRoeeLUI+wLYPsiZDWEftOwgY3VlcWNymhGtU3s5NUFeGWx5GB41gHqxJVnz9mfw9aoa6SJrbFnEjwsfVb7b/K58OhuCa4ZtzPQVB493lHXft7xUZO2TD2dDqIKLpWqAqb2JT6rx8nPPbm7MCWrv4twoSLdH54w84vsv6A5ZVtbs+xAqNuY7GYJhZu37Bwa0S8AqNq8IXwGSmmrvWtoH2M0y4F+zKRnyDX4ZtsK9GgLlbnzng/uqeEEwt+8jzBHY6W0lCj6BVt4/UBxU5oN/57xdeYidDt61HlK+j4Ff44+z2OWuy5+yg2NOnpo+UuRiFmPWoJpjnfpsuWrcMib9389h8XvLJmCpW3hmw4ldq/h2wnynrmGb8AB11EijFkHul+VlVL9qATTmsSulki258U4ozIMaQoua1X8lB14Ln7UoMx5Q3AOAPpeHODXNxQqswgm7TpjH6BoKjMwg5ibh5nSTLv82N69tIuS4wa/9FBlEks4H8yiGAEOPQeoQ8OwMI0ng1jEF37afbUtS7Y/b7g8S54t5cEMwNUYbWIURKAzZscIgIHUsVKuNbUpP7r3OEJh3NQQ7kDqGCW0966pm8yvA+xZ6if9UhF99bHLgwdM29w0JvkQGK9vXhH83XcS2Z2Pcj46581wf6VRD9P6B1gu09/Czudr35f9thIt4wHNan0gpHXJoLGDl0vU3SN1dr7JQZjTMybtcR9hp1gAIqtllv09tEViyVMW/4f6mvpUko3MP4XhQ2frCfCdDWKNGGUBt925NujFGpypbcL5CGMITWgZ1FXJJLrosWev0m6psUfhyHiNQ82v9I9Cq/zqx8fg1bZzEqehUZQvJF3mVYwT84TaGQuZljSuVa6qy+fDFs6UhdJ0+kkex7HpJt1BMa4qI/6eJailG7ZArSdli8eVlCGj6EfReRePWP8ZNpooYhTdV2aLw5WsIeknDrssLPU0tn9OcT2EamapsUfjymjWg8HOaqgtIRqAuQnv25TqDIsfqUpUtFl/5gMtYLTUA31RgDRbzAZcAWE4mC98xwmSA2HdBNEZ4F/rTXJ3xoFnRRsxGR38W9UBLs3B40jjlRTN+SzjZAy5NviNg/gFlSPfFfyi95fo+pKuz61vJLozBk2Yr3q/Opyib5I/A1x9kK93X/kxL4U4OS9V7cnKDHgJZshDu18glmbSHojUStR5JTI1jYH6MgGrmy/MICFvbU9gU5jm90FnsNygP6HZT+rHcQw69d/HY8mvCSeyjwEsQzMkn2BLzGZlpS1Qvdsppcelp4ZvGj1oC6wrwbWCoKy1BNMw7twiAcQud2loz5StyqoUYvEUIibl0B6G2+j+vErez0yg9+UMgO0vVeIBaKDkfQR83bg0RjIJ55yVmwNF4QG6Mh0CWnDv8mqb4U5I+peYM+iZWTMx9FpTsZtOoDoE46L6WBOU/gzQlPUIRpay1RP0FomHe2RAAZ9M0yrYW/UEUvyRNkXVGQ/3qIC4m5p0NwSChprJN6WM8x1AqF2UfQreQBtJNp/WEIV0UzH0NQYdXHgow6BpSbXmF5drp5JCAdS7b1Dw1w3dWNvwnCkM/QEO6KJh3HiwaJGyf9R3hDyG1BCmMrA17XpcvePon6B/oI0hGcH9L80yS6C4a5p3XEQClWIFEdNWWlX9QsffHfE1VthB8kUe5juDbNZBXdlNAIBvCFLQYQIZsCAFAnEIW2RCmoMUAMmRDCADiFLLIhjAFLQaQIZQhFDtiAfhJMYtUZQvKl5chsOqWv+DSs+nGwjwfcOlZkSkVZ7YE8gGXlJQyNC9eXcPQzOfywyHgvOlEM7L1QMs2tthziHrgBR41gJrMAZc+MXc2hNhK3GZELvfh8Y50kzng0ifmzoZQVQSWqh3L/AWXKiiR/bEx38kQDDP2/w/mTdm5r9e5gnxNpJn3EGHk0VmG/AUXmqbWE0+AcwKtvH+gOKg8MYXf+12HdTzEjod3HSTJX3ABhK6u2lKoz5arxi1j0v/NX3DpqiNqafWImLJRiyA3xj/XKpZw6R7UMshJtuQ+QhIK853GCAUcu/10+prIbkWET43ybQs2xsM7nTC/Fx7GZY5m4Ji/4BIL4JZ8fTCPYggwlL/g0qKomFG+mAc3BBjSjCJ/wSWm1ht5h8C8agivyXDRIH3Jtc0VA6nmDZ4tFpqIn8IXXIqDJFZGZBv68I4X5vC/ol9kK/WrtX8VoKa8zdVOLZHZnEQCSLVe/h+hl4xciw9jcl//l7DCMPej7jFQZhRnlK/vX99Bkqv3WUMozMmn6K7bgJL+/gPjpkZCh5wrBAAAAABJRU5ErkJggg==\n",
      "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⎤\n",
       "⎢                ⎥\n",
       "⎢5/6   -5/3  5/6 ⎥\n",
       "⎢                ⎥\n",
       "⎣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": "iVBORw0KGgoAAAANSUhEUgAAAI8AAACYCAYAAADdsLqwAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAATtklEQVR4nO2de3Bc1X2Av7Mr7UMP63Ul27JsMH4VQwK4YEwoxbQFvB4D42BioG54FM+E4iHBZIoJThMIEIqJIBNK4wlJoIXGjAdDi+01Q9w2gMeY4W2aGGEL23rZ0kpa67l67D39Y72ytFrt7t17V3vv1X4zDF7t3nPP2d+395577rm/I6SU5MiRDo5sVyCHdcmb6A3hD8wEigChoTwJ9Emf0qS3YplG+APVQCHa2hcPCbRJnxLUX6vMIfyBfKAGyNe4qQROSp/SNa7M2NOW8AcuBB4D5qdZT4AG4GHpU97RUUZGEP7AUuAnwNkGFqsCB4D7pU9pN7BcQxD+wD3A7cC0NIsIA/uA70ufcmqk3NHyCH/ADbyrYyejGQSWm+nLFP5AEfAOUJChXbwrfcrfZ6jstBD+wErgaYOKe0v6lA3RF7F9nsswRhwAF7DcoLKM4goyJw7AZcIfMOr7M4oVBpZ1pfAHPNEXsfJUGrgjgCqDy9OL0e2LxQmUZ3gfWjGyzS6gNPoiVh6jr770dkaNZjKuLs12BWt0DBzj/pES3Z0O7lk+h9VzFnD4U5fBlTI3dm27jnZpk8dTqPKTV5pYek23pu3sgF3brqNd2uTJd0H59LDWndgCu7ZdR7vMdn7OYSFy8uRIm5w8OdJmwntbE7Jp9SyOHfLQXO/i2nVBVt057p6HbbFr29Nsl3Z5nnjN9Dc9M4Zd255mu3KnrRxpo/3IY0c2rpjNoQ+9cd9bcGE/P3+rYZJrZAmSy7OycqGmEne31aVbmaxRuye5HIm+B6u12aCYJpfnlbrDbFpdQ3O9my1vHGP+BYN8vt/D1oeqcOZLyqcP8+DzLeTbZ8Q+Lrvb6vhgr5dttRWoquC6uzq56saebFcrLQyKafI+T7zh6+lzhtiys4Fn3mygqmaIt18v0t0gsxPqE+x4rpzHdzRS62+wrDhgWEyTyxNv+LpyVhhPQWQWWV6+RNi0333y+BxOtVcAcHCfF5dHZfNNNfxwbTWBZmeWa5c+BsVUX9Sbv8rj03cKueL6HgZDbtpPzGAw5NZVplkYDLkZGvTQe6oUKaGj1cnJ4y4e3d7IinWneOExJdtV1ISU0NddSPuJGajhiadpjI5pEtK/2uoJOnjyOzO5/9nIubG/10l/zzRCvcXk5Q9SVNaO223dQ1JPMDLpSUpBf28hRSUqi5b043LDJVf3sf0XZpv0FZ+BUD4DfWX0BMuR8nQ8qk7G/WxsTJOQXnCHh+DRO2Zy6/3tnL14CABPQR/OvGGkFAwNugm2zuDXD/9UCFGW1j6yiaoK+nsj00mldNATLGPx0hCNR1xIFb74yE3V7KEs1zIpQoiVHHjzBro7FVTViQQKik/hcIx/WC9eTJOQ2pEndvja4YQjn3nZVutgW20FvtuCXH1LN0Wl7XS1VyGlQEoHea5uwPRf8jj6uooZ/fUOhgoony5YtqKH+3yzEcDGZ09kq3oaaEIQOXpCZE5hcWkHkHpME5CaPPGGr1feNv7+R+G0Lrraz8xbHgqVSSmtd1XizB/G5QoxOBAZOHR7+hAOlTUbgqzZYOrns2I4iqqe6djnu/vJcw0Dqcc0Acb2SRwOSWFJB/muASpmNBIOe4QQ8c+vZsZb2EfV7MjAodvbS2VNI06nmuVaaUIIUQJERFeqj+NwDjOt3NDHoGLl0X+KKVXamT7nGN6iPq644adAlYkEmoxTaNZP02PEuXzV7/AUhKieW4+noN+A4gej/4iV548GFH6Gr1++C1iGeQQytn3jCQLNGd5HQsaIAwXk5R80sPiTwMjRa4w80qccIvJEpRF8BrwnpTyAeQT6BHg/g+X/RvqUrM1zjhVHStkPbAOMmnf0vPQpI5cS8Z5VdwNrgb8gvadHe4H9wDbpU0Y6y0KIS4H3gFYp5fQ0yjUE4Q94OdO+xEPwf3htLcWlLSy56u0En5JEfpF+6VPeNKyiGplAnMh7/sBc4GZgMdoTHajACWCX9Cl7x+xzMvPzmEWgVBFCSGCPlNKX7bokIpE4mWRSR4BNdgqzBdkSB7IwkzAnkHFkUxzI0jTUnED6ybY4kMU5zDmB0scM4kCWJ8DnBNKOWcQBEzw9kRModcwkDphAHsgJlApmEwdMIg/kBEqEGcUBE8kDOYHiYVZxwGTyQE6g0ZhZHDChPJATCMwvDphUHpjaAllBHDCxPDA1BbKKOGByeWBqCWQlccAC8sDUEMhq4oBF5AF7C2RFccBC8oA9BbKqOGAxecBeAllZHLCgPGAPgawuDlhUHrC2QHYQBywsD1hTILuIAxaXB6wlkJ3EARvIA9YQyG7igE3kAXMLZEdxwEbygDkFsqs4YDN5wFwC2VkcsKE8YA6B7C4O2FQeyK5AU0EcSJLoQPgDi4HbgfOJLIs8EUPAF8C/S5/yoZEV1EtscgXhD5QC64FvEMmSMXFa2Y//8NcUTmtn4UWfJNiFCrQCfuA/WFk5DROLI/yBNcB1wEwmPnhIItlO9gG/kj4lbiq9CeUR/sA8YDtQqKFug8C3pU/5WMM2GWdEIIezlZ0n9gPnprRh4+GFkbRys1JbUqgn+BLfWvDI6VdmFOdu4HsaNzsErJY+ZVxavUSnrW+iTRyIHJ3Watwm44ycwhZcUEVz/aqM7EQNO2g5+gjOPDCnOAJYl8amfwZcHO+NRPIsSGNHANpWVJkkpJQHWPfAj1FVJ8318wwtXA07aP5qPt4ieGrXQrOJc5oyIN2s9XFjmkiedLPDm3cNrz//q2NUzjpuqEBRcQCqz/mSRUvMuny2nrjEzSamrcCm+jzuu/YsauZFMmJufrHZcmuNu70hKmcdp61pDs3186g+50hK28Vre6kix4gTL7O6mdEZT+02nntxPw//LqsZP3WTrkCj2x57xLGaOFF0xFP7OE/dx16+d81stj6kIC2V13osUYG0nMKibf/lDyppqre+OKArntrkqawO8+sP6nl6TwOnAk7+51VrL9KmRaBo23+2q4mOE2V89L/WF0dnPLXJ4/JICooiC3ldvqqH+oPWX1srVYFcHonHK2g5Np+LlkNrQ4elxQHd8dQmT2/XmdHYg/u9VJ+T9VT5hpCKQN2dzpE+ztE/dTLrnMG4n7MSOuOprcP8ydsFvPTPCm6vSlXNEHc9HNC0vZlJ1IlWww7e3TmP/9wK7oJ+ptc4Wf+IUVnVs4fOeGqT5/JVvVy+qjfhZ3q7CoSofAqokVLerKn8bBMrEETWqmr+aj5fuwyu/VvL9XGEENXAy8w9bw+/2OvAMWr1nlTimQBjBvQi61cW0XOqnKN/mg/8JWCt8Z8oowWCyEJtYOXO8VnApXS2LqXlqBd3QQ/FpZ24vSG9BRsjT6ClmoG+SE89crnnAPJPp983D39zM3zru+P//tQ/QN0E93LnfR0e2ArN9andrvmnm4+IE8d0VDJDSBn5L9RbTKi3mJr5dXqL1CbPysqJ71u9/H9B+rpLcDhUILq63DfSr1oGWH7jSpTqb4/7+xOvj30tpYP2lpqR18Khsv7SiS8udhw7PvLvJcs3svtFMw2iLgGeBLwIAQ5HmMKSTiBxPHe3JZVLmzy72+r4YK+XbbUVqKrgurs6xyxOX6K04XBG16UsllLu11R+hhH+wCIg8eF69MgxQL47xNCAh5/tDvOvm4Zw5kkcTsmmX7VQOWv8qXnDUx/LXS/UG133dBFCtAPDuD27qJh5AW5vP+L0RVayeCZBmzyhPsGO58p5fEcjrjhDAg6HZObZTVLK64VIYVV3sxF7y6G5fgEORzhyGR+ewz9uhZr5R9j5m2ns+m0Jt2/uyHKNkyKlrBNClPHCxwqxa6kli2cStAX44D4vLo/K5ptq+OHaagLNzok+KqXF7l0kulfl9oaYPuc4EBkH6u9xcNa5lhnnmTAWGuIZD23ydLQ6OXncxaPbG1mx7hQvPJbu/BBzkcpNTrc3RDDQwqN3ONn9YhWLlui+Wsk6OuOpTZ6iEpVFS/pxueGSq/to/DLRvGZroOXu+PnLunl6z3FuWA//9vjcyapixtAZT23yLF4aovGIC6nCFx+5qZpt7dsTWsQZHIj83+0NUVHdissjDJ+RONnojKe2DnNZVZhlK3q4zzcbAWx89oSm7c2E1vk4dR95eP7HlTgckO9Wube2CVWdpWk+kNnQGU/tg4RrNgRZsyHuoxiWIZ2JXOdfFuKZNxvG/G2gX/uEMrOhI57Wu5zWi5EzANOZUGYjEsmT7qW2ee9pdba6DZ86Ol4gs7Zfz9BJ3DYlkifdIfaWNLfLKEKIErbc/UvA+JucUYGGBp188+z3DSvXWIJAX5rbxnUhkTw7iTx2qpX/SmObjDLy7Pgf34fque9n5O642xsi372dUG/Ws3PEQ/qUYSKPRGslSOSx43FMKI/0Ke8DDwKp3iJuBh6TPmWP5uplkDFJBwZDBTictwMfoO8wHks/8AbVc+/EJOldJuAR4FUiz6EnQyXyPd0ufUrchxgTJjoY+ZA/UEySRAfSp5huZl2ibBXCH/AABSRKdLB2YYB8115e+jzRI9Qq0CV9yki/IDa5go4mZAThDziBYs7MfohFAn3SpyQcRU9JHitiRJqT0/OR9kgpfWlsa2qBjMCWl+pmyI9jhgRTmcZ28phBnCh2F8hW8phJnCh2Fsg28phRnCh2FcgW8phZnCh2FMjy8lhBnCh2E8jS8lhJnCh2Esiy8lhRnCh2EciS8lhZnCh2EMhy8thBnChWF8hS8thJnChWFsgy8thRnChWFcgS8thZnChWFMj08kwFcaJYTSBTyzOVxIliJYFMK89UFCeKVQQypTxTWZwoVhDIdPLkxDmD2QUylTw5ccZjZoFMI09OnIkxq0CmkCcnTnLMKFDW5cmJkzpmEyir8uTE0Y6ZBMqaPDlx0scsAmVFnpw4+jGDQJMuT04c48i2QJMqT04c48mmQHHTygl/4GvAFUQehp84EcB4JJEcMPuBD6VPGXkQ3kziCH/gIiJLGxSRqH033QslFecIf2BTguJU4CTwe+lTmgytaIpIKQ8IIZYB7wkhTsY+Gy/8AS9wJXAeE6xUnACVSM6l30ufMib30rhEB8If+D6wXuMO4vGq9Ck/ANOJsxn4u5Q+3Hh4IW5vL5WzUpFiGLhX+pS9euqnh3jJFYQ/UAq8BKS28MrEDAJ3S5/ybvQPY05bwh+YCdylcydRbhT+wHkmE2cuqYqjnTwi+YyyxgSnsFvRLw5EUuw8MPoPsX2ei9F2mkpMsO1KTCLOaZZmuPzZp3+AWSOOQJcYWPxC4Q+URF/EylOgu/jermI6WysJhx28vOVfouWaQBwA7yTswzMJ+0jIaIHYt/MWwsNOAs0zGRrQ2t+Jx4gjxizWFkVKOBWoQkoHvV1l0Z2ZRJzUGRrMp6+7GICB/kK6OsooLuscWWooMcYduXUw0okeHnqPlqORNL/CIamYoTfx+kj7tF2qd3c6uGf5HFbPWcDhT8enmevvKURKgZSRHVTM+BAY0FfXLNDfU0R355lFPLo6KulsdSdsuzmpx+k8syRAqLeYcHhszJPFNAHa5PEUqvzklSaWXtMd9/3uzgqkPFNm+4klQFb7AGlRMK2L0Utc5uUPUlw2kLDt5uQ6VHXs2aX3VMmY18limgBt8uS7oHx6/CTVgwMuhgY9CCERQlI4Lcj1638kpczK2Icu8vLCuDyRnMVCSIpKOhK23bz8lnMv+W/c3l6EkEgp6AmWM3p4Rke7jOvzhIfycOYNUVTaQeG0LhwOSVlVu2HlTzZFpZ0jKxsXTLPS0WYEKaUU/kAAaGJ4KI/uYBmh3iLUsBNnnu4fgnHyeIv68BZ9ZVh52cZT0IdwqLi9vRZdEnssefnDlFW2QWWbYUUaVZDtEAKq51pzJZtJQrs8m1bP4tghD831Lq5dF2TVnaZL3p0x7Nr2NNulXZ4nXrNeB9go7Nr2NNuV9TnMOaxLrs8DsHHFbA59GP/WxYIL+/n5Ww1x35viJJdnZeVCTSXubqtLtzJZo3ZPcjkSfQ9Wa7NBMU0uzyt1h9m0uobmejdb3jjG/AsG+Xy/h60PVeHMl5RPH+bB51vIt8qIfZrsbqvjg71ettVWoKqC6+7q5Kobe7JdrbQwKKbJ+zzxhq+nzxliy84GnnmzgaqaId5+vUh3g8xOqE+w47lyHt/RSK2/wbLigGExTX7kiTd8XTnrzOu8fImYAv3ug/u8uDwqm2+qwe1V+e7TJ1GqrXa7IoJBMdUX9eav8vj0nUKuuN66v8JU6Wh1cvK4i0e3N7Ji3SleeExJvpEF0RDT9OXpCTp48jszuf9Z+/d3AIpKVBYt6cflhkuu7qPxS/s1WmNM05NneAgevWMmt97fztmLh5JvYAMWLw3ReMSFVOGLj9xUzbZXu9OIaWrjPLHD1w4nHPnMy7ZaB9tqK/DdFuTqWyx55zllyqrCLFvRw32+2Qhg47N6Z+RlFwNimpo88YavV95mj/s6WlizIciaDcHkH7QABsQ09rRl/akHiZmM9k2Z7zBWHqNPPWY7Ok3GVaHZTt9G1keOLi9Wnv2AkWMX7xhYlhG8S2aPDIekTwlksPx0MDIGn0ifMvIDHCOP9ClB4EfoF0gFaqVPOaazHEORPuUk8DjG/kCiBIHNGShXL69gjEDtRNwYYdyz6jDyfPMyYBrpJTp4z4S/wBGEP1BOpH2JEx2kRjTRwQHpU0z7mJHwBxaRXqKDMHCCSPvGXMLHlSdHjlSYAjelcmSKnDw50ub/AXdP1nrcWkF9AAAAAElFTkSuQmCC\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.Matrix([[1, -2, 1], [10, -20, 10], [1, -2, 1]]) / 12\n",
    "assert expected_result[:] == isotropic_2d_00_res.as_array()[:]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "sympy.tensor.array.dense_ndim_array.MutableDenseNDimArray"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "type(isotropic_2d_00_res.as_array())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "sympy.matrices.dense.MutableDenseMatrix"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "type(expected_result)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Isotropic laplacian"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAHoAAABNCAYAAABpJnDxAAAACXBIWXMAAA7EAAAOxAGVKw4bAAAINklEQVR4Ae2dvY4cRRSFdy3HaDESIiCaJSBey0/g5Q0WExBjvwEWEaG15gkWiYAUWyIh8/IExg7IcYyQzFqIgGz5Tk9Vq7qnd6aqu3r61kxf6U7/d997Tv3X7Z7D6+vrg1Q5PDw84Zpzrv0s9dopz8fuc55/hC7QN+hjfHjH0rwMxfx2iofuYY/cNQKrGHEkX0CsCD5g+xmLX9G72rYqOTFX7lS27tJnyvFt5dxT9I/2fsvbshddeBtZ937X+/wxi0vsXYu586+LQ9y5Pghz9PfuZBa1VKm/3ip/xRfZ8sQX2eG+kj301VLowz02zrQjJFp17q4RWzuNb8f1xnJFOVry23JR9i/+KaM2hGJfJFdE32oc2a8N5YBiGmNDqdlLol3D7JJc8HQogKVcv3dEQ/JDyDmCZN97KIWrQXbuFdGQrJbrsSdZ26ivqwcBaf3ivkRr0KEocYSqXn7J+pmU9cdoKQ3QQZiHre6NxAHOwoGjnLFgW4MOLwup6zQ4IrBkcy3Y7rtZ9T5LK7kwTyIaUJT6i6zbsP19SwTG2pIL875Fd6yd83lGEJiJNkLE2GZEF93UFRunuShmDsc2uM/9S7U9p93RRFslMYb4Um3PaXc00SGgpDR1UyRv0XcYtDLOWh019uPsls0a976DfoXtplvdHsIcmCdN1/HgV+gpAMkGdbeutG5dsfMCfejtdNsv/LblJbb2wpzrNFaAa9cHSY0xUtXXXPiGCy9ZSv5Gn1RrZfyEM1gVeNbNzol5dI4GlCv0zHLqj7UNP5TDzefoIZiLKxRImoEH7Nso1TCcS2U6WTmkuKk+7JcfD9D7qHXJgnl0YwxwVB9L7pFCNEasuKtTFioCwyJRh8yKS6RfYOAT/Hht1lAMy4l5Sh2t+ljyYrmoygTV1RrzFuFFCOQ+RRUQeIzdtS9Gjc+GeTTRgLOuG+Jzu1G8Os1SF1HTlFWoTecZE+/MiXk00c5nFXVdpJqe6oPMI/QK7Zp7Vn/asmTBPLqOdkiobv4BwFT0ecJ/J+X57pZJwJQzsFnFYFgqqbr5F/0SXRnw4Xx1JT9A32rJPap2CevbliyYp+ZovZnxnfNUKe2TbXs94Hmfc+0jEYj+yPo36C/ox2hDOK5wo4pcCFZcmYIVGvPYjQvG3ciGeUo/uuggeEhrjODBj3L1yosI2qdj4flsayhQsWaNe4y97WypXzJgO4Wvuh+dmqN5Tl1ka90Xhb4Y176ihZyrfqv8UVEfinxVwphCQnx7YX47xWpSb7u/7Bs3OxEE77DwoHpAPUQifusNt1yY3/Je9FzuchB8O0cLIuX2qaUX5r2JpojTA3c5CL6de9vbWyd8COa9iOaBuxwE78cE2rlX2+ppTCJDMU8mmgeqQbKzQfDUiaqbRXZXDp6kLZID8ySieeAJAKjILjUIHtMb0s61/qCmMNXvrgS/1U157hKB27udRS7MFcwn8qoZKBzxxVanFzxU89Er4HCdyaDATifYiR8LFhpxUumk9edo40UEzvEjYxya7l2tIZhzrRKoPmZwmNq9KjIIXkyF4hL02hcROMfEm5bYkQXzpKI7BGteLwuBmeiy+OptrYru99C/3PLGG1HeTx7Ajw1qH+hluRRRSG/dLYrxY9PNud9W2iQxtm6w5SN8EbfVN0z+YfkhquWNsuGGN16X8wA2qOsz6HNRFvyIxSSDrX/yLHHb+FhN1PNJZepeaY5W49530GKC4LF1rVj1LYddSXU0D1T/UtN6irtSq1XjwVPN064lLfWgVd9y2ZVEtAMvnMFS/1t90V0Rq74Ntiu1H93ue6q+NB1GFJsCXQkVnm7Ct1x2JREdokCRohbwA/R+uH8X1q36NsSuPkW3hhA1PKhujvkg+NSEZ9W3oXb1Ito1xlS0lRAEn8S1Vd+G2tWL6AA5dbVMB8EHtqauWvWtl13RRKt+QEsNgl9LslXfctoVTTRFh0al1G9WMPwJqveW1LXS/p9Q84LN5+iFbHdLNSgVvlv7FjgxuW+hXdg7CPPUVve3APEz+h/6KSqA7jqDWLUrAKUib90X+H2Af/VmBucuUAu+ZcNcgQeasKiDxCFubZA65yq1rwS+b7puyuOyN/SR9WS/J7Y/GXN8HBTAz/XFinKpF5VGknDfcs8O/qYW3cVCQG4MhxHlh3K0ZJKAv+Wjt/cb3Rjbnklbe5Lq7OI+y9EXnb0k2jXMdvnlg5X0sHdEQ/Iuv3ywQrDfsVdEQ7Jarjv78oEntWvZl+hqoKHrhlb3QbIaX6qXS335YBDmSa1uwFoAVB34zraiSxqB72xbFc22CaxGRAytcd/NMml3LsyTiAYUvcnRDj4wCVDbKGzPEgjfvu/Y27kw71t0j+3ffP/MCMxEZwbU6u2ii27qiskD+PuCWKrtOe2OJpq6YitvJ/Qlc911pdqe0+5ookMgSWnqpkg0pTd/gb+CYtyfHJgnTdfhThXLTWqTZ+puzV/g3zCtK6yGaF/Mua7fNCWpStGf8xf4AWFbkhPz6ByNc/MX+AfmztScPQRzrq1zdGodXQ3DuVSmRK053uKm+rBffpTy8kEWzKOJBhzVx5L5C/xLHEb/zYl5yoCJIkAl9VfrKYbmL/AvMRnrNxvm0URD6rrBf5/bx3J4jPuqi2j65YOcmEcT7ZB+zbKL1LWfrRqDpZR7qk5GS335IBvmvtWtjnFb9Y2qug/I8VO0/q8o1tWqexWeY3UdO9vhvoo0qb6bZtVm2YVEYe78a/NXbes+aowpN970NwJKTbVwwaXLHfrygYpytQhLeW3WaoB+jW/XSgLmqoqqFnrXff4H1m8Da4MfCKsAAAAASUVORK5CYII=\n",
      "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⎤\n",
       "⎢               ⎥\n",
       "⎢2/3  -10/3  2/3⎥\n",
       "⎢               ⎥\n",
       "⎣1/6   2/3   1/6⎦"
      ]
     },
     "execution_count": 11,
     "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": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "expected_result = sp.Matrix([[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": 13,
   "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": 14,
   "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": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWwAAAFmCAYAAABJBKDfAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAS20lEQVR4nO3da3Bc5X3H8d8jW5ItW7JsycT3G5bSGAYYoGQmAZtbKALiwZca2uLiCTATIC2QdkKCmwFaUwghdQIFkmmJS6GUYqCZhGFpS6Eww5QZsAtOWzARYOqYDLCSdbNulvbpi7VkXVbalbR7zvmf/X5mmNGuzu7zvDj+cvY556yc914AgOgrCXsCAIDcTA97Aig+LpFcJWmFCrf/eUldkvb6htrOAo2BIuYSyUWSfkv524e9pCZJb/uG2v4xx2VJBEFxiWSlpIcl/XZAQ/ZIuts31P5jQOMh5lwiWSrpB5J+p0BDJCXd6Btq38r0S5ZEEKSbFVysJalc0h0ukVwW4JiIt6tUuFhLUq2knS6RdJl+SbARpHNDGnddSOMifs4NYIxFkuoz/YJgI0jVRTYu4ieofSnjOAQbQcr4MS/G4yJ+gtqXWBJBxDW+XaZvrV+iLatP1CXzM34kHGb/nnJdf85yXb6kTtefs1z795QHMEtguN0PVOuGtcu1flGd7r52wbjbPvlXc3Xl50/UxuWrdc91C9TbPaH/ARBsRMf0Uq+zv9quG+/9JOu2vT3Sjm2LtW5Dm55qbNR5m1u1Y9ti9fYEMFFgiJqFfdpyU5PWbWwbd7v/fL5CP/vJPO3YfVC79n6gTw+W6pE7ayYyFMFGdKxYc1Trr2vVypOyV3fvSxXq75euuPmwymZ4bbmpRd5Lb75YEcBMgePO39yhczd1qHLumNdPS5JefHKOztvcqtWn9GpOTUpXfrNJrzw7ZyJDEWzYdOCdci2t75Ebsgsvre/RgXdYFkE0/bqxTKtOPn4wUndaj9qap6nls5w7TLBhU9eRElVUpoY9V1GZUlcH+zSiqbuzRLOqjh+Fz65O/3yknWDDgBceq9SGZXXasKxOt16+eEKvnTlrdJy7Oko0c3ZqjFcA4ZpRkdKR9mmDj4+0pvffWZU577N8lwjCc/HWdl28tX1Sr13xhR794pG58ikNLoscfK9cl33tcB5nCOTPktW9+uC/yyWl9/lfvTVDVfP6VT0/52BzhI3o8Cmpp8vpaG/6UqeeLjfmZU+nn9+pkhLpqfur1dvttPuB9I0GZ17Ilz0hWH1H0/tqql9K9ad/7js6ersLrmjVy7vn6P1flqmtuURP7pyndRtbJzIUR9iIjo8PTNd1X1w1+HjDsjrVLOzTY/s+kCTdevlirTmrS1ff1qyycmn7rkP60S0L9MR987VoZa+27zqkMs45ImCP3lWjZx48fnnea89VadONTbpkW6tuWLtSD736oRau6NOXLu3UR+82a/vmpertdjrrog5dc3vTRIbi2/oQGJdI7pU0K4ShH/QNtfeHMC5ixiWSP5f0+QCGuto31L4+8kmWRADACIKNII1/Y0Hh9IU0LuInqH0p4zgEG0H6VUjjNoY0LuIniH2pX9KHmX5BsBGkx0MY8/8kvRrCuIinf1Lhj7Kf9w21GU9GctIRgXKJ5FckbZG0UtK0LJtPxRFJr0v6G99Qm/3LpIAcuUTyLElblT75WJqnt/WSmiX9u6S/9Q21Ga4LJNgAYAZLIgBgBMFG5DnnHnfOeefcmrDnAoSJYCPSnHNVkjYqfeb8GyFPBwgVwUbU/YGklNInKLc652aGPB8gNAQbkeWcc5L+RMNvZ98U0nRQZJxzpc65jc652WHPZQDBRpSdIWnoHzWdLelPQ5oLioRz7gTn3B2SPpX0jKSvhDuj4/i2PkTZzZJGLoHUO+fWeO//N4wJIb6cc2dI+raky449NUPSYUnPhTapETjCRiQNOdk4ch8tlfRHwc8IcXRs2eNK59w+pe+I3ah0qGdI6pb0Q+99xptYwsCNM4gk59zXJd2nzF/HekTSfO99V7CzQlw4506QdIOkm5Reaci0Tt0taZn3/rMg5zYejrAROWOcbBzKi5OPmATn3BnOud2SPpJ0q6RqZY51StLPoxRriSNsRJBz7kxJ/6Hx/9jB297704KZEaxzzm2SdLukE5Ve7sh2sNopaa33fk+h5zYRBBuR45x7XNLvafx/VF2SzuTkI7I59omtWemj6Vzt896fWqApTRpLIoiUcU42jsTJR+TEp49KL5LUkeNL2iX9ZeFmNHkEG1Hz+0qvH2YzXdz5iBx579+QdL5yi3afpGcLO6PJIdiIjBxONmayuUDTQczkGO3IXco3FMFGlJRJqlT6CGfofwNGPl8uqS7gOcKwY9F+IctmDwcxl8ngTkdEhve+R8NvRZckOedulrTTe5+vv+6BIuWc+wulP5U9LeliDb+kL5KX8g3FETaAonAs1n8maYf3/nc1enmkW9K9YcwtVwQbQOyNiPV3pYxr2o1Ru+56JIININYyxXrAkGi3SLozhOlNCGvYAGJrvFgP8N6/4Zyb5w3cRcgRNoBYyiXWAyzEWiLYAGJoIrG2hGADiJW4xloi2ABiJM6xlgg2gJiIe6wlgg0gBooh1hLBBmBcscRaItgADCumWEsEG4BRxRZriWADMKgYYy0RbADGFGusJYINwJBijrVEsAEYUeyxlgg2AAOIdRrBBhBpxPo4gg0gsoj1cAQbQCQR69EINoDIIdaZEWwAkUKsx0awAUQGsR4fwQYQCcQ6O4INIHTEOjcEG0CoiHXuCDaA0BDriSHYAEJBrCeOYAMIHLGeHIINIFDEevIINoDAEOupIdgAAkGsp45gAyg4Yp0fBBtAQRHr/CHYAAqGWOcXwQZQEMQ6/wg2gLwj1oVBsAHkFbEuHIINIG+IdWERbAB5QawLj2ADmDJiHQyCDWBKiHVwCDaASSPWwSLYACaFWAePYAOYMGIdDoINYEKIdXgINoCcEetwEWwAOSHW4SPYALIi1tFAsAGMi1hHB8EGMCZiHS0EG0BGxDp6CDaAUYh1NBFsAMMQ6+gi2AAGEetoI9gAJBFrCwg2AGJtBMEGihyxtoNgA0WMWNtCsIEiRaztIdhAESLWNhFsoMgQa7sINlBEiLVtBBsoEsTaPoINFAFiHQ8EG4g5Yh0fBBuIMWIdLwQbiCliHT8EG4ghYh1PBBuIGWIdXwQbiBFiHW8EG4gJYh1/BBuIAWJdHAg2YByxLh4EGzCMWBcXgg0YRayLD8EGDCLWxYlgA8YQ6+JFsAFDiHVxI9iAEcQaBBswgFhDIthA5BFrDCDYQIQRawxFsIGIItYYiWADEUSskQnBBiKGWGMsBBuIEGKN8RBsICKINbIh2EAEEGvkgmADISPWyBXBBkJErDERBBsICbHGRBFsIATEGpNBsIGAEWtMFsEGAkSsMRUEGwgIscZUEWwgAMQa+UCwgQIj1sgXgg0UELFGPhFsoECINfKNYAMFQKxRCAQbyDNijUIh2EAeEWsUEsEG8oRYo9AINpAHxBpBINjAFBFrBIVgA1NArBEkgg1MErFG0Ag2MAnEGmEg2MAEEWuEhWADE0CsESaCDeSIWCNsBBvIAbFGFBBsIAtijagg2MA4iDWihGADYyDWiBqCDWRArBFFBBsYgVgjqgg2MASxRpQRbOAYYo2oI9iAiDVsINgoesQaVhBsFDViDUsINooWsYY1BBtFiVjDIoKNokOsYRXBRlEh1rCMYKNoEGtYR7BRFIg14oBgI/aINeKCYCPWiDXihGAjtog14oZgI5aINeKIYCN2iDXiimAjVog14mx62BNA8XKJ5CpJJ0kqG3fDq287RY375BLJTeNu9+hdV+rUszfrg/+5x7c1E2vklUskSySdJmmZpGkFGsZLOiLpdd9Q2zpqDt77Ao0LZOYSyVJJP5R0YU4vaGueq7bm+Vqy+r0xt2lJ1qijpUazq5tUXXtI0i2+ofalvEwYRc8lkgsk/Z2klQEN2S9ph2+ofWLokyyJIAyblWusczE81k2SZkj6vkskxz9yB3J3q4KLtZQ+gv+uSyRPGPokwUYY1ubtnUbHesBspT++AvlwTghjlkg6e+QTQNAq8/IuY8d6QFVexkHsOOf+3Dn3h865mTm+JD/77MQN24cJNsLgpvwO2WOdn3EQV9+R9KCkz5xzP3bOnTTWhi6RDHM/GjY2wUY07H6gWjesXa71i+p097ULxt32sXsW6+tfrtEfX+D14++UqrebMGMyZkuaJekaSW845952zm2dwFH3cY1vl+lb65doy+oTdcn8+qzb799TruvPWa7Ll9Tp+nOWa/+e8lyGIdiIhpqFfdpyU5PWbWwbd7uXn16g5346S7f9tEW79r6vTw+W6pE7awKaJeJpuqSZkk6R9JDSR90PO+fW5P4OpV5nf7VdN977SdZte3ukHdsWa92GNj3V2KjzNrdqx7bF6u3J+lKCjWg4f3OHzt3Uocq5/WNu05Ks0SvPVmnthi6d8uVPNacmpSu/2aRXnp0T4EwRbwNH3ddKetM595Yu/dxVSqXG/xS3Ys1Rrb+uVStPyl7dvS9VqL9fuuLmwyqb4bXlphZ5L735YkW2lxJs2NHRUqOPP+xX/Wktg8/VndajtuZpavmMfRn5NHDUfap86iH95sPVav7kc+rtmfqlogfeKdfS+h65Ibvs0voeHXgn67IIdzrClu6uaervW6hfNy6UJPX1pZ8/8O5qjbwx7OFv/6u75JWAJ4gYmi3vpc72OepsnzPuDVy56DpSoorK1LDnKipT6urIetBBsBF9lXMPD/48Y+Zcpfo7VTUv/dGzrdlJqlXtwqSq5g2/bXfVyU/rv16Z2j8uxNXOCW3tXHrfmjGrXZL0wmOV+sn29Mnx+tM79b2fHcr5vWbOGh3nro4SzZydGuMVgwg2os85qWpeOtpL62boYGPf4OP9eypUNa9fS1Y3j3rdNXf8s3/6r/8t0LnCBOfc95W9f/2SuiX3G1XVVGpWZZtKpqXDffHWdl28tX1Sg6/4Qo9+8chc+ZQGl0UOvleuy752ePwXsoaNqOg7KvV0OaX6pVR/+ue+o6O3u+CKVr28e47e/2WZ2ppL9OTOeVq3cdSX5ABTcERSl6QnJK3T85/Wq7K6dTDWmfhUep892ps+OdnT5ca83PT08ztVUiI9dX+1eruddj9QLUk688LObBPjCBvR8OhdNXrmweOX5732XJU23dikS7a16oa1K/XQqx9q4Yo+fenSTn30brO2b16q3m6nsy7q0DW3j3XjDJCrY0fT+ljSDyQ94b1vl3K8cebjA9N13RdXDT7esKxONQv79Ni+DyRJt16+WGvO6tLVtzWrrFzavuuQfnTLAj1x33wtWtmr7bsOqSz7pdh8Wx8C5xLJf5B0ZgBDfcM31LIkglGcc0eVPmA9ovRKw25J93vv94zaNh3sd4Od4aB7fEPtroEHHGEDKEZNkto04mg66gg2wtAX0DgZFsEBSdJy7332m1wk+YZa7xLJfhXujxaMZ9g+zElHhCGoj5dc0oeMco31EPsLMpEJjkuwEYanJHUUeIx/8Q21Hxd4DBSPXdk3ybt9koatqXPSEaFwiWS9pKsknSypNE9v6yW1SHpZ0uO+oZYlEeSNSyTPlbRRhf+bjh2SXpP0976hdtjaOsEGACNYEgEAIwg2ABhBsAHACIINAEYQbAAwgmADgBEEGwCMINgAYATBBgAjCDYAGEGwAcAIgg0ARhBsADCCYAOAEQQbAIwg2ABgBMEGACMINgAYQbABwAiCDQBGEGwAMIJgA4ARBBsAjCDYAGAEwQYAIwg2ABhBsAHACIINAEYQbAAwgmADgBEEGwCMINgAYATBBgAjCDYAGEGwAcAIgg0ARhBsADCCYAOAEQQbAIwg2ABgBMEGACMINgAYQbABwAiCDQBGEGwAMIJgA4ARBBsAjCDYAGAEwQYAIwg2ABhBsAHACIINAEYQbAAwgmADgBEEGwCMINgAYATBBgAjCDYAGEGwAcAIgg0ARhBsADCCYAOAEQQbAIwg2ABgBMEGACMINgAYQbABwAiCDQBGEGwAMIJgA4ARBBsAjCDYAGAEwQYAIwg2ABhBsAHACIINAEYQbAAwgmADgBEEGwCMINgAYATBBgAjCDYAGEGwAcAIgg0ARhBsADCCYAOAEQQbAIwg2ABgBMEGACMINgAYQbABwAiCDQBGEGwAMIJgA4ARBBsAjCDYAGAEwQYAIwg2ABhBsAHACIINAEYQbAAwgmADgBEEGwCMINgAYATBBgAjCDYAGEGwAcAIgg0ARhBsADCCYAOAEQQbAIwg2ABgBMEGACMINgAYQbABwAiCDQBGEGwAMIJgA4ARBBsAjCDYAGAEwQYAIwg2ABhBsAHACIINAEYQbAAwgmADgBEEGwCMINgAYATBBgAjCDYAGEGwAcAIgg0ARhBsADCCYAOAEQQbAIwg2ABgBMEGACMINgAYQbABwAiCDQBGEGwAMIJgA4ARBBsAjCDYAGAEwQYAIwg2ABhBsAHACIINAEYQbAAwgmADgBEEGwCMINgAYATBBgAjCDYAGEGwAcAIgg0ARhBsADCCYAOAEQQbAIz4fxjyRMz2Xz2bAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 1152x432 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": 16,
   "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",
   "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.7.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}