demo_create_method_from_scratch.ipynb 65.6 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "nbsphinx": "hidden"
   },
   "outputs": [],
   "source": [
11
    "from lbmpy.session import *"
12
13
14
15
16
17
18
19
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Demo: Create lbmpy Method from Scratch\n",
    "\n",
20
    "<img src='../img/collision_space.svg' width=\"90%\">\n",
21
22
23
24
25
26
27
28
29
30
31
32
    "\n",
    "\n",
    "### Defining transformation to collision space"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
33
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAoIAAAAVCAYAAADSDI/HAAAACXBIWXMAAA7EAAAOxAGVKw4bAAAI60lEQVR4Ae2d7XHcNhCGaU0KUJQO5A6kuILIHUjuIFYHzk/pX0buwEkFHrsDJxVkog7kDuKog8v7ng4XEgcsPgiCgADM8EgCIHbxYJdckkdyuLm5eYfpYTedbTabwTSh/NSUX0JeCbqVoIM0FiXoV4IOnZHZvyUuqqyE8StBB8XDNC9BvxJ0MLFReSXoV4IOioc+L0G3EnTQuYzXS9CvBB3GTMbLLt1QPon7joZh+AHT9e3t7UtM91g+SMh/h8yzg4JyMk53Oq6iUQV8yKUzcltHZyQz6nxkPiztjDojNwG5RrchmQ9LOyOZkcgHMct7TC/RxGdMJy8QGd5h4Qsy/zC1i/xL5L/C/Be9HHnclukfTGz0DnlfmZEy+chBHQarj5j/llK2qy3IW52P0hG6MFj/hOkcy48qX82R1ywj9P0UHJQN/4jlb1xH/uTkp0RGHD/oJY4t68xNkOH055b5KL6usWiZEfpepZ/56q1sYM4csqrzs5x8FFvIFPd5KC/qeJaLka8cHz6os43/xEAQlY4xKH9ifq4GR82R9zeWf8WcEeWAOesy7zWWkwWDaMtbzq7uT5gfBEHUMXWCnBL4UIffMTGwYYBD5/nexqBRRjw4fUDfX2O+TVimA3BHQnudnASVwgh6BI3tU8/ifnd99vLnRvkEjUWjjGr1syC94zzsaatK/Swnn1r9LAsj2E+QHNd+COXbQPDIYdSs9EGvg43fIu8Y820QyHIsP2LG9YP6LI9JEXIomzrnSqvyYSfJHdMVpmusfvToeHOMwITjRD77BF68Okib5RVUPZXCKHRs9X54rYNFqD83xYcQwSh0LJpjBExV+lmE3jSJ4FSrn6GjoeMazEZtUKufZWQUOhZe+yFXIPgGA2O61XqFjk9uqe0G8i/ML7ANo/oUKUjOTlfqnEq+qw9r83Hpd1DeKKMLgHgw2AWvBPKEhmdZ+1QQo71OCy/U6mcLY4lvvkEbIqxa/SxI73irGGr1s1x8gtEW5Ge5GAXJ8eVjDQTRAP/7ZrvFS2V4K1JPqj7LU6QYOdThTQrhUhuF8JFUlMpaY8SA7yvGjFcATcl04lACI5OuS+TV6mdLsEjZZks2RG61+lmM3jF2Uquf5eITw5TblOBnuRjFyHHy+U4gz/9TUegk4WBqOmhO6mDlRM8IXZ8hhzpTd9OVzFA1pPqr8pEU8yhrihFsiWfipsT/Uw4oN13dXpWRSdkl8mr1syVYLNBmEzakuNXqZ5F6q257zWv2sxx8vCDaK63uZ7kYRcpx8rFeEQRzPnjwYGCvgjzb1RVu4hMsGpqeZMXKoc6TW32TVtOtrM1nTk+aZwSHYhBIO1FPEus812ak67PUeq1+thSPlO22YkNWZrX6mYfe1j5bCp6Vny3Ax4LNK7tIP8vFyEOOk48UCDKYM93+9RkZvpswRzLJoc45AsEa+NjGoDN6ekjkM5zovQVSDYwsqifPrtXPkoMIbLDbUL1+xofIpP1DoCl4Va/Jz9bgY4NYqp/lYuSS4+QjBYI8gzFd9WOjtqTOevhewbkpVg7vh6e4IunSf20+Lv2k8qYZIfjjk1T8z6DtljHZrc1IGr+UZbX6WUoGS7XVig0Z+dXqZ556G/ssZD4bP1uIj4DOWVScn+Vi5CnHyUcKBI30IVgFh6ZgS+VR8Kw0Q44tQJulj+/GM/T2FZGiXrOMMD58VcoJ5vt3ClqArsrIolPy7Bn22gSfmcCbZVSrnwXoHWQaz8XPluITBPOwclF+lotRgBwnHykQ5BmMCux09Pzzoen2KwUysTxFipFDnaWzrxR6sY0S+MT2pUlGcBw+Cc9PKe6vBGKZn+Ix2XIJjGLHN3S7Wv0stJ+567dkQ3u2tfpZoN77/gYsVO1nGfgEoJxULcbPcjEKlOPkIwWCvKpnOkByBHhPmg9L6OkcGfdQ8lEVYJlKxCZvOSMBDEadVyRn6kVxJfAZdTtosTlGGO8zEDJ9KpHBoenEoQRG3oM6055r9bNcfLzlaBVbsqFt15fws5m2rYZE2l8PvnrP1KVaP8vER41V6LwIP8vFyFfOCKKTjxQI3qOhV6PG9otQhK9m+YY5D6LbhGUGfHx/389POdvXcjDvX5TxM3HByVeO1vA2GNXyJqs7XaP12jW2Op9Jp55W1J+NOfBSaooRxpsnNNwJ8+XR/NTcfkLeNdb3Jy4jaKszGunCRevYQv8m/SwXH02OOBZa3WZsiP2GHSb3s7m2PRoPaX/tpfdcXbC913FzpDMXV7ch33Gdy0frN1et+zytbjOMfMcilI8UCH5EY2dag+NVwud3Wu8w8Zut/N4tv/NLh9smLPMAyzMxHoB5sIpJTjlao3xp5xctb7KaSK9S+AzozydO6CT//8a0XUeeWn/K/f+3NUa0B+7syUOfTEEgSRXByGdsUadZP8vIZ/CRRcMZpdZsKLmfJbJtDom0v/bSO5EuNR7PcvKp1c9yMfKSM9oHcdG5Hxpubm7uMF1sNptBn5D/gOlMzw9dRxuXmI5DtwutDxmnmB58t5urF2VhqoYPuUDfzshg62Ob6YwO9wWdj8xkzKf7mR+rED9D3dnHELRR1f46hM/O5mYxqo1P9zO3n7lsCOXb+E+6IshI8uADx8yMSPxvlu3KS0Rz1k34cmDq7Jvm6lUbH3LpjNzW0RnJjDofmQ9LO6O0jObuq6lNbfvrbkNpbYitzbWjZ2lDYiCI4I3/abA9VekeItRAG7wlnOK9gqI8yOGtP+rq9Wm5FHrtZFXBh/Cgb2ckWlFn5MDTbcgFqPuZB6EwP8N+K8kxBO3045kwOjXxYTegbz+eyePpzedIaEcV8VUbfAFvbHqLAbN9vSG2TdN21PHaVGDJS6VXLXyIoTOyGMMouzMawTAsdj4GKFpWZ6QBMayGMEq1r6YateyvQ/iwX6kY1cKHfe6MSMGevPm84D1itKOe/r1C0LZ/2EO1v4u8LzMFdEqs9xx68WEVfg7I+doY70YDKpbOh13pjNwD2hnJjDofmQ9LO6POyE1ArtFtSObD0s5IZuTisyvnhTO+YeTqP+82TfEYFbz0AAAAAElFTkSuQmCC\n",
34
      "text/latex": [
35
       "$\\displaystyle \\left[ \\left( 0, \\  0\\right), \\  \\left( 0, \\  1\\right), \\  \\left( 0, \\  2\\right), \\  \\left( 1, \\  0\\right), \\  \\left( 1, \\  1\\right), \\  \\left( 1, \\  2\\right), \\  \\left( 2, \\  0\\right), \\  \\left( 2, \\  1\\right), \\  \\left( 2, \\  2\\right)\\right]$"
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
      ],
      "text/plain": [
       "[(0, 0), (0, 1), (0, 2), (1, 0), (1, 1), (1, 2), (2, 0), (2, 1), (2, 2)]"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from lbmpy.moments import moment_matrix, moments_up_to_component_order, exponents_to_polynomial_representations\n",
    "moment_exponents = list(moments_up_to_component_order(2, 2))\n",
    "moment_exponents"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
59
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAVoAAAAaCAYAAAD2fpSuAAAACXBIWXMAAA7EAAAOxAGVKw4bAAAHXUlEQVR4Ae2d71UcNxTFFw4FYNIB6SC2O8Ad4KQCkw7s40/wLQd3YFxBgjuwU4ENHdgdBG8H5P4GaSxmd2alGa14y0rnaPV33ru60r7RvNHCzu3t7ayGfgbOzs721frW9Th06SvVz/uvqi0WGbA6lxZxWcSUsqZy458qb68LXgJPVPdd6edu25aWz8XFn37syr9X/krxV19X041hwOpcWsRlEVPKQsuNP0oe9iG0Fx7wrs+QqsOxkqdKq5H9ScyJ+Dj6WZydK3+out+CuprdDAaszqVFXBYxpayy3Phj5WGQP3WBtoZWjc0jstJ299btvKVl+Pi6pWN/bMO2OpcWcVnElLIec+OPkif7+V0gL5W+DsHueB+tGrDCdLgIO9T8fQbEDzvaY6XVdXCfmo0rWZ1Li7gsYkpZcLnxr5Kn9m/Ch3egeZfT+GhV4CXPM6UvhsCrncfly1DAUP/H1ubG37hXHtvYtm08VufSIi6LmFLWa278kfJ4l/NB8SVY/cuwN8ov3clKKC4FLrhRfKbo37wruz1BPDBudrPtXWp7Rv+4Rmp1Li3isogpZTXmxp8gD3uKv3Zfce4NLScNlu5m6aS2xiorj99h614CadwY2TdKG45ceaYUf0wNG8SAmztzc2kRl0VMKUstN/4Ueeo7V7wW3t8VL/ZU8IazvvBZMouOXB4D+HJ6rnCM8xRQwwYxYHUuLeKyiCllqeXGP1IeNpVN6gU72iNFzs02Tlvla7jPAGdmcZ+QtkF8YWxr2CwGrM6lRVwWMaWsttz4x8jjhRg72tnO6ekpL7c4F/qUiqGgPrgO8FM+yW2YJQ9jhguDt/lXKrc+Y9f2QWnjwlD7pFBSVyzQABOXPFd8pYjL4g9Fwhf1+XiXnf6Zoi/ou/a5YWSBPoqDXAR9i2ADUEwIcNF9cAwx8nL0yY0pkFeE+0BfFk4DeWvBL/lsYjnN9WRXHxg4XnQ9dHgrYO8EAmAY8zBwV+Btf65QUlcsZhzn7xwHX3QRLyCPVMZFwRx1OVHVpJCirzRflrHFkp4yhliZU/vlxmR5XcRwtW783q4eYmgPFOcxqNbVR8YE3yfGhcALJw+wqXB1WX6tVlKXB78qFSZ28qEhZT64sfhdPXMUtqs4PqToK82XZWyxjKeMIVbm1H65MVleFzFcFcLv7eoBPlp2S13DFoM1Zx98xNdOILvXvzrC2YJ36zpdoosldcWC+qrxhycYcONcq66ZKKVZXCYBmBR9pfmyjC2gcDCbMoZBQRkbc2OyvC5iaCuBv7Wr7Gj5MrNjerAQGBQMKobf7+RmamO3S12uHa03XmvXJcxRQWP0Nxnfn5vN376QO03Rp75F+bKMLXYeUsYQK3Nqv9yYLK+LGK4K4fd29QZDaymwc2t3cg4YBtGfScuJtaSuaNxaAP4G0L74Uh2HnrnZZA8J+orzZRlb7EQkjCFW5OR+mTFZXhcxXK0Tv//OzjG0PLL6ihhgvX00gVPlHDo8oQ58tgu72ZK6QjB9+bF4uE7xkyIGlsCxMW4soSsBp/28aXUfKo/imusUk/VJ7drnpgS2kMO+/FhukTdhDH1w2vqxuNaJSeAsr4uWu4FMNP4BGX1N93a0vITip7Ux4RfXyQtor3GL4IdSzpuNDaFxmUkWxofISYQ2lNTVKh3ITMTjx3jj5LR+HVSqjnb/orBBUVpfo/TuhuyyLa7cc5PMRQq2FvxAZiK3SB47hgFUdwZcHcZ+v9aCyQEu8Z19cPyDk9PfiBFn0zTnHC2DaM56UbHsGtVfunr6spPCpwjB7IxCfyoHdAmj/h6AZCGbY03I/k+R823N2Vq1dSe0mC5hWBmEbxQeN2ZOFDTXq8wRL04csJPnpoUBbt0IKjdBdaX1rX1uNCZ0jOEiGpujbzAZyy1Cx45hEJBrHItrzZiiud90/DFzFPbRePlF6YHSlzP+TKKM7Q8MLvmpUXKOFfenynG4ziXrqk9WSV19GML6nHhCuX350vpCHNJdbG5CvTH5VdgiZWRbxzH6Yvs85JzHYFzF/abjj+HA99FYrxRPKO86C/yPUnZQOcJzWfB5qiBdw2Hq1u2gPHdKdrNDx7pK6ooZ0ig8MYJ7+hTRV3JuesbZWz0SW6+8oKEIt4G+2KwZXCO533T8UfPk7BcnprCtraFli4tRmxSccB75xwQel8MjTbgQLiRz4bEZ4SV1oW9VmIhnlfiF9sL6Ss7NwlhXVCRhWyGraS7MbQwkq7iSuDfIaxL+6Im668gRzY8ac7Pp3KNOBY5UfVbkPwcsNWx31678PNH1/Ix2TOCnpvgzXivFN/te+YXTBoHgkroCtb3ZKXh6hQ40lNRXcm4Ghry0KRXbUiGdypLcdlQPFq3hSuV+0/EPTk6nEW5aL0H4r2x4VP9Xxm3lH5fpCKzFykBloDJQGXAMyIbiHeAYZbvpbA0tfdTAVhofCta4hspAZaAyUBlIYEC2kyNdPI23u1ku9y/DGlFqxG3wTSnHuGqoDFQGKgOVgTQG2KTya7N74X9F01De26WPawAAAABJRU5ErkJggg==\n",
60
      "text/latex": [
61
       "$\\displaystyle \\left( 1, \\  y, \\  y^{2}, \\  x, \\  x y, \\  x y^{2}, \\  x^{2}, \\  x^{2} y, \\  x^{2} y^{2}\\right)$"
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
      ],
      "text/plain": [
       "⎛       2             2   2   2     2  2⎞\n",
       "⎝1, y, y , x, x⋅y, x⋅y , x , x ⋅y, x ⋅y ⎠"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "moments = exponents_to_polynomial_representations(moment_exponents)\n",
    "moments"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
85
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXMAAAF5CAYAAABtIcr0AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAAsTAAALEwEAmpwYAAAg5ElEQVR4nO3deXRlVYGo8W8nNc9AiqEoQApKVJRBWkGQQRww4mucQG2ccVZEQV83jd363lNpVKAFBBHHtsHXOLc2QW1BQUSeouKMIMhoUaSoEWpKst8fN7cquZWb3PGM328tlubmnGTrOvk4d599zwkxRiRJ+daT9gAkSe0z5pJUANNqXwgDgwF4DnAssBAIdfaNwDrgBuD7sb9vpFuDlBoRBganA/3AkcBcPHaVE2FgcBZwInA4MJvJj901wPWxv++6sd/YIebAPwGnNjGOlwFXj+4npeli4FlNbP8y4CvA+7szHGlqYWCwF7gCeHoTu50SBgY/F/v7zqu+MG6aJQwMLgZe2cJ4Tg4Dg3u0sJ/UEWFg8CCaC3nVyWFgcEmnxyM14Rk0F/Kq14aBwUXVL2rnzJ8ywWuNCMDBLewndcohbezrsas0HdTifr3Ak6tf1IZ7ZsvDaW9fqV0z2tjXY1dpauf423bcTzRnvqMzn78Xd/5mFr29la932nWIz/387jYGIHXfVy5exPVfWcj9d87gGS9Yz9mfWZH2kKSmfO+q+fzHhbuwasV0Fu4yxBn/uoKnHrdxok0biznAaR9YyUlvXtuxQUrdtsseQ5xyxipuvX4uWzbVWx0gZdNPr53Dl85dzPs+9SBPPmITDz84aa8bj7mUN8e/bAMAf/rVLFb91WNd+XLVx/o4+V2rOOioTQDsttfQZJs3frHzyo/1ccr++/Gu5+zNz38wu71RSpLqGh6Cu38/i7Wrenndofty6oHLuPBdu7LpsbrvMBuL+ev/+WE+f+td/Ptv7+KEU9fw4Tcs5b47pnds4JKk7Vat6GV4CG6+Zj4f+869XHL9Pdz9+1l88SO71NulsZg/5chNzF0QmTErcuLr1/H4Qzby04G5HRu4JGm7mbMrd0A88fWrWbznMDvtOsyL3vIIv/xh3e62eG+WEPFui5LUHQt3GWGnXYcIY2dVJr+GP3XM1z3Sw83XzGHzxsDQVrj2S/O5/dY5PP15j7Y5XKm7hrbC5o2BkWEYGWbbMSzlwbNeupbvfH4Rq1b0snZVD//56Z047PgN9Taf+gr/0NbAl87r46Nvm0lPT2SPfbfwD1c8wOOe6F+Fsu2LH96Fr31y+xzjTd9ZwEvfsYrTPrgqxVFJjXnt+1exbnUvb37GvkyfETmifz2vOfuReptPHfOddxvm0h/d29FBSkk47YOGW/k1fQacefFKzrx4ZSObez9zSSqA2pi3c1XTK6LKK49dpakjx19tzCf8zH+DHmtnIFKb2rkg77GrNHWku7Ux/wXQyoXNYeDnbQxIatctLe43DPyskwORmtTqsbsRuK36xbiYx/6+9cDHW/ihF8T+vjUtDkhqW+zvuwv4Ygu7XuixqzTF/r5fAV9rdjfgvNjft+2sPsQJPvwTBgaX0cQzQGN/351NDkTqijAweCBNPAPUY1dZEQYGDwaOYOpngK4Gfhj7++4Zt/9EMZck5YtLE1V4IYSvhRAuTnscUjd5Zq5CCyHsAdxN5ULnrjFGb0OhQvLMXEX3RirzjCPAKSmPReoaz8xVWCGEHmAFsHj0pd/HGA9McUhS1yR2Zh5C2DmE8N0QwmtDCLOS+r0qtecCY4+1x4UQDk5rMCqPEMK8EMI7QgjfTqp3SU6zbASOAi4HVoYQzgsh7Jng71f5vAeYP+brmcDpKY1FJRBC2C+EcCnwEHAh8FRgSxK/O7GYxxg3Ap+isn5yPnAGcOfov7mOCiH49HR1zOiFz+NqXu4FXhlC8ClZ6phQ8bwQwvXAb6lcp5lD5dP058UYR5IYR9IXQC+iciEKKmdJs4ATge8CtzsFow6qXvis5YVQdUR1KgW4l8onOI+j0rTq85ED8IXExpP0BdAQwgBwAhN/wmkDlT/Ay4CLYowPJDk2FcMEFz5reSFULQsh7AecBbyWSq8meqc3BHwuxviWpMaVxtLEc6l/h7t5jJ+C+Y5TMGpB7YXPWl4IVVMmmUqpN2W3FTg/qfFBOmfmAbgD2K+BzSOVWzw+ADwzxvhwN8emYgghXEvl3V89w8AXYoxvTGhIyrEQwv7A9cAiKiecjfhxjPHorg1qAomfmcfKvz0+QmVKZSqByk1ndsYPOKkBdS581vJCqJoRqNx0sNGQb6DSuESlFcgvN7jdCLAKODzG+FAXx6PiqHfhs5YXQtWQGOMdwNFU7rTZiHVUFnUkKpWYjy5TvJzJ119WQ35EjPGuRAamXBu98Hk6k8+XV80D3tvdEakoYoy3AccwddAfI8HliGOl9nH+EMLewO1M/Ic3AqwBnmbI1agQwvOA/wKmNbjLCHBojPHX3RuVimT0wvlPqX/CsBHYPcbY6Fl8xzR60HdcjPHe0SvDz2f8MsURKu8YdgY2pzE25dYfgKsneP3vRv/zqprXh6msEZYaNUT9kA8BX0oj5JDyjbZCCMdQOZOqXlioTq0cQ+UPE2Cp683VjhBCBB6KMe6e9liUXyGEA6ksS4TKx/R/CCwYs8km4OAY458SHhqQ/gqRG6ncwwDGz5H/ke3/9rvfe7hISlNNyHtjjL9kxzn0n6cVckg55mOWKW6l5mJnjHEzBl1SyiYI+QjscFF0MyksRxwr7TNzqCxT/DwTrFox6JLSVC/kVWOC/mlSWI44Vi4eThFCmEllPgqcQ1eTnDNXK6YKedZk4cx8Sp6hS0pS3kIOOYk5GHRJychjyCFHMQeDLqm78hpyyFnMwaBL6o48hxxyGHMw6JI6K+8hh5zGHAy6pM4oQsghxzEHgy6pPUUJOeQ85mDQJbWmSCGHAsQcDLqk5hQt5FCQmINBl9SYIoYcChRzMOiSJlfUkEPBYg4GXdLEihxyKGDMwaBLGq/oIYeCxhwMuqSKMoQcChxzMOhS2ZUl5FDwmINBl8qqTCGHEsQcDLpUNmULOZQk5mDQpbIoY8ihRDEHgy4VXVlDDiWLORh0qajKHHIoYczBoEtFU/aQQ0ljDgZdKgpDXlHamINBl/LOkG9X6piDQZfyypCPV/qYg0GX8saQ78iYjzLoUj4Y8okZ8zEMupRthrw+Y17DoEvZZMgnZ8wnYNClbDHkUzPmdRh0KRsMeWOM+SQMupQuQ944Yz4Fgy6lw5A3x5g3wKBLyTLkzTPmDTLoUjIMeWuMeRMMutRdhrx1xrxJBl3qDkPeHmPeAoMudZYhb58xb5FBlzrDkHeGMW+DQZfaY8g7x5i3yaBLrTHknWXMO8CgS80x5J1nzDvEoEuNMeTdYcw7yKBLkzPk3WPMO8ygSxMz5N1lzLvAoEvjGfLuM+ZdYtClCkOeDGPeRQZdZWfIk2PMu8ygq6wMebKMeQIMusrGkCfPmCfEoKssDHk6jHmCDLqKzpCnx5gnzKCrqAx5uox5Cgy6isaQp8+Yp8SgqygMeTYY8xQZdOWdIc8OY54yg668MuTZYswzwKArbwx59hjzjDDoygtDnk3GPEMMurLOkGeXMc8Yg66sMuTZZswzyKArawx59hnzjDLoygpDng/GPMMMutJmyPPDmGecQVdaDHm+GPMcMOhKmiHPH2OeEwZdSTHk+WTMc8Sgq9sMeX4Z85wx6OoWQ55vxjyHDLo6zZDnnzHPKYOuTjHkxWDMc8ygq12GvDiMec4ZdLXKkBeLMS8Ag65mGfLiMeYFYdDVKENeTMa8QAy6pmLIi8uYF4xBVz2GvNiMeQEZdNUy5MVnzAvKoKvKkJeDMS8wgy5DXh7GvOAMenkZ8nIx5iVg0MvHkJePMS8Jg14ehrycjHmJGPTiM+TlZcxLxqAXlyEvN2NeQga9eAy5jHlJGfTiMOQCY15qBj3/DLmqjHnJGfT8MuQay5jLoOeQIVctYy7AoOeJIddEjLm2MejZZ8hVjzHXOAY9uwy5JmPMtQODnj2GXFMx5pqQQc8OQ65GGHPVZdDTZ8jVKGOuSRn09BhyNcOYa0oGPXmGXM0y5mqIQU+OIVcrjLkaZtC7z5CrVcZcTTHo3WPI1Q5jrqYZ9M4z5GqXMVdLDHrnGHJ1gjFXywx6+wy5OsWYqy0GvXWGXJ1kzNU2g948Q65OM+bqCIPeOEOubjDm6hiDPjVDrm4x5uoog16fIVc3GXN1nEHfkSFXtxlzdYVB386QKwnGXF1j0A25kmPM1VVlDrohV5KMubqujEE35EqaMVciyhR0Q640GHMlpgxBN+RKizFXooocdEOuNBlzJa6IQTfkSpsxVyqKFHRDriww5kpNEYJuyJUVxlypynPQDbmyxJgrdXkMuiFX1hhzZUKegm7IlUXGXJmRh6AbcmWVMVemZDnohlxZZsyVOVkMuiFX1hlzZVKWgm7IlQfGXJmVhaAbcuWFMVempRl0Q648MebKvDSCbsiVN8ZcuZBk0A258siYKzeSCLohV14Zc+VKN4NuyJVn02pfCAOD04GXAM8EFgChzr4RWA/cBHw99vdt7tYgpbFijJtDCLOATVSCvjTG+EDdY/esS2DajIVhYPDfqj8CWA1cD3w79veNGHKlKQwMzqZy7B4JzKN+d0eAVcB/x/6+gbHf2CHmwEXA8U2M47mj27+piX2ktkwUdK55+Bzg2TtsfMBh0NM7HTi85jv9wFNDCFdjyJWSMDAYgE8BRzSx2wvDwOCBsb/v49UXxk2zhIHB5TQX8qpjwsDgk1rYT2rZuCmXJcvuZ2jrCU3/kC2bXsW8RYZcaTqU5kJe9ZowMDin+kXtnPmT2xjQU9rYd5wQQr23GNI424K+zxNgxT3LGNo60bvNiW3ZPIOV9y9j7wPAkKsJHW5Uq92dCSyvflEb8xktD6e9fQkVzwwhfBVYFULYt52fp/KIMW7mrE++Cmg86Fs2z2DlfY8D4NyvP8eQq1EhhMOA1SGEL4QQntqBH9mR7ja2muVDr9udVz5hP16yz/68/rB9+danF7bxy8cJIewcQng38BdggMpFgOnAzp36HSqBOfOG2HPZHcD4oK8d7OHSv4e3H9PLq5+yjO/++/xxId9zvz8xfUZMa9jKpd2BXuBVwI0hhD+GEN4UQpjfld92zx+n87d7LudDr9t9ss0ae0v6ijMfYe/HP8SMWZG7fzeDs1+yF8sP2cSTnt7SCpbRtyhHAe8GTqRyhXbOmE2GW/m5KrnQE9lz2R08cNdyVtyzjN33uYtPnLmYadPhgmuHWTv4Vz702qXstGtgybJKyJ3RU2uGqQR9DnAAcAHwidGL6RfFGH/Rsd/0yfftxrIDN021WWNn5vsftIUZsypnLyFEQoAH7mr6rUGds/BZjA+51Lpq0AHu+eMyfvb9+Zz0Zpg1Bw48YpiDjg789FpDrk6bB8xmx7P1eW391O9dNZ85C4Z5ypGPTbVp4x8auuD0XXnR0uW8/dh9WbR4iKNO3NDIbjVz4Q8AHwb2ZvK1lFLrqkF/6F7o6YXd9gYirLzvcey1P6y451FDri6pPVtfOTq3fmjTP2nDmh6+fH4fbz334UY2b/zK/5kXr+SMC1fy65tmc9uNs5k+c8p5xtGz8PdQmf+eS+Px7gGOCSEsbnh8Krf3XXYwRzx//Du83mkPM3tu5RgaGekFYGHfGjZumM3GDdu3Hfi3p4UXfKCtC/gqlSNorGXVs/JXASeHEO4Dzo8xXtHQb/ns/+rj+FPWsvveQ41s3njMAXqnwaHHbuQHVy/gG5ct4pQz1kyxx4VU5sObvW3AfCr/VpMac8t3Yfkh418bHoKNj45/be2qRcycDatWLN322s++/y/dH6BKrHq2vhy4HJg65rffOpPf/GQOl97wl0Z/SXMxrxoZgr/+pZEzmScA7wReN/p1o/NHa4FnxxhvbWF0KqEwMPhy4H+Pe3H+TjMZGd6Hh+6D3faqvPbXuzewzxO3sHT/wW3bnfuNt8X+vusSHK5yLIRwInAl0OiqvvVULpheRiXmU/vlDXMYfHA6rzloPwA2b+xhZATedvRMLrvxnol2mfqMedWKXr531XweWx8YHoKbr5nDT65ZwCHHPDrVrjHG22OMpwOLgTcDvwA2Alsb+h8ktWrL5hmsX70Phx4H/3kFbN0yzJ23wa3Xz+PZp0x57Ept2kLlVhM3AK8GFscY/zHGOGGId3DSm9fwmVvu4pLr/8Il1/+F57x8DYcc/Sgf/ur99XaZ+sw8BBj44iIuP2c34gjssscQr3v/So59ccN/EDHGTcCXgS+HEA6gtbN1qTFj15GfedGdnPfW/XnPCb3MWzjEqe+bxpz5ezG09S6mTW9oLlJqwriz8IbjXWv23MjsuduXaM+aO8L0mSPsvFvdZdtTx3zn3Ya58Lv3tTSgCcQYbwdODyG8D3gx8F7giaNjmd6p36OSqv1AUAjwjvOgp3eYJfveRRwJ49ahG3S1bwuVa4M/A84H/ivG2Nnj6rQPrppqk9TuZx5j3BRj/HKM8TAqN5q5HNgw+s+sSXeWJjJRyGuNXYfe7L1cpO1mUTkLX0NlscYTYozHxBi/1fGQN6g25u18rLnlfSeYW/8WlQ8WSY359U17Thnyqtqg/+6WvgRGqOL4DfBtts+Fn93ydEpFR7pbG/N2Lgw19CGiyYw5W395jHHKtxUSjD4h6KuXnA80/snOsUG/4p/+I4mHRKsYYoz3xhhP7uBZeDvd3bZvbcxvoTL306wI3NzGgKSWbHtC0O23wp77/aGpT3aGnsiSZb/irt9Blx8SLU3iJy3uNwjcXv1iXMxjf98g8NEWfugFsb/voRYHJLVk3KPe1q/uJYSPNfkjhunp+QBDW7r6kGhpMrG/717gkiZ32wr8c+zv23byHWLccbomDAwuYftzFCezDrgp9vc90ORApLbUe2bn6LF7FGM/0PGJ91zBtOnreMdHzxrzIx4Bbhg9gSGEMJPKumCApTFGj2klKgwM7s32Z4DWE6mckd8Q+/tWj9t/ophLWdbsw5dDCBF4KMY46f2gDbryLLWliVIrmg15M8Y9U9QpF+WMMVdudDPkVQZdeWXMlQtJhLzKoCuPjLkyL8mQVxl05Y0xV6alEfIqg648MebKrDRDXmXQlRfGXJmUhZBXGXTlgTFX5mQp5FUGXVlnzJUpWQx5lUFXlhlzZUaWQ15l0JVVxlyZkIeQVxl0ZZExV+ryFPIqg66sMeZKVR5DXmXQlSXGXKnJc8irDLqywpgrFUUIeZVBVxYYcyWuSCGvMuhKmzFXoooY8iqDrjQZcyWmyCGvMuhKizFXIsoQ8iqDrjQYc3VdmUJeZdCVNGOuripjyKsMupJkzNU1ZQ55lUFXUoy5usKQb2fQlQRjro4z5Dsy6Oo2Y66OMuT1GXR1kzFXxxjyqRl0dYsxV0cY8sYZdHWDMVfbDHnzDLo6zZirLYa8dQZdnWTM1TJD3j6Drk4x5mqJIe8cg65OMOZqmiHvPIOudhlzNcWQd49BVzuMuRpmyLvPoKtVxlwNMeTJMehqhTHXlAx58gy6mmXMNSlDnh6DrmYYc9VlyNNn0NUoY64JGfLsMOhqhDHXDgx59hh0TcWYaxxDnl0GXZMx5trGkGefQVc9xlyAIc8Tg66JGHMZ8hwy6KplzEvOkOeXQddYxrzEDHn+GXRVGfOSMuTFYdAFxryUDHnxGHQZ85Ix5MVl0MvNmJeIIS8+g15exrwkDHl5GPRyMuYlYMjLx6CXjzEvOENeXga9XIx5gRlyGfTyMOYFZchVZdDLwZgXkCFXLYNefMa8YAy56jHoxWbMC8SQayoGvbiMeUEYcjXKoBeTMS8AQ65mGfTiMeY5Z8jVKoNeLMY8xwy52mXQi8OY55QhV6cY9GIw5jlkyNVpBj3/jHnOGHJ1i0HPN2OeI4Zc3WbQ88uY54QhV1IMej4Z8xww5EqaQc8fY55xhlxpMej5YswzzJArbQY9P4x5RhlyZYVBzwdjnkGGXFlj0LPPmGeMIVdWGfRsM+YZYsiVdQY9u4x5Rhhy5YVBzyZjngGGXHlj0LPHmKfMkCuvDHq2GPMUGXLlnUHPDmOeEkOuojDo2WDMU2DIVTQGPX3GPGGGXEVl0NNlzBNkyFV0Bj09xjwhhlxlYdDTYcwTYMhVNgY9eca8ywy5ysqgJ8uYd5EhV9kZ9OQY8y4x5FKFQU+GMe8CQy6NZ9C7z5h3mCGXJmbQu8uYd5AhlyZn0LvHmHeIIZcaY9C7w5h3gCGXmmPQO8+Yt8mQS60x6J1lzNtgyKX2GPTOMeYtMuRSZxj0zjDmLTDkUmcZ9PYZ8yYZcqk7DHp7jHkTDLnUXQa9dca8QYZcSoZBb40xb4Ahl5Jl0JtnzKdgyKV0GPTmGPNJGHIpXQa9cca8DkMuZYNBb4wxn4Ahl7LFoE/NmNcw5FI2GfTJGfMxDLmUbQa9PmM+ypBL+WDQJ2bMMeRS3hj0HZU+5oZcyieDPl6pY27IpXwz6NuVNuaGXCoGg15RypgbcqlYDHoJY27IpWIqe9BLFXNDLhVbmYNempgbcqkcyhr0UsTckEvlUsagFz7mhlwqp7IFvdAxN+RSuZUp6IWNuSGXBOUJeiFjbsgljVWGoBcu5oZc0kSKHvRCxdyQS5pMkYNemJgbckmNKGrQCxFzQy6pGUUMeu5jbsgltaJoQc91zA25pHYUKei5jbkhl9QJRQl6LmNuyCV1UhGCnruYG3JJ3ZD3oOcq5oZcUjflOei5ibkhl5SEvAY9FzE35JKSlMegpx7zEMIeIYT7QggvrvN9Qy4pcY0GPYTwlhDCn0MI85Ib3Y5SjznwTmA34MoQwovGfsOQS0rTVEEPIbwFuABYArw64eGNE2KM6f3yEGYCK4EFoy9tBP4uxvhNQ65OCSFE4KEY4+5pj0X5NNqqTaNfLo0xPjAm5HNGX78XeFxMKappn5mfXDOG2cBVIYT/iSGXlBETnKGfw/iQA+wMHJ/02KpSOzMPIQTgD8ABk2xmyNWwEMKRwLeBUPOtnUb/c3XN60PA0THG27s9NhVDzRn6RK6LMT47qfGMleaZ+eHA0km+vxH424TGomJYC8yjEu+x/1TVvr4IWJPoCJV3rwMem+T7R4YQ9k1oLOOkGfN/YPxblFrVKZcXJTMc5V2M8XfAH5vY5b9jjA91azwqlgnmyCfSC7wnmRGNl0rMQwhLgBPY8e1wLYOuZn0c2NDAdhuo/GFKU2ow5ADTgTeEEOZ2f1TjpXVm/g6mDnnVbOBrIYT9ujgeFcdXaezYehS4rstjUQGEEI4ALmPqkI+V+DLFxGM+egHhncDMBjbfAKyncrb1YDfHpWKIMW4EvkTl4mY9G4GLvLiuBt0OfIrKXPmjDWw/Fzh7dJFHYtI4M69djlgrUvk/7A7gdGDXGOPfj/6RSo24BNg6yfd7gM8mNBblXIxxdYzx7cDuwN8DDzD1VF7iyxQTXZo4xXLEzVRC/gPgXOAnaS2+V/6FEG4DDqrz7YEY4wuSHI+KI4TQAzwH+Efg6cA0KnPltRJdppj0mflEyxGrUymfAJbHGF8YY7zJkKtN9S6EeuFTbYkxjsQYvxdjPI7KCcNnmHgKJtFlikmfmX+T7WvHH6MyD/4R4P/GGCdbiC81JYQwG3iYyvzlWA8BS5wvVyeFEOYDr6Wy5Hohlc87bAU+FWN8VxJjSOzMPISwO/A/gGHgGipLEw+IMX7BkKvT6lwI9cKnuiLGuD7GeAmwN/BS4EdU1py/MYTQzCqYliV2Zh5C6AXeBHwnxnh/Ir9UpTZ6s7afUVneCpXrMvv4QSElIYSwHDga+HwS08ap3jVR6raaC6Fe+FRhpX3XRKnbqhdCvfCpQvPMXIU2eiF0kMqKKS98qrCMuSQVgNMsklQA02pfCAODs4FTgWey/XFu9awHfgxcGfv7JrvHr9R1YWBwPpVj90gq63wnswa4Ebgq9vdt7vLQpEmFgcHFwKuAv2H76qt61jDBsbtDzKncHewZTYzjCODY0YFIqQgDg9X7rRzcxG5HUflU8lu7MiipAWFgcB5wJbBPE7sdRaW9b6m+MG6aJQwMPonmQl71tDAwWO8+GFIS/obmQl71rDAwmMqTYaRRz6W5kFcdFwYGl1W/qJ0zf0IbA3piG/tK7Wrn+PPYVZo6cuzWxnyiO381akYb+0rt8thVXrVz7G7bd6I58/FevPfycV9v2Rx47ivW8O5PrGxjAFL3bdkUuPBdu/Lbm+eyYV0vuy7dwmvOHuSoFzbygAEpXQ/cNY2Lz9qNO2+bzbTpkcNPWM/p569k2sTtnzrm37j3jm3//bH1gVMP3J9jXrS+YwOWumVoCPqWDPEv37qXPfYZ4qbvzOXj71jC4550N3sum+xJRFL6Lj5rNxbuMsyVv/sz61f3cPZL9uLrly7ilDPWTLR5c+vMr//qfObvPMShx/rUH2XfnHmR0z64ij2XDdHTC0ef9CiLl2zh9ltnpT00aUoP3z+dY05az8zZkb4lwxxyzKPce3vdx202F/PrvrKAY1+0juBnjZRDq/7ay4p7Z7Dvk7akPRRpSi98w2p+9I35bHw08NB90/jVDXM57Nl1pwgbr/KDd0/jj7fO4YRXr+vIQKUkbd0C575pD445aR37HmjMlX0HH72R++6Yycn7Lef1T13Gsidv4rgX1332aOMx/96VC3j8IRtZut9kD8qVsmdkGD5y2h5Mmx4541+9l7myb2QYPvDKpRzx/PV8/S93cNXv72TD2h4uO3txvV0aj/mPvrGQ409e25GBSkmJI/DRt+7O2sFpfPDKB5nuKkTlwNpVvaxaMY2XvH0NM2ZFFi0e4bmvWMcvf1T7GMRtGov5bTfOYvXKaRx/iqtYlC/nn74b9/95Bh+6+n5mzfEWocqHnXYdpm/PrXzz8kUMbYV1j/Twg6sXsPfj695HaOqliQDf//JCnvac9cxd4B+D8uPBu6dx3dULmTYjcuqB+297/S0fXsHzX+2JibLtnM8+yOXn7Mq3Pr0zPT2RJx3+GG8/r+7nexqL+XsvdZ5R+bNk3yGuefhPaQ9DaskBh23mgmvva3Rz1xhKUgHUxrydaRSnYJRXHrtKU0e6WxvzduYRXX+uNLVz/Dl/rjS1c/xt27c25jcDwy38wBHgJ20MSGrXj1vcbwtwSycHIjXpxhb32wL8v+oX42Ie+/vWAB+guaAPA/8n9vcNtjggqW2xv+9B4DyaO3aHgHNif593UVRqYn/fz4EvNrnbEPD+2N+37ROhIcYdp2vCwOAiKo8kWgCEemOgcor/09jft7rJgUhdEQYGd6Fy7M5l8mN3LXBz7O9zelCZEAYGl1B5YtYsWjh2J4y5JClfXJooSQVgzCWpAP4/Vte8hUm3Po0AAAAASUVORK5CYII=\n",
86
87
88
89
      "text/plain": [
       "<Figure size 1152x432 with 1 Axes>"
      ]
     },
90
91
92
     "metadata": {
      "needs_background": "light"
     },
93
94
95
96
97
98
     "output_type": "display_data"
    }
   ],
   "source": [
    "from lbmpy.stencils import get_stencil\n",
    "d2q9 = get_stencil(\"D2Q9\", ordering='walberla')\n",
99
    "ps.stencil.plot(d2q9)"
100
101
102
103
104
105
106
107
108
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
109
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWQAAADhCAYAAAD/Ec//AAAACXBIWXMAAA7EAAAOxAGVKw4bAAAb1UlEQVR4Ae2dUW7cRraGpQs/G8YEyAKkHVjjFcTage0d+GYHee15C3x3YHsFQbwDJyswrB3MXUCAawRZwf3/nq6ZNs3uU93N6jolfQVQRVaRPP/5qnVEFcnTl6vV6unFxcVnLXPlwz/+8Y+Xcx20QQACEIDAYQQUT/+pI67mjlLf5aOtjv/RunfeLv+7vcE6BCAAAQicRODNzNG3anvh9u2A/FYRmgA8Q4smCEAAAksQUIx9Nz2P2tz0TUCe7je7rYM9xfGrlhut/zm7U6PGnrYjlzJrG117ZrZoiz5d8/2Zuc0r/rq1lf7tK+SvLW5tyfgTbb7X8kXL37XMzoGoffHS03bkTGZto2vPzBZt0adrvj8zt3nFX7eeQ39tQPaV8PrmnkT9pHVfJZ+lyF4325GDmbWNrj0zW7RFn675/szc5hV/3XoO/f/1tUm2IAABCECgFwECci/y2IUABCAwIUBAngBhEwIQgEAvAgTkXuSxCwEIQGBCgIA8AcImBCAAgV4ECMi9yGMXAhCAwIQAAXkChE0IQAACvQgQkHuRxy4EIACBCQEC8gQImxCAAAR6ETgmIH+3Efu3DqJ72o7czaxtdO2Z2aIt+nTN92fmNq/469Ym+i+38iFf69XAndne1OeEQi7PtTi3xZ0W7/9Rfd9kMFL7YqWn7ciJzNpG156ZLdqiT9d8f2Zu84q/bm2hX+f8b1lxts3L6oD8tSy2IAABCEBgCQLbAfmYKYslNHAOCEAAAhCYECAgT4CwCQEIQKAXAQJyL/LYhQAEIDAhQECeAGETAhCAQC8CBORe5LELAQhAYEKAgDwBwiYEIACBXgQIyL3IYxcCEIDAhAABeQKETQhAAAK9CBCQe5HHLgQgAIEJAQLyBAibEIAABHoReHSIYb3i92az//+pvtbyRm07818ccu6afWXrqfZzTo0brf9Zc0yWfbJql66uYxqNzwD60n4ms37mojF3/0PVXh2QBeizOP2s+sMGmBMMfdb2rZZmQVnntp33Wr5o+buWKy1DlOzapa/LmNYOXlZ9mcc1s7Zo3NF+cVEVkAXK2YieqF4HY4PV+p+b7bfavHVbi2I7Ou9Ln1vrP6nyFckQJbN2aes2pjWDl1lf8nHl96XmA7bwPkt9JmrnkB0QnW5zWj6p4bnE+CqWMhaB7GOaXd9Yo43aIQjUBmTnQPaUwbSUqQr3U8YikH1Ms+sba7RROwSBMCBXXv32+PaQIQBnFJl9TLPryzimaLofBMKALDdLsPXc1K7ClMUuMjnbs49pdn05RxVVwxOoCcg1Tpbvl6rZl33GIJB9TLPrG2OUUZmKwCOpebxRVOqpwLm547JPuZLxc8n3rmz+df5djh3yH8BLHXeXHEb2Mc2uL/nw9pE38u9LZ+3flxFzQN5bJNSPt3mfuaBU2srNvb3nGq3TvkvzzWi6I73ZxzS7vojvQ+0f+fcli3ZPWfy1+QCVeu7z9Jsar2Y6yhWy+yljEcg+ptn1jTXaqM1M4I8irnYO2a8r+y25afHV493mr8u0j+3cBLKPaXZ9uUcXdUMSqArICrjv5N0X1S+Kl1r3dMUrLa9L2xnqciOnXJmfweRiJlJpTzSms4Cz69sSnWpct3R5NbO2idRvNh+k9svVavVUKJzT4Fq/BDvngjcB2IloPK/qm3jPtDi3RfMbWLLhqyUXvyzgPwS2aa0f1ec/FmlLZu3SZpZdxrRmwDLrSz6u/L7UfMAW3ufYz4SOcxqDt6ovqwPywto5HQQgAAEIiMB2QK6asoAaBCAAAQi0J0BAbs8YCxCAAASqCBCQqzCxEwQgAIH2BAjI7RljAQIQgEAVAQJyFSZ2ggAEINCeAAG5PWMsQAACEKgi4IDs53l/1LIvoUvVydgJAhCAAAQOJuA0AY7BFw7IV1r8vXgjvv0m2RQIQAACQxPwC2+OweuAPLQniIcABCBwXwgwh3xfRhI/IACB4QmE+ZC3PdQrfs574OJcFtda3qhtZ/4L77hkkS3n3fB7+jdad06NNCWztn2QpLvrmO7T5r4B9PGZjAZxpn/U35fiSiv91QFZApyAyMmEPliUaiem+az6VkuzoLyx8162fNPRKUA9552iZNZWA0j6u4xpjTbvk1Vf5nFHW+2n6/D9zsG2KiBLiLMRPVG9DsZ2Rev+JhFvezL61m0tiu3ovC99bq3/pMpXJClKZm0RIGnvNqaRNvdn1pd53NFW8+k6bp9zsK2dQ3ZAnEuz+UntzyXUV8uUsQhkH9Ps+sYabdQOQaA2IPuxjLnnlMtUhfspYxHIPqbZ9Y012qgdgkAYkCuvfnmGeYjh/pfI7GOaXd9AQ43UwQiEAVn+lGC776kGpizGGvjsY5pd31ijjdphCNQE5Bpnyvdf1ezLPmMQyD6m2fWNMcqoTEXAAfnxRlGppwLn5o7LPuVKxs8lU8YhkH1Ms+sbZ6RROgKB74tIB+S/NhulLn3rWvN5ZapiblqitJWbe18dy0ZOAtnHNLu+nKOKqoEJ/FG0105ZOBvRVTloqy5XyO6njEUg+5hm1zfWaKN2CAK1AdmvK/stuWm5UcPd1hXNtJ/tvASyj2l2fXlHFmXDEqgKyAq47+ThF9Uviqda93TFKy2vS9sZ6nIjp1yZn8FktYnM2r5xItGYfqPNDdn1bYnOPO5o2xqohVebsL1crVZPJdQ5Da71S7BzLlh9DsBORPOnFt/Ee6bFuS3uVDctsuGrJRe/LGAdtmmtH9XnPxbdSmZtERRp7zamkTb3Z9aXedzRVvPpOm6fFmx1TqcxeKv6sjogHyefoyAAAQhAYB+B7YBcNWWx72T0QQACEIDAMgQIyMtw5CwQgAAETiZAQD4ZISeAAAQgsAwBAvIyHDkLBCAAgZMJEJBPRsgJIAABCCxDgIC8DEfOAgEIQOBkAg7Ifp73Ry37ErqcbIgTQAACEIDALAGnCXAMvnBAvtLi78XL+PabZFEgAAEI3GsCfuHNMXgdkO+1pzgHAQhAYBQCzCGPMlLohAAE7j2BR4d4qFf8nMvCxbksrrW8UdvO/BfecanS03bkQ2ZtkfaafvnnfCfOJ3KjdecyOWvpbf8UZ3tq72n7FGY+Nrv2VvqqA7IEOAGRkwl92ABzYprP2r7V0jQo97RtX/eVzNr26Y765JfH970W3+x16lXfazhb6W3/FEd7au9p+xRmPja79nPoqwrIEuJsRE9Ur4PxBt6fm21PRt+6rUXpaTvyJ7O2SHvUL998JfzS+2n9J1W+Sj5b6W3/FEd7au9p+xRmPja79nPoq51D9i/mXJrNT2p/LqG+mmpVetqOfMqsLdJOPwQgkIxAbUD2YxlzzymXqQr3tyo9bUc+ZdYWaacfAhBIRiAMyJVXv02eYe5pOxqnzNoi7fRDAAI5CYQBWbJLsPWc4q7Sasqip+1dvpb2zNqKRmoIQGAgAjUBucad8v1SNfsuvU9P25EvmbVF2umHAATOTOCR7D3e2Cz1VMLc3HHZp1wl+rnkFqWn7cifzNouNlMqv8uJQ/57eanj7iLH73s/7I4bYbgdx01HfV+OdED+a7NR6tK3rgXZj7d5fe4Xu7SVm3vrY5b60dN25ENmbdZufapuIj/o/5YA7L5lUtMCtxpKs/v8UVprpyycjeiqHLRVlytk97cqPW1HPmXWFmmnHwIQSEagNiD7tVm/rTUtvgK72/xlnPYttd3TduRDZm2RdvohAIFkBKoCsgLuO+n+ovpF0a91T1e80vK6tLWoe9qO/MmsLdJ+YH+5OVn+Izrw8JN3723/FAd6au9p+xRmPja79ib6Ller1VM57zwV1wowO+eCNwHYyYU8N+mbeM+0OLdF85tAPW3Lx70ls7a9wis65Zv/A3DxCzD+A+yx9mfko/r8R7pp6W3/FOd6au9p+xRmPja79hb6dE6npnir+rI6IJ8KmuMhAAEIQOBbAtsBuWrK4ttT0AIBCEAAAksTICAvTZTzQQACEDiSAAH5SHAcBgEIQGBpAgTkpYlyPghAAAJHEiAgHwmOwyAAAQgsTYCAvDRRzgcBCEDgSAIOyH6u9Ect+5LlHHl6DoMABCAAgYCAUzA4Bl84IF9p8ffi9XoLS6YpEIAABB4sAb945Ri8DsgPlgKOQwACEMhEgDnkTKOBFghA4EETcD7k6qJX/JzLwsW5LK61vFHbzvwX3nGp0tN25ENmbaNrH5mt2Uu/c8U4J8iN1p0H5mylp+3IyczaIu1R/ym+VQdkGXECIicT+mBBqp1s5rPqWy1Ng7LO3822fd1XMmvbp9t92bVn17eLr3T7d+O9Ft8od9pa36c5S+lpO3Iws7ZIe9S/lG9VUxYy5mxET1Svg7HFad1/7b29nox2W4vS03bkT2Zto2sfnK2/Zcdfh+U7579EY7Fkv2x2sx35kVlbpD3qX8q3qoAsMS+1zKXZ/KT25xLjK4JWpaftyKfM2kbXPjLbiD39EJglUBuQ/VjG3HPKZarC/a1KT9uRT5m1ja59ZLYRe/ohMEsgDMiVV79NnmHuaXuW1lZjZm1bMmdXs2vPrm8WKo0QWIBAGJBlowRbzxnvKq2mLHra3uVrac+srWjcVWfXnl3fLq60Q+AkAjUBucZA+X6pmn2X3qen7ciXzNpG1z4y24g9/Q+UwCP5/Xjje6mnKObmjss+5UrGzyW3KD1tR/5k1ja69q5sN1MmvwviIf/5+amKuwg8/TkJdB7z7wsVB+S/NhulLn3rWkL9GI3X5z6cpa3c3Fsfs9SPnrYjHzJrG117b7a2L4Y3EUf67w+BzmP+RyFZO2XhbERX5aCtulwhu79V6Wk78imzttG1j8w2Yk8/BGYJ1AZkv/rpN46mxVcRd5u/LtO+pbZ72o58yKxtdO0js43Y0w+BWQJVAVkB952O/qL6RTmL1j1d8UrL69LWou5pO/Ins7bRtY/MdsK+3Hws/01Ouptu9rQdOZZZW6Q96j/aN88h15Yb7ehkQs9U+yae6x+0fY4bGT1ty829JbO2vcLVmV17dn07+er3wlf4Ls//VV38qjbfa/mo2hc4zUpP25FTmbVF2qP+JXy7XK1WT2XIyXuudcImN+ciR+iHAAQg8FAJKO46V9Bb1ZdVUxYPFRR+QwACEDgnAQLyOWljCwIQgMAeAgTkPXDoggAEIHBOAgTkc9LGFgQgAIE9BAjIe+DQBQEIQOCcBAjI56SNLQhAAAJ7CDgg+1E3f9XMvoQue05BFwQgAAEInEDAaQIcgy8ckK+0+HvxerxJJLMUCEAAAg+agF8eWn83KVMWD/pzgPMQgEAmAgTkTKOBFghA4EETOCSXxYVe7XuzoeVcFtdanNviLK9b97QdfUIyaxtd+8hszV76nZrAeS1utO48y2crPW2f6mR27a30VQdkCXC+i59VfzBs1c729ln1rZamQVnn72bbvu4rmbXt0+2+7Nqz69vFV7r9u/Fei2+UO22t79OcpfS0faqD2bWfQ1/VlIWEOPnFE9XrYGzwWvdfe2+vJ6Pd1qL0tB35k1nb6NoHZ+tv2fFXOvnO+S/RWCzZL5vdbJ/qR3bt59BXFZAF+qWWuTSbn9T+XEJ9RdCq9LQd+ZRZ2+jaR2YbsacfArMEagOyH8uYe065TFW4v1XpaTvyKbO20bWPzDZiTz8EZgmEAbny6rfJM8w9bc/S2mrMrG1L5uxqdu3Z9c1CpRECCxAIA7JslGDrOeNdpdWURU/bu3wt7Zm1FY276uzas+vbxZV2CJxEoCYg1xj4rmanRvv0tB25lFnb6NpHZhuxp/+BEngkvx9vfC/1FMXc3HHZp1zJ+LnkFqWn7cifzNpG196V7WbK5HdBPOQ/Pz9VcReBv8/9cDt6dL8vRzog/7XZKHXpW9eC7MdovD734Sxt5ebe+pilfvS0HfmQWdvo2nuztX0xvIk40v81Abh9zeOArT/KvrVTFs5GdFUO2qrLFbL7W5WetiOfMmsbXfvIbCP29ENglkBtQParn37jaFp8FXG3+cs47Vtqu6ftyIfM2kbXPjLbiD39EJglUBWQFXDf6egvql+Us2jd0xWvtLwubS3qnrYjfzJrG137yGwn7MvNx/Lf5KS76WZP26c6ll17E32eQ64tN9rRyYSeqfZNPNc/aPscNzJ62pabe0tmbXuFqzO79uz6dvLV74Wv8F2e/6u6+FVtvtfyUbUvcJqVnrZPdSq79tb6Ller1VNBdPKeaxlrcnPu1EHieAhAAAL3lYDirnMFvVV9WTVlcV9B4BcEIACBTAQIyJlGAy0QgMCDJkBAftDDj/MQgEAmAgTkTKOBFghA4EETICA/6OHHeQhAIBMBB+THG0GlzqQPLRCAAATuO4F/57JwQC45LEp9353HPwhAAAKZCBycyyKTeLRAAAIQuJcEmEO+l8OKUxCAwIgECMgjjhqaIQCBe0ngkFwWF3q1782GgnNZXGtxbouzvG7d03Y08pm1RdrdL/1+fd65F2607lzAqUpmfWg77qOSmVuNR630VwdkCXC+i59Vf7Bg1c729ln1rZamQVnn72bbvu4rmbUFuj1+77V80eLUqlda0hRxTasPbcd9TDJzq/HoHPqrArKEOPnFE9XrYGzxWvc3iXj7rZZbt7UostHNduRPZm0V2n0l/NL7yY+fVPkqOU2RprT60HbcxyQztxqPzqG/dg7Zv7hzaTY/qf25hPpqplXpaTvyKbO2SDv9EIBAMgK1Adk5Xf2v7bSUqQr3tyo9bUc+ZdYWaacfAhBIRiAMyJVXv39r4VdP25E/mbVF2umHAARyEggDsmSXYOs5vV2l1ZRFT9u7fC3tmbUVjdQQgMBABGoCco0739Xs1GifnrYjlzJri7TTDwEInJmAA7LngX/UMjdHbDm72t1XrhL9XHKL0tN25E9mbZF2+iEAgTwEfpMUx+ALB+QrLX50rQRXrf6naK60TFXMTUuUtnJz7z8HLrDW03YkP7O2SDv9EIBAKgJ+OMAxeB2Qa5Q5gjtwT0sJ4u5vVXrajnzKrC3STj8EIJCMgK+Qa4pfq/XbXNPir2m/27panPYvsd3TdqQ/s7ZIO/0QgEAyAlUBWQH3nXR/Uf2i6Ne6pyteaXld2lrUPW1H/mTWFmmf9Jebj+U/nkl3983M+tB23McjM7caj5rov1ytVk9l3bkirhVgds4FbwKwkwt5Ttk38Z5pcW6LO9VNS0/bkWOZtVVo9xW+i+ew/AfWY+nPwEf55T/CXYs0pNWHtuM+Gpm51XjUQr/O6fQQb1VfVgfkGrHsAwEIQAAChxHYDshVUxaHnZ69IQABCEDgGAIE5GOocQwEIACBBgQIyA2gckoIQAACxxAgIB9DjWMgAAEINCBAQG4AlVNCAAIQOIaAA/LjzYGlPuY8HAMBCEAAAscR+L4c5oD812aj1KWPGgIQgAAE2hP4o5hgyqKQoIYABCDQmQABufMAYB4CEIBAIUBALiSoIQABCHQm8OgQ+3rFz7ksXJzL4lrLG7XtzH/hHZcqPW1HPmTWNrr27Gylz7lgnHPjRuvO85KmZNYWQXqo2qsDsgA5AZGTCX0wTNVORvNZ9a2WpkFZ5+9m277uK5m17dPtvuzas+qTLn/232v5osVpaa+0pCiZtUWA0F6ZoF6gnI3oiep1MDZYrftqwNvrTPdua1F62o78yaxtdO2Z2fqzr+WlFn/tzi8R63P2Z9YWcUB7ZUAWyJdanJpxWj6p4blA+oqhVelpO/Ips7bRtY/MNmJPPwRmCdTe1HO+XP97Ni1lqsL9rUpP25FPmbWNrn1kthF7+iEwSyAMyJVXv3+bPfuJjT1tR9Izaxtd+8hsI/b0Q2AfgTAg6+ASbD1nvKu0mrLoaXuXr6U9s7aicVedXXt2fbu40g6BkwjUBOQaA+X7pWr2XXqfnrYjXzJrG137yGwj9vQ/UAKP9O+hb9Zd7vF/bu647F6uZPxccovS03bkT2Zto2sfmW3E/t72b6aafpeDh/zH7KdV7npD6aldtv39levvsHRAfqoNP+c7+yWn6vcjPuqehVzAl5t73m+x0tN25ERmbaNrH5ltxP4+93vc5N/NiD721C7b//6S09opi98E+WoGdLlCdn+r0tN25FNmbaNrH5ltxJ5+CMwSqA3IfjXUbyRNi/8a3m3+ukz7ltruaTvyIbO20bWPzDZiTz8EZglUBWQFXM9vfFH9opxF656ueKXldWlrUfe0HfmTWdvo2gdiW24ulv8WI/Tn7M+sLeLwILU/iqhs9ftq2MmEnqn2TTzXP2j7HBPyPW3Lzb0ls7a9wtWZXXtaffrc+wrexS+wuPyqNt9L+ah6fYNm3drhR2ZtEY6Hrv1ytVrtvakXAaQfAhCAAASOJ6A/Qgff1DveGkdCAAIQgEAVgao55KozsRMEIAABCJxEgIB8Ej4OhgAEILAcAQLyciw5EwQgAIGTCBCQT8LHwRCAAASWI0BAXo4lZ4IABCBwEgEC8kn4OBgCEIDAcgQIyMux5EwQgAAETiJAQD4JHwdDAAIQWI4AAXk5lpwJAhCAwEkEDsllcaFX/N5srDmXxbUW57Zokgt56lVP21Mt0+3M2qZap9vZtQ+gz6kHnNfiRlqdDzhNkR60NRqNVmyrA7IEOIn9z6o/2EfVzvb2WfWtlqZBWefvZtu+7iuZte3T7b7s2rPqky5/9t9r8TebOC3tlZYUBW3thuEcbKumLCTEyS+eqF4HY7usdV8NePutt1uVnrYjnzJrG117Zrb+7GvxVw/9KM6/RKzP2Y+2drTPwbYqIMvFl1rm0mx+UvtzCfUVQ6vS03bkU2Zto2sfmW3Enn4IzBKoDcjO+ep/z6alTFWUnLDT/iW2e9qO9GfWNrr2kdlG7OmHwCyBMCBXXv02+baEnrZnaW01Zta2JXN2Nbv27PpmodIIgQUIhAFZNkqw9ZzxrtJqyqKn7V2+lvbM2orGXXV27dn17eJKOwROIlATkGsMfFezU6N9etqOXMqsbXTtI7ON2NP/QAnUBOS5ueOCq1zJ+LnkFqWn7cifzNpG1z4y24g9/RDYSSAMyJrPK1MVc9MSpa3c3Ntp6JiOnrYjvZm1ja59ZLYRe/ohsI9AGJA3B/+m+mrmROUK2f2tSk/bkU+ZtY2ufWS2EXv6ITBLoDYg+9VQv5E0LTdquNu6opn2L7Hd03akP7O20bWPzDZiTz8EZglUBWQF3Hc6+ovqF+UsWvd0xSstr0tbi7qn7cifzNpG1z4Q23Jzsfy3GKE/Zz/a2tFuwvbRAXp9NexkQs9U+yae6x+0fae6delpO/Its7bRtadlq8+9r+Bd/AKLy69q872Uj6p9AdOtoK0d+tZsL1er1VPJd/Ke680Hqp03nBkCEIAABL4ioLjrXEFvVV9WTVl8dTQbEIAABCDQhAABuQlWTgoBCEDgcAIE5MOZcQQEIACBJgQIyE2wclIIQAAChxMgIB/OjCMgAAEINCGw/djbP3WXb2rkg9qcKJwCAQhAAAInElA8/adOcbXrNA7IfnbSX0UzV5rkqJgzRBsEIACBB0CgfFH0rKv/Dwa+71bTrHEnAAAAAElFTkSuQmCC\n",
110
      "text/latex": [
111
       "$\\displaystyle \\left[\\begin{matrix}1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1\\\\0 & 1 & -1 & 0 & 0 & 1 & 1 & -1 & -1\\\\0 & 1 & 1 & 0 & 0 & 1 & 1 & 1 & 1\\\\0 & 0 & 0 & -1 & 1 & -1 & 1 & -1 & 1\\\\0 & 0 & 0 & 0 & 0 & -1 & 1 & 1 & -1\\\\0 & 0 & 0 & 0 & 0 & -1 & 1 & -1 & 1\\\\0 & 0 & 0 & 1 & 1 & 1 & 1 & 1 & 1\\\\0 & 0 & 0 & 0 & 0 & 1 & 1 & -1 & -1\\\\0 & 0 & 0 & 0 & 0 & 1 & 1 & 1 & 1\\end{matrix}\\right]$"
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
      ],
      "text/plain": [
       "⎡1  1  1   1   1  1   1  1   1 ⎤\n",
       "⎢                              ⎥\n",
       "⎢0  1  -1  0   0  1   1  -1  -1⎥\n",
       "⎢                              ⎥\n",
       "⎢0  1  1   0   0  1   1  1   1 ⎥\n",
       "⎢                              ⎥\n",
       "⎢0  0  0   -1  1  -1  1  -1  1 ⎥\n",
       "⎢                              ⎥\n",
       "⎢0  0  0   0   0  -1  1  1   -1⎥\n",
       "⎢                              ⎥\n",
       "⎢0  0  0   0   0  -1  1  -1  1 ⎥\n",
       "⎢                              ⎥\n",
       "⎢0  0  0   1   1  1   1  1   1 ⎥\n",
       "⎢                              ⎥\n",
       "⎢0  0  0   0   0  1   1  -1  -1⎥\n",
       "⎢                              ⎥\n",
       "⎣0  0  0   0   0  1   1  1   1 ⎦"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "M = moment_matrix(moments, stencil=d2q9)\n",
    "M"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
150
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABbsAAAA0CAYAAACjIeUlAAAACXBIWXMAAA7EAAAOxAGVKw4bAAAgAElEQVR4Ae2d7bXkNNaFL706gAYieCEDYCKgJwN6JoKGDGYWv3r+9YIMgAj4yAAmAqbJgJkIgM6Adz91pWqXrz8kW7Llqq21XP6Sjs7ZZ1vHlmXVW3/++eedkxEwAkbACBgBI2AEjIARMAJGwAgYASNgBIyAETACRiAHgX/9619PlP/zUOa9sH6u469z5Bwx7y3b3rK/HgfHPNDxFkj5wGgfMAJGwAgYASNgBIyAETACRsAIGAEjYASMgBEwAkYgFYEv1If4Wcys7a+0/UrL+/HYFa9v2fbd3TrWp/1ImkHAP3rLF7trbAWMgBEwAkbACBgBI2AEjIARMAJGwAgYASNgBIxAAQTUMfajlg8KiLKISwQ+Fa5PO4foU3zvRrC+Zds7Ll++KZ480cK1yRcCuel7Fej3aX//OEj5q4T+lCvR+Y2AETACRsAIGAEjYASMgBEwAkbACBgBI2AEjEDLCKjPi06x77X+pWU9D6obo7r/c1Dd16p9y7avxe5UXtfkay3/1M4rrT9kP1Ww8v61m1f7n2r/Wezs7p7zthEwAkbACBgBI2AEjIARMAJGwAgYASNgBIyAETg8AuoAo6P7v1p/fXhjGjRgAFc6gMH76l8s3LLtJakIV7ScOrwld/X0N0xj4mQEjIARMAJGwAgYASNgBIyAETACRsAIGAEjYASuCgF1oDHSkyk16EhzqoyAcGaamE+0fFi5qubE37LtJZwh/H6QnB+05uXUquTO7lXwubARMAJGwAgYASNgBIyAETACRsAIGAEjYASMQGsIqNPsPenE/NHPWtPtGvXp4J01FcU1YHHLtpf0n3DkpdQHWvOSanHyNCaLoXNBI2AEjIARMAJGwAgYASNgBIyAETACRsAIGIFGEfhRer1Ux9l/G9XvatQKnb3/1Po0h3LYv9P66rEPtt6k7ZUITIf3N8L1Oy3J83d3dXFndxcNbxuBAyOgRuArqc8/2PLph5MRSEJAfOEfj/+t5eOlgSSpImcyAkbACBwAAcfSAzgpU8WjxDlzL9Oxzn5GwBw/Q+GNjRFonXvSj5Gh72j95cbQ3Fx1wpgR9PRH0OHLNCYk5u2++qljbtn2k5cr/AhTpjL5XKL5KgMeZae3Xrx48atKfSZBP2WXdgEjYASaQEDXL4GFQO7Ps5rwyEMl5JunOop/eDPJzQA3Ak285ZYe3JAwL9bNfW4mm52MgBEwAicE1BbeTCydikk6x0tQOgj4c6BX2r/4My/tEy+ea71opI3Kbp6ka9NxTvrdDPc2d37BCuUn38stxPOWOD7FE51z+7qQQ0uKCe8/VO5rra++w3UhPqNtWi5XA9bw+yLp+FsXBxrZkV43afuU3Y245qxG0JUvM97XdnK/ifJyD/vs8VmSN4yAETgkAuFi/kjrm/sDiKM4TL7hjeQTrU9vJbUmuL7S8nYLNkgf/vkYHb/R4hcmLTjFOhgBI7ApAmoDuTG+iVga2vupmPS58vBClj+XIi6cO7sDTp9ofahYIX2bjXMB05vg3qYXdeHK5Cffyy3E9JY4nsATt68LeZRbTL74h8rQ+foyt+wt5C/NVclr4rk2xXe3anuC3SnwbZZH+v6k5RdVSPzNvu/0H1Ru5ipXZATKI6CLn5FKXPwfl5duiSUQkI/oLKBjoPv5zX90jI4GOr2bSNKFzgx0Ql8nI2AEjMDNIKB272ZiaWjjR2NSwOLn4Hzm3Py9RwSO8eBxkVTui7D8g/XFyUZ2pFdzcU463Qz3GqHBIjXkJ9/LLULu7u6WOD7Hk4DFkvaV+/NPtTBQpskk3ZprXwUUz15MhXCYr5C2cq4wmWzTdJ7YlM3VrfRfU8+t2j5n9xpMK5flZRX3rQ++Gpir153dcwj5vBFoG4HTiCtd/A7i7foJH/Uf/JnGhJTdaN8Xq/bLG1P+CKI1vaoZbMFGwAgYASFwS7F0Lib9VzEg/vfH34QN02t0Ey9pv+0eUH5Gxf+mNaPBmReV/w/5sZunoe3W4twtca8hGmSrMnfdZAusWMAcrwjujOg5nixpX+l0pI3l3rz1+/NmuKcYBG48b3lUt0AYSMW5OlBHq4du1fY5u1v1V5xum3YwK3kakyy4nNkItIOAgjhvZAnkq0d1SxY3T58H62JH7GHm42xVf+lFpwDYfhewjas4ovvB6LiYYY+19H2tBV0JhtywOhkBI2AErhoBtXnFYmkOUKp387irOmdjkvK8xo5O3u4UJtwfoHd88CAr6eILM5Xls1M6vN/TkjzH4r2our/Sp5k4J12Kck/yNudUSW+1qr/0mr1uSuKwVpb0vUqOt8qP6K8UnuAb8nfyzravysuzAtMw0V40naRjM9wTUDzXok/RZy3JK97O1pA5RRTVN9umKU82V6lT5eIAr9+0+672m5orXfrUtB1uMEiA/0NsanrZFLuld5NJunMdc9/JlxpZfzT7uEmLrJQRMAIpCBBMvqQBSMk8k4fPj2lATknbjOTiUzn+nOoIqVX96TDmBrXvo7/rOKM7muoECI6GV79Kt+Y6KYJ+XhkBI2AESiJQMpbm6LVH3MqJSeQlTnXj16mzRcfOnQfECuXjAa8fzyjHC/n+cR3aPbUS50pzbw9OlXRmq/rnXDcl8Vgj6xo53io/op9yeJLUvkbBB1u3wj06NfsvZktAWYOHNWRO2VqFq7ofiF95nToktf9UCy++mf6slVTLdu534ByJe6LWUo7dremOPvwx+lfiUlb/xKMWLbFORsAITCOgC730p1mfSmZsoKmcGxUaE+o5QmpVfzBlfu5zCpiCK0GnuST96JiIb0+b088KGQEjYARKIRDaYzpr9/jMeY+4lROTwOXcqR0w54H11Hkg7NCfB7qxh7rfde6dUK6pVQtxrhL39uBUSd+2qn/OdVMSj8WyrpTjrfIj+imHJ6nta5R9mHVD3CM2/VgBuBo8rCFzyvRaXKX/4PyCQVxgmw5v+N5KqmK7bGRwG538Lb7gB/scu1vxVVcPvjwnYUdyepScM2RsjKy56lfL3wIuLehQDeAGBe+MN59m8Znw60LQMKr7olO2kNytxDSnv3zDTRbB/fywH47xZpLPm/qdCFthlVIPOmbPi5UiODWP8GnpxihV7UPn2xvzves/tPMOqnwDPi8dS3M8sWncEta5MenigU3leUn7kZYYu97Xsal7EGIfdbaa9o5zNbi3KacqOLY5/RdcNxVgWSzy2jjeHD+iZxbwZG37Gqtudb0392KH2LnjtSBQNXhYQ+agybW4Krk8txHzL7itfe4TuH/YPdWyfXfDZhRYYPeMxO1PywZ4xJL1lcDjHFVVyT+UHwL3SZwj5lrzMgqWfwnNmkemFBj2TSkks+Ts6XM+JSYwFkniz3nOuCAQ2Xy+HB9qi9RTS0ij+scbLf6wi7aTxLQwzw6AKzeHT6QnbVr8ozL03yQFvBxrNkH7opLd2jT7/MIPt7SzG+cCyEVjaY7jxPmt425uTGKeTf6wmJFazL1Jm8wclHxuTUz7VguJh4+h9EQHW35e2DXOCZvi3NuBU0N+X3ysUf1zr5vF9lcoeFUcb5Qf0W25PEltX6P8o6335h4dYszzWzwG1eBhDZkThKnFVWL+UGrpK69atg/Z3dKxXLtb0r2rCwMzs16cPO6WntrWRchNGZPMD3bm6jgV8xbvQ22P3fhOVXHoc7KZUbYfaOEzlP4DTFXbVN+kb6pWfsPChfsuPle98SLnRqJ4CvLhVFN/rJBqaEP6c6PFJ03caA22m6k2bZ0PnYPe2LBpZ7fqdXu2tcNDfcJ+rzbNPt/J53tXuxfnsFt1V42lOdgGXWrH3ayYJJ24lx+abuvimPLFeb0Z1dV/Qd7fz4Glat6gN/F5jzhXnXsbcaqajxrSP+u6qQbIAsHXzPGG+BE9k8UT6Z/UvkbhR1vvyb2A1UdaV/9iuQYPa8js8Wdrrr6j+sc6wnuqVd/d2vbqBiVWkGV3osw9snFPGTvuk+p/lJJLFx0E/Vzri39T5biW77XwZ3bfaLnpz86FAx1aTE+w2QUd6nrgG+nhtAECe/hcZhHA71R3jbfVXMOM4jrkSyth0pL+NMZVXkhI7hYJ3U9c26Iy6pD/BmPNVvW7npMPNo1j9rlZJw5syrkO4tViaaeO2U3Zv1XcqhmTXsrQ8wOIbGL7B62L36fMApqXYfM4F9Sryr0NOZWHdmLuxvSved0kIrIq29VxvDF+ROccnSfRjpLrXbgnfvAswVI1/tTgYQ2ZAw6txdXXA3VxqLovRuodOlzL9qG6Wjp2LXb/Cqi6TrAnKT1OynXf+UWH9kVSRZD6NMJD23zSGEcqXOS7sR1worOw2BQTM/hR1wPfzJTx6bIIbO3z09u5siacGg4euJlyA/l3WrPPuurNAnWUSEHfIvpLFo0on9cv+kpD5Qjs4PdzCdt2kkFA2XrebrdnOzm7V+2WbZp93gP/Rne35FyEuEosjcJT1iXj1lR9tWOS5H+phelNuJ6Z8oQpu55P6dTIuT3iHKZX45584Hu5QC5h4Xu5u7ur4niL/JZOVe/5g830K0Q+n9pZHW/9q9G9uHd6fhVe1F8l1eBhDZl941VHNa5KdrWvvCQ7cn/Rczk41LS9j3PJ/bW2H9XuEQzj1xr0OfMybTY9ns1xn+FvAmqrzttEldrMJpy+1vKHFjrdeBlQO9k3tRGekb+Dz2nwky7wGdXPp2UDNwZ0NsDb+NKKa/7ia45zgcY2KujPzQDL0hRHRBf101JlFpY7veSAD1q2+hTd7dlCZ5UsJn9vGcfs85LOO6isjTkXUSoeS6PglLVs3jLuVo9JsucQ9ws93+wR51ChCvc25lQPyvW7FfT3vVwY3Spst7yXq8LxCvxYT9p7CVXbV9lNO+X2Nd1bsbO7ymCtGjysIXMErqpcVZ0vtRDfTs+NsovtEl95rW3Lpcb5a+WjPZuvtb22z8F2qxSvaQZUJKXHc7lE0k+UJwqey+7z9wiA19+0LH77lAKkfZOC0mZ5NvF5sIZGr/Q1+Uoykcv6nMSxo7zkak1/3jz+Vfi9PoN5vI3IMW4aq3d2uz1rjiDV2zT7vDmf761Qdc71DKwRS3tVTO5uGbeuISZNgrnw5KZxrqNjLe5tyamOOcU2W9P/Gq6ba+J4a/yIxL8GnkRbSq734h5zRJNi/fd75X5r8LCGzCGLq3JV9/Utf+VV23ael+k3oYOfr8Nb+QKjqt1DJKt1TJjyp7OIj9f4bFWznd2SwKd2zb0BkaEQiSlU6EyCXIxILdKoSc5XWkY7+XSO+n7U8kzbQ51A4AVuVTu7Qx2DvpFe3EgzBQFvPl5p/6xLOPeN1hd/MqR8RZLkjvom1D2oF5XrPH9y+lxr/LpZUn1H8XnEZPbTLNmUzAHlfTsK3mqtOovxZA/9p3CSPvB38NqcKlfinOpey+WoRmxPkwNKLLhwPRlrZFfkM+L/ouW5Ftriv2sh/aw8P9xvlv2V3GJcLavZuLQCPNgijhXzeYcfVx/zpmyFETq/SxwNda9pf7bgXP+iyYmllJ1se6Z8E86d7720v1ncVV27xaQ+4KX2ZdMarkU1to5zsV7Ws9wjk+xMjn3Kuxmn0I2kOovFxz30v7di+Ff67HbdqO4S/MawpjkuO4vyuyNvs3uBPXkyzNx1R6+Ae3CKxPWblDq8If9cnE9qZzsyZ7movEkyk4yZyKR6qrdpqqPJrxBq2y75tLXN2V7b7gm63anuUnGsX028xvvHH+w/enDk4QGGvifdkD0sWueIgONNCR3N/BkkpGL6Bd6IrU6Sx0h2OrLnEp0s3OANJfDifO005ZvPZQtzeWELeHUTo86xs3hSnXO+GdVLZekE/0Tr5OBUwgDVdxifS9fIq98TbB/FWmWrcSBBrzvZcTiepNi1d55CXD6ZIVnxOkwOKCvtn2rPEM18sIwYoF1jLnT+FPmp9okB6Nhv53RofZL8w3FVOh+lTSvp813au534MWqr9NkljnKlFeDdVvdO6JoTS3PanlHfCKJd4+761rAdCQW4djJGcraOc7ncQ88c/m3qJOF3uPi4KUALKyvFb6qXrNY5XprfboMX8i7wpcT940mDPbgXTH83rFOelUPWKu2suRjR9frmECgZx3rgEdOSB+I97hUe2qUTIaexGJJR7FgAjg5R3pLFxPD8JzpGx8fa0ZRMPTA6qpsKdZ4J+OlwiTcQHO4m8IoPUt3jpbcHfSPdPlBF8Y/xGDnX99/kaLqlSqpeAuSob4Iec3r90q9fcrGTh0RebnzYP19g/2g+nzV5Lw7MKqYM0q0KT1LqvoE8JbjchyneNPaPl94fbM+oRJyhA6/bmU3bC48Y3U0i6HXPnw6u/anFVcmNuv4mHd/VfumRACV4sEUcK+Jz4XdtMW+UH4m27hFHudzW8m4LzsVmAe7NJuGd3PYk+mbtPeqszjeSYS3X+jBtFeeoN4l7ZMzhH/m3TNLN93L1AC/NbzRtjuOl+S15u9wL1KPBLpKPzj1AS25jyVyah0GmuQgQTreMQI22BDx5Vki+xlM6u+lEeI3kRhKj+fodA7FjOdnwIVvU2DFSm89/z0nHkPm91nQQdxMjkMYeWviMYZUu3Yomtsd8Q2d8fNilk/hlTwZ29o/1sizanfPNf3L1Un6CBfqSimMq+Uf0OVi85mci7cWBCZXOp4rz5Cw5c0P+56uQyK9uaa6tO50fevH1i45PTgGk8392hU1tK+9bU+dTz0lOKS73qyx+3fUrCPtj7RmnaTvip7js89ILP5yuA60n/UGBhak4V6UrnWe/ac0L0zutn2r5UUs/xixSGXkqeJQ4Vsrne7V3e/Aj21ZxomocjTwuwLut7p1Q+dTGaz0XS3PanmzfoMgWSRxIjkk19FH9ReIcuklWqTaua+pWcY46U7lH3hz+kX/LVLz9W6q8OOF7uXnwWuR4aX7v0gaLf1fRvlZqW2HmltyjvlMbK3vm4jt5SaV5iMxduEjFU6l1rkq/Km05mNyq7Tl2K2+RezXJqXGf1qV2vI/qHhvcfjx4tNGDATgazO96KgIoKXbw3u/l//IGot+R/pHEDDWW7yvveR7sXlU4YKhML1udXel1qruD11lPHftAtYLhWEf9IqU6dY36ZkYvXlg80Etl8CmdWowgqZGO5nPeZpHAajTNYF2FA6PKdE5IL65VdC/KE6qQbORmfQGgMkOd2ciCb/y5xKlDEvk5SeWKBIucOpW3FJczq62fXXj22/ahl3hFFVGdtbjKqN2Po7Kq5yctdHbDt26HfsySuy7Fg73jWLLPhdtrQNI6+uywMU9mTPJjxta94ijwl+DdlpxLjaVN8xDgU5J4s0dMSlFtSZ4SXFtSb6kySdyjMvktmX+llEuRI71iW1vjXm70y5Yx3aSP7+XGwNnneBLHS/Nb8ja/FwBe1Xst7evR29ZFbC/Nw8CJXbg4B0DrXJV+Vdry4JOmr9Natu/k82bakkdzF4XOE7DoRGohMXrvPJqvo9Dftc0btLUdBTwo9hMdg992D6qeoXzdLOAVA333eOntOd8M4cUNKv9k+kthZYbqoooh35AXf50CQdDj1JldQa8gfnQ15MuWfT5qyMiJIb/U4sCICheHh/QhwyqeiDf4jBGzXHuttFdSZdNUist9pbf6z4a59uykl3wdH7LPf0apY0xjVdrvxbkqHfERevZj1Wsdg8MlUikeoGftOFba50M+q9XeDdWF/xa3ZZn8oP5W4ih2l+DdFpxD10VJ/klte4a4UYuHi2w5eKESXOtDsFWc69ebvJ/Bv2SZKzIOcRxxi9s/CsvG+OXTP7XNYANeBP/IuRtKNfgNfE1zXH5ObV9TqDDET7fB88hdC/dO967iFPcUWakwD6nbXMzygDNfCQK12pIIT/Lz6aNYYmLNQ/mQwhNFHp5a0uA8lHJHoGJ+7nOSXDoIWGhMzmlhfXSYnDscwjYjDM6flwe5fCLOZxZj6R2d6HdmPMi7UMeunDnf4Le+HtjyYFR3AV2SfaP60esXLd101ku6fJqrT27+TsWb+XyFjh117+gUI8GxlJTMgRRhBWyowhPpxUswHoz6fE8xq6k8KzAuxeU+HpFz/eMX+yv0jnIG2zPkauGBF+6QeOv/WvtdX/MnMBd6aj/7Jvck/c1PDa6O6UTQvrimV+hfigezcWyFjhHloj6X0OT2roDuu/Ij2NpKHMWfJXg3yzkqKuA7xMT24uK640RM1KMlu+1R+WQexrqm1oXsnapi83MrbSrBtb7NkQ/94xf7K/WOsmJdo9wjI3VpWcK/WM/ouoAdNdo/9OWZ6/yMIj3Zfqo119Sh0gqMa/Ab7CLvJnFcoXdXbqxrlOPUo6UKv6VIchtcyN6u7btvr7DpGrgH/rP8i06qzEOqSeZi1GlqvcK3U2J9zggMIrCCb7XaEvQkrsRrfFDv7sFH3Z2RbR6m/jJyrns4/vHFg8AWgPpD61fdAjnbQQYNxll+OEbHM39ceH7oC8eX1IeMf6s8/wyN3H9rYW5YOl8JyHRw/0/LV9rudrTo0EU6zSd7caS3o/JPdGiJjl1Jc7650FF1cnPKcjFKYq0uoXySb4Lyfb0+0HGmi4k+ZIqYZBKv1H8Tn6/UMcB2Wv0eduBPSupjPciBFEFrbQjld+NJio1751mJcSkun2CQLviKFDl3vzfwu1LvKBH9h2JN5OzvoZ4LfXSM8z9HIazX6hPKb8lV4tr5ml6pfykeTMaxlTpGdxXzeRCY1N6t1X1vfozYumccRaUSvJvkHJWs9R0yQortyPm6iyc66+y2J5RN4mGnntHNgvaO1rH1iQI2leDayWzpsnWco94U7pFvKf8oO5rW4h/KF4+PwRdcjxfXj/Zfa6F9O0xaiXExfgNYwJXNyDu2B9NKvbsyY1012tduPWPbFxySXfFaKvr8O1b5nsdX+vAauAf8/Pk7aYp/9zkqtbNRuNZJXOzkH91c6dtRuT5hBIYQWMm3om3JgH4X19XA+fOhR+et8Q2m8Bi9yRAQ/HkjHcN8ekY67etY3L/TNjcqKEUvf0rDg5x+IlCR+LTtHyzaZgTAM22f5+ckg/aX1vdSxRk5jmxu5D6ULJzFSEI6Y0/zxPbr0/F+QteLgNrPsELHrqhJ3yjjP8msuui8x6Y4+v08aiKcX4oXxUnJvrnP/kCvE9Y6x5y16IldyUll1ui/ic9X6njGIsjB3vfPB6c3kjgwLeL+bAEbduVJio1751mJcSkuRxiehA34NplW6h1lj7VntFe08fDnU9VFe/xca146nqau0fo8pYnO3Wl/TZuAiFpcHcMSrM+Be6X+pXgABqNxbKWOYEwq5vN7cQ9iy9FiXhI/RmzdM46iUgneTXKOSgrxLsoB76lYmt32oKNSS3H3XqOGfgv4sATXIiK0vaSxa+/+rH4L6H2SFeTMcY+8S/l3qmfsp4AdteJj9EVfdTpOeSF8mLQS45L8BrOIa2scr8LvQJKkNniln5rk40qbDs+94JRZrnecV5OHVJPExY4+o5srfTsq1yeMwBACK/lWui3pqkhM474gKb314sUL5vBiZDQX+2DSOfLQqUzH7+Kk8szLzJ9x5TRCp/pUhlHVH2nNyJ+ktKa+pAoGMqlOHjgZBT71AHUuuVZHlU/2jfJ+oYqfaj2I4VJdVC7bN2cAEjeCbnTaj+K6VP9EFUazqd5kn5fQUTL4QoJRrufpdUaV651QmUkO9LIP7i61QeWa4MmgUb2DwcbFf1AZxUkO3Pgs7NM4sw+PR9tb8ob6F7WVlC+Rgg68yHxb20lt9lq9VT65PUuxcak+KleNq5L9h3T/WOtzPNX2nzrGFy3nDm/s0/7imEn5pUn1btamqa6iPu/aLNmT7d1SfFWuCX50bU3ZDvYePo5i61LfdXGSjMWxtCtnblv1TPJwrjznS9ibUk83j+pcFL+6Mqa297Cpr0/QYdM4hw6qdxPu9e3t7i/FX+WqtH+Sy8AqcLm459BxYuZLrbP/MFxliKG+l9NAtD6u2h9NAbdV96CSsTvHo4HSZbINLmFvrCt1rTrdvg6AVcIXQQac50/yJp+3BlSoekj6THIxpfISGKXUE/OovqpcjfVMrYPNq9vyqTr651qwG51K2i5ZxO+YeInM4LHJZ/xQ/6p4ECtcu5Yu9KVwT/Cltk8vksZk6jwD4p49GsvQO86FGTtseqeydv+iiicBnZDGSILcBmtNfROqTJ4CePBKTWt1HPSNcOaBlhuNU9I25MDpvGkZS0t1WeKbMR3WHF+q/5o6KZvj8xI68vUBXxpMpoUcmJQZTi61oRWepNhIO7W0rTrJD9ccX6LEhTaUIMPLMB7AptJSjKdk5p5jSpH+H+DNyVir92B7NlfpxPml+tTkKm0w8k9JXGD7B60vOrrD6aX6h+KLV1u2aUV8LvyuJebl8GOxgycKHoFzqF9Cz6RYOoHVg1MLefhAzsCBEvYOiB0+JDu4Z4yxi3VO/BoW+vDopjY9rP50ZI84R8XFuTdi39ThpfjXio9j91xwcSg+TtkWzyFzTG7MM7leeS0sxXhSp8yTN8XxhW3wpn5ayalU929q04hSe3Evthd00u6WFnIxRd/NfLsRV1NsXt2Wp1QS8zRkNyqtth17tDC46FetGeTMPR3PG//T9tx1shnfpM9cirr+Npcxnk/q7BYIfELO25RYQSyfvFbZJ8qcrFhXcChL3T93j09tr6lvSu7UuYAPOF1MqzJWpoSOoa4h39CZ9m2n7m+0/bXyX3zuH88v1SWUy/JNrDNlLfnYRmfI51pO29pnqpOLFPRYxK8LQZk76KciST4vqCMvMVKmBMriQIrpS20I5arxJEX3nDzSlzeYSdfxhFxeLjHtRrdjO76wg8+DKWC1OZcHlGGU1Xn08cD5i0Ml9A6YD7VnF3Wl7CzVJ5SrxlXJZ3Tau1rTORunl3retynosTkPVO+mbZrqW31/EbDLau+W4hvK7c6PPl+m9vGplquIo1AB/VMAAAoJSURBVNgZfFDi2kiNpVPw9s9l8bBfeGi/oL1D4seOLYpfY8L6x3eyqa8G+5vHuaBEDe4N2Td4bCn+oVyV9k+y6aDioR75/ZR8L9ItKJm+l7s9jme1wUuvhS7PFmy7fR0ArZQvJCe2F4Nfsg9UXetQFhdTlCiFUUpdIU9VrqbqIbtLtOWp1ZGvCbtRpJDt9AO+I1nnL6S0zXXCi/fuaG/tvknKs7j/9o2Uolvx/iC+0JoV/ng2x5sMz7R5+iTkzaGsLTp9zgBnlXwzijV2FKUUX1NfivyhPJCFNyWpqZSOQ75hZB6kpjOFqT+Y33YKv6W6fBSMnZIdsuSvpDNkxpa5tFT/Oblz53N8XkrH76QU9Z5GhE4omMuBCVHnU0ttqMqTs3ZtbRBELt7Gis+vtcxpuRTjObm55+FXK+1Zru7kX4pjda6KA27TLj06FMMuc8zv5bZ3R+fHPCIhx5XFUaxa6rs+ZqmxtF9uaj+Xh1Oy4rlS9kZ5Keul8StFNnn2sGlIt73iXA3uDdk3dmwp/rXjIyPN8Mmps0ptF9tjXz6N2Vb6+NJrYSnGpfW/NY7ntsF7+Gkpp1K5sYdNQ7rtxT10of8gtldDum1xLJeLKTpt7dvaXE2xeY8812b3JwIRm/qJY/wXIoMo6bPop6351q+/vx87u4ds6ec97SfN2R1LCgQq+ETrpZ3WUVTWGgeoAPN1V+lQzVJmJLN0o1N5txsy1W/fjPim1uE9fa66+eM4ppjI6YysBcWsXOlZ9RqWfPgPFtxYMVKLNuo3Hd+0rVKdk0n6EGx4acgnRGtHjk/WteakdANHOPZgHuk1clPLqv5d2jP0U91VuZqKwR75ZPtucWxPn+dgfcv8yMEpNe+enENH1X+oWJqKa418wuoQ8SvVdtmzd5w7HPeEWfX4qDr4CoXE1xsM1mEandccaCVJn0NcC9LTHG+FNDN6HIVTM2acTzfAPZ61nkqPt89KeaMIAtfG1VRQjmq39CZuM881o+Mv/u9N+zz3EXObm99eOj1I0pfrmr7otx6c7B1QHkbnP3vcOz65q0K8Jdu880j1cpPTbEc3oEnHzXGh3pjsm4jEduudfc7FTuN0lM7uqtdw4D9v0JtN0pGHjpPPtN1sR3cAkGD4S8B1c0xDvbu0qaq7Klc3BzOjQtm+C+aouKfPMyBCz5vlRw5OqXmF526cCzoeKpam4lo6n/x0pPiVav6ucU5KHo57W7R/qsP3cqkMns9njs9jtHsOt69VXMDLRDrFxkasVqn02oVeKVdn3XZku6V7/KqcP6Tsp3fDAQaZHSFxL5rVJ5zV2X0EBKyjEbgRBL6TnUxNw1vrrIv+RvBpxkz5h5HmNM78UQuf3TA/VuuJt6HPW1fS+hkBI2AEViLgWDoB4EHj14RFF6f2jnPm3oU72t456LVgjjdMq4NyKhXRvbkXn7V4/hr8v7JUQ5zvNCDliM+yq113Rdco1wDXQj/hV9KT+1W7v/IFHfLoyYus5PQoOaczGgEj0AwCuuAZXciIuPi5ZzO6WZFLBOQrRkh/qYV5ib/V8krbfALbZJJufNJ0p7VvDpv0kJUyAkagFAJq5xxLJ8AUPoeKXxOmXJySXbvHOXPvwiXN7xztWjDHm6fU3dE4lYpoI9xjcBHx/e+pejvfOALXytVxi+/PXJHdpwFssufc4a1tOrq5RkjJf/h4n32X36h7Vv+EO7t38ZUrNQLrEVAjxaee74XGar1AS6iOgHxFA01g+V7bT6pXuKyCz1XMo7qXYedSRsAIHAwBx9I0hx0kfqUZc3fXRJwz91Ld1Va+g1wL5nhbtJnU5iCcmrShc7IJ7kkfvp6JHWQd9by5BoEr42oyFEe2W7rT9/B/Wp5pmz+kZLAkI6V/1kI6Qmc303Lxf3VZurqz++Rf/xiBwyJAh/c3h9X+ihVXY/wBy4CJ/wnHmrsBk76MdiOQZL01HbDRh4yAETACR0LAsbTjrSPGr476k5sNxjlzb9Jj+5484rVgju/Lmbnaj8ipOZvi+ca4x/8iMGd3c89bEa/W19fM1Snsr9Fu2cTc3Z9p4Wtz/viZZ/04Z3dWB/IUdhXP8VX8V7ny3dmdi5jzG4GGEFBDdfqjQ62ZG82pLQReSR2mLGl1BPcFWtKTN7yMxmC6FScjYASMwM0g4Fj6wNWHil8PtB850GKcM/dGnNXO4UNdC+Z4O8SZ0ORQnJqw4+JUa9yTPvynFZ14fq658FTWzlVyNQGBW7GbQXk/6Vph5HezSfrFfq5Tv1eOou7szkHLeY1Amwh8LLV4Q0dnpVM7CBA4fhgIIB8FFZv5Y1HpSIc8f/jwXNtHeLsbIPTKCBgBI1AMAcfSN1AeJn69UXl6q/E4Z+5Nu2/Ps4e5FszxPWmSVfdhOJVqVcPcY7qGT4N+qeY43xsEro6rb0yb3Loqu8X/T7T80b0OwjZfPfB1WesJHYf6VGb1dmf3LETOYATaRkCNFQ0y8xi1PA902yDW0Y6G+eIfg+UrPsF5ooXPiPBbK4lP/b6QTp6+pBWPWA8jYAQ2RcCx9ALuI8WvC8UndpqNc+behNf2P3Wka8Ec358vKRociVMp9pCnSe6pbWUkKM9bfLnqlI/ANXI1BYVrs5sBkb/3DOeapT+CP3NtNkk/OuTRH59kp7devHjxq0phaDOjDLOtcAEjYATudA3TEDA3WdON1i25Sr6gge5+PoeP6FRupr2VLnS+f9SSTrfEEdtqBIxAWwioLXQslUuEQ/PxK5U5suUQcc7cS/XotvmOcC2Y49tyYm1tR+BUqo2tc0/6MQUCc/2+re2WBhqlQrxrvmviag6Q12a37OErBxL3QyQGSTbTH3Gv0sNf6cjAQeYb7/anPMzYO6L8XPfP3NndA8a7RsAIGAEjYASMgBEwAkbACBgBI2AEjIARMALHRkAdX8zBzNzEi0aHHtt6a28EjomArlcGXdDZnf2iKnZ2Pzqm6dbaCBgBI2AEjIARMAJGwAgYASNgBIyAETACRsAIjCLwXGc+VQcYX245GQEjcAwE+CLjma7bxV9kuLP7GI62lkbACBgBI2AEjIARMAJGwAgYASNgBIyAETACiQios4wpPl9qofPMyQgYgcYR0DXLtCu/aL3q/8Tc2d24o62eETACRsAIGAEjYASMgBEwAkbACBgBI2AEjEA+Auo0+1KlmPs3zl2cL8QljIARqI6ArlGmL3mqddY83UOKxTm7uei/62aQ8MXDxbtyvG0EjIARMAJGwAgYASNgBIyAETACRsAIGAEjYAT2QkB9XMzf/VLrVSNG99Lf9RqBa0ZA1+UHsu8bLR/n9Ecrb/zjzQgPf1D51ziym086/ugsfuMVYfLaCBgBI2AEjIARMAJGwAgYASNgBIyAETACRuCwCKhT7EMp/1noVDusHVbcCFwbAqHDOrujO+DwvdYP+rP/HwCKEnakCOFAAAAAAElFTkSuQmCC\n",
151
      "text/latex": [
152
       "$\\displaystyle \\left[ \\left( 1, \\  \\rho, \\  \\omega\\right), \\  \\left( y, \\  \\rho u_{1}, \\  \\omega\\right), \\  \\left( y^{2}, \\  \\rho u_{1}^{2} + \\frac{\\rho}{3}, \\  \\omega\\right), \\  \\left( x, \\  \\rho u_{0}, \\  \\omega\\right), \\  \\left( x y, \\  \\rho u_{0} u_{1}, \\  \\omega\\right), \\  \\left( x y^{2}, \\  \\frac{\\rho u_{0}}{3}, \\  \\omega\\right), \\  \\left( x^{2}, \\  \\rho u_{0}^{2} + \\frac{\\rho}{3}, \\  \\omega\\right), \\  \\left( x^{2} y, \\  \\frac{\\rho u_{1}}{3}, \\  \\omega\\right), \\  \\left( x^{2} y^{2}, \\  \\frac{\\rho u_{0}^{2}}{3} + \\frac{\\rho u_{1}^{2}}{3} + \\frac{\\rho}{9}, \\  \\omega\\right)\\right]$"
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
      ],
      "text/plain": [
       "⎡                                                                             \n",
       "⎢                         ⎛ 2      2   ρ   ⎞                                  \n",
       "⎢(1, ρ, ω), (y, ρ⋅u₁, ω), ⎜y , ρ⋅u₁  + ─, ω⎟, (x, ρ⋅u₀, ω), (x⋅y, ρ⋅u₀⋅u₁, ω),\n",
       "⎣                         ⎝            3   ⎠                                  \n",
       "\n",
       "                                                       ⎛           2       2  \n",
       " ⎛   2  ρ⋅u₀   ⎞  ⎛ 2      2   ρ   ⎞  ⎛ 2    ρ⋅u₁   ⎞  ⎜ 2  2  ρ⋅u₀    ρ⋅u₁   \n",
       " ⎜x⋅y , ────, ω⎟, ⎜x , ρ⋅u₀  + ─, ω⎟, ⎜x ⋅y, ────, ω⎟, ⎜x ⋅y , ───── + ───── +\n",
       " ⎝       3     ⎠  ⎝            3   ⎠  ⎝       3     ⎠  ⎝         3       3    \n",
       "\n",
       "     ⎞⎤\n",
       " ρ   ⎟⎥\n",
       " ─, ω⎟⎥\n",
       " 9   ⎠⎦"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
Markus Holzer's avatar
Markus Holzer committed
177
    "from lbmpy.maxwellian_equilibrium import get_equilibrium_values_of_maxwell_boltzmann_function\n",
178
    "\n",
Markus Holzer's avatar
Markus Holzer committed
179
180
181
    "eq_moments = get_equilibrium_values_of_maxwell_boltzmann_function(moments, order=2, dim=2, \n",
    "                                                                  c_s_sq=sp.Rational(1, 3),\n",
    "                                                                  space=\"moment\")\n",
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
    "omega = sp.symbols(\"omega\")\n",
    "relaxation_info = [(moment, eq_value, omega) for moment, eq_value in zip(moments, eq_moments)]\n",
    "relaxation_info"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "\n",
       "        <table style=\"border:none; width: 100%\">\n",
       "            <tr style=\"border:none\">\n",
       "                <th style=\"border:none\" >Moment</th>\n",
       "                <th style=\"border:none\" >Eq. Value </th>\n",
       "                <th style=\"border:none\" >Relaxation Rate</th>\n",
       "            </tr>\n",
       "            <tr style=\"border:none\">\n",
       "                            <td style=\"border:none\">$1$</td>\n",
       "                            <td style=\"border:none\">$\\rho$</td>\n",
       "                            <td style=\"border:none\">$\\omega$</td>\n",
       "                         </tr>\n",
       "<tr style=\"border:none\">\n",
       "                            <td style=\"border:none\">$y$</td>\n",
       "                            <td style=\"border:none\">$\\rho u_{1}$</td>\n",
       "                            <td style=\"border:none\">$\\omega$</td>\n",
       "                         </tr>\n",
       "<tr style=\"border:none\">\n",
       "                            <td style=\"border:none\">$y^{2}$</td>\n",
       "                            <td style=\"border:none\">$\\rho u_{1}^{2} + \\frac{\\rho}{3}$</td>\n",
       "                            <td style=\"border:none\">$\\omega$</td>\n",
       "                         </tr>\n",
       "<tr style=\"border:none\">\n",
       "                            <td style=\"border:none\">$x$</td>\n",
       "                            <td style=\"border:none\">$\\rho u_{0}$</td>\n",
       "                            <td style=\"border:none\">$\\omega$</td>\n",
       "                         </tr>\n",
       "<tr style=\"border:none\">\n",
       "                            <td style=\"border:none\">$x y$</td>\n",
       "                            <td style=\"border:none\">$\\rho u_{0} u_{1}$</td>\n",
       "                            <td style=\"border:none\">$\\omega$</td>\n",
       "                         </tr>\n",
       "<tr style=\"border:none\">\n",
       "                            <td style=\"border:none\">$x y^{2}$</td>\n",
       "                            <td style=\"border:none\">$\\frac{\\rho u_{0}}{3}$</td>\n",
       "                            <td style=\"border:none\">$\\omega$</td>\n",
       "                         </tr>\n",
       "<tr style=\"border:none\">\n",
       "                            <td style=\"border:none\">$x^{2}$</td>\n",
       "                            <td style=\"border:none\">$\\rho u_{0}^{2} + \\frac{\\rho}{3}$</td>\n",
       "                            <td style=\"border:none\">$\\omega$</td>\n",
       "                         </tr>\n",
       "<tr style=\"border:none\">\n",
       "                            <td style=\"border:none\">$x^{2} y$</td>\n",
       "                            <td style=\"border:none\">$\\frac{\\rho u_{1}}{3}$</td>\n",
       "                            <td style=\"border:none\">$\\omega$</td>\n",
       "                         </tr>\n",
       "<tr style=\"border:none\">\n",
       "                            <td style=\"border:none\">$x^{2} y^{2}$</td>\n",
       "                            <td style=\"border:none\">$\\frac{\\rho u_{0}^{2}}{3} + \\frac{\\rho u_{1}^{2}}{3} + \\frac{\\rho}{9}$</td>\n",
       "                            <td style=\"border:none\">$\\omega$</td>\n",
       "                         </tr>\n",
       "\n",
       "        </table>\n",
       "        "
      ],
      "text/plain": [
Markus Holzer's avatar
Markus Holzer committed
252
       "<lbmpy.methods.momentbased.momentbasedmethod.MomentBasedLbMethod at 0x7f3aca401520>"
253
254
255
256
257
258
259
260
261
262
263
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from lbmpy.methods.creationfunctions import create_generic_mrt\n",
    "\n",
    "force_model = forcemodels.Guo(sp.symbols(\"F_:2\"))\n",
264
    "method = create_generic_mrt(d2q9, relaxation_info, compressible=False, force_model=force_model)\n",
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
    "method"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Example of a update equation without simplifications"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
283
       "<div>Subexpressions:</div><table style=\"border:none; width: 100%; \"><tr style=\"border:none\"> <td style=\"border:none\">$$vel0Term \\leftarrow f_{4} + f_{6} + f_{8}$$</td>  </tr> <tr style=\"border:none\"> <td style=\"border:none\">$$vel1Term \\leftarrow f_{1} + f_{5}$$</td>  </tr> <tr style=\"border:none\"> <td style=\"border:none\">$$\\rho \\leftarrow f_{0} + f_{2} + f_{3} + f_{7} + vel0Term + vel1Term$$</td>  </tr> <tr style=\"border:none\"> <td style=\"border:none\">$$u_{0} \\leftarrow \\frac{F_{0}}{2} - f_{3} - f_{5} - f_{7} + vel0Term$$</td>  </tr> <tr style=\"border:none\"> <td style=\"border:none\">$$u_{1} \\leftarrow \\frac{F_{1}}{2} - f_{2} + f_{6} - f_{7} - f_{8} + vel1Term$$</td>  </tr> <tr style=\"border:none\"> <td style=\"border:none\">$$forceTerm_{0} \\leftarrow \\left(1 - \\frac{\\omega}{2}\\right) \\left(- \\frac{4 F_{0} u_{0}}{3} - \\frac{4 F_{1} u_{1}}{3}\\right)$$</td>  </tr> <tr style=\"border:none\"> <td style=\"border:none\">$$forceTerm_{1} \\leftarrow \\left(1 - \\frac{\\omega}{2}\\right) \\left(- \\frac{F_{0} u_{0}}{3} + \\frac{F_{1} \\left(2 u_{1} + 1\\right)}{3}\\right)$$</td>  </tr> <tr style=\"border:none\"> <td style=\"border:none\">$$forceTerm_{2} \\leftarrow \\left(1 - \\frac{\\omega}{2}\\right) \\left(- \\frac{F_{0} u_{0}}{3} + \\frac{F_{1} \\left(2 u_{1} - 1\\right)}{3}\\right)$$</td>  </tr> <tr style=\"border:none\"> <td style=\"border:none\">$$forceTerm_{3} \\leftarrow \\left(1 - \\frac{\\omega}{2}\\right) \\left(\\frac{F_{0} \\left(2 u_{0} - 1\\right)}{3} - \\frac{F_{1} u_{1}}{3}\\right)$$</td>  </tr> <tr style=\"border:none\"> <td style=\"border:none\">$$forceTerm_{4} \\leftarrow \\left(1 - \\frac{\\omega}{2}\\right) \\left(\\frac{F_{0} \\left(2 u_{0} + 1\\right)}{3} - \\frac{F_{1} u_{1}}{3}\\right)$$</td>  </tr> <tr style=\"border:none\"> <td style=\"border:none\">$$forceTerm_{5} \\leftarrow \\left(1 - \\frac{\\omega}{2}\\right) \\left(\\frac{F_{0} \\left(2 u_{0} - 3 u_{1} - 1\\right)}{12} + \\frac{F_{1} \\left(- 3 u_{0} + 2 u_{1} + 1\\right)}{12}\\right)$$</td>  </tr> <tr style=\"border:none\"> <td style=\"border:none\">$$forceTerm_{6} \\leftarrow \\left(1 - \\frac{\\omega}{2}\\right) \\left(\\frac{F_{0} \\left(2 u_{0} + 3 u_{1} + 1\\right)}{12} + \\frac{F_{1} \\left(3 u_{0} + 2 u_{1} + 1\\right)}{12}\\right)$$</td>  </tr> <tr style=\"border:none\"> <td style=\"border:none\">$$forceTerm_{7} \\leftarrow \\left(1 - \\frac{\\omega}{2}\\right) \\left(\\frac{F_{0} \\left(2 u_{0} + 3 u_{1} - 1\\right)}{12} + \\frac{F_{1} \\left(3 u_{0} + 2 u_{1} - 1\\right)}{12}\\right)$$</td>  </tr> <tr style=\"border:none\"> <td style=\"border:none\">$$forceTerm_{8} \\leftarrow \\left(1 - \\frac{\\omega}{2}\\right) \\left(\\frac{F_{0} \\left(2 u_{0} - 3 u_{1} + 1\\right)}{12} + \\frac{F_{1} \\left(- 3 u_{0} + 2 u_{1} - 1\\right)}{12}\\right)$$</td>  </tr> </table><div>Main Assignments:</div><table style=\"border:none; width: 100%; \"><tr style=\"border:none\"> <td style=\"border:none\">$$d_{0} \\leftarrow f_{0} + forceTerm_{0} + \\omega \\left(- f_{5} - f_{6} - f_{7} - f_{8} + \\frac{\\rho u_{0}^{2}}{3} + \\frac{\\rho u_{1}^{2}}{3} + \\frac{\\rho}{9}\\right) - \\omega \\left(- f_{1} - f_{2} - f_{5} - f_{6} - f_{7} - f_{8} + \\rho u_{1}^{2} + \\frac{\\rho}{3}\\right) - \\omega \\left(- f_{3} - f_{4} - f_{5} - f_{6} - f_{7} - f_{8} + \\rho u_{0}^{2} + \\frac{\\rho}{3}\\right) + \\omega \\left(- f_{0} - f_{1} - f_{2} - f_{3} - f_{4} - f_{5} - f_{6} - f_{7} - f_{8} + \\rho\\right)$$</td>  </tr> <tr style=\"border:none\"> <td style=\"border:none\">$$d_{1} \\leftarrow f_{1} + forceTerm_{1} - \\frac{\\omega \\left(- f_{5} - f_{6} + f_{7} + f_{8} + \\frac{\\rho u_{1}}{3}\\right)}{2} + \\frac{\\omega \\left(- f_{1} + f_{2} - f_{5} - f_{6} + f_{7} + f_{8} + \\rho u_{1}\\right)}{2} - \\frac{\\omega \\left(- f_{5} - f_{6} - f_{7} - f_{8} + \\frac{\\rho u_{0}^{2}}{3} + \\frac{\\rho u_{1}^{2}}{3} + \\frac{\\rho}{9}\\right)}{2} + \\frac{\\omega \\left(- f_{1} - f_{2} - f_{5} - f_{6} - f_{7} - f_{8} + \\rho u_{1}^{2} + \\frac{\\rho}{3}\\right)}{2}$$</td>  </tr> <tr style=\"border:none\"> <td style=\"border:none\">$$d_{2} \\leftarrow f_{2} + forceTerm_{2} + \\frac{\\omega \\left(- f_{5} - f_{6} + f_{7} + f_{8} + \\frac{\\rho u_{1}}{3}\\right)}{2} - \\frac{\\omega \\left(- f_{1} + f_{2} - f_{5} - f_{6} + f_{7} + f_{8} + \\rho u_{1}\\right)}{2} - \\frac{\\omega \\left(- f_{5} - f_{6} - f_{7} - f_{8} + \\frac{\\rho u_{0}^{2}}{3} + \\frac{\\rho u_{1}^{2}}{3} + \\frac{\\rho}{9}\\right)}{2} + \\frac{\\omega \\left(- f_{1} - f_{2} - f_{5} - f_{6} - f_{7} - f_{8} + \\rho u_{1}^{2} + \\frac{\\rho}{3}\\right)}{2}$$</td>  </tr> <tr style=\"border:none\"> <td style=\"border:none\">$$d_{3} \\leftarrow f_{3} + forceTerm_{3} + \\frac{\\omega \\left(f_{5} - f_{6} + f_{7} - f_{8} + \\frac{\\rho u_{0}}{3}\\right)}{2} - \\frac{\\omega \\left(f_{3} - f_{4} + f_{5} - f_{6} + f_{7} - f_{8} + \\rho u_{0}\\right)}{2} - \\frac{\\omega \\left(- f_{5} - f_{6} - f_{7} - f_{8} + \\frac{\\rho u_{0}^{2}}{3} + \\frac{\\rho u_{1}^{2}}{3} + \\frac{\\rho}{9}\\right)}{2} + \\frac{\\omega \\left(- f_{3} - f_{4} - f_{5} - f_{6} - f_{7} - f_{8} + \\rho u_{0}^{2} + \\frac{\\rho}{3}\\right)}{2}$$</td>  </tr> <tr style=\"border:none\"> <td style=\"border:none\">$$d_{4} \\leftarrow f_{4} + forceTerm_{4} - \\frac{\\omega \\left(f_{5} - f_{6} + f_{7} - f_{8} + \\frac{\\rho u_{0}}{3}\\right)}{2} + \\frac{\\omega \\left(f_{3} - f_{4} + f_{5} - f_{6} + f_{7} - f_{8} + \\rho u_{0}\\right)}{2} - \\frac{\\omega \\left(- f_{5} - f_{6} - f_{7} - f_{8} + \\frac{\\rho u_{0}^{2}}{3} + \\frac{\\rho u_{1}^{2}}{3} + \\frac{\\rho}{9}\\right)}{2} + \\frac{\\omega \\left(- f_{3} - f_{4} - f_{5} - f_{6} - f_{7} - f_{8} + \\rho u_{0}^{2} + \\frac{\\rho}{3}\\right)}{2}$$</td>  </tr> <tr style=\"border:none\"> <td style=\"border:none\">$$d_{5} \\leftarrow f_{5} + forceTerm_{5} + \\frac{\\omega \\left(- f_{5} - f_{6} + f_{7} + f_{8} + \\frac{\\rho u_{1}}{3}\\right)}{4} - \\frac{\\omega \\left(f_{5} - f_{6} - f_{7} + f_{8} + \\rho u_{0} u_{1}\\right)}{4} - \\frac{\\omega \\left(f_{5} - f_{6} + f_{7} - f_{8} + \\frac{\\rho u_{0}}{3}\\right)}{4} + \\frac{\\omega \\left(- f_{5} - f_{6} - f_{7} - f_{8} + \\frac{\\rho u_{0}^{2}}{3} + \\frac{\\rho u_{1}^{2}}{3} + \\frac{\\rho}{9}\\right)}{4}$$</td>  </tr> <tr style=\"border:none\"> <td style=\"border:none\">$$d_{6} \\leftarrow f_{6} + forceTerm_{6} + \\frac{\\omega \\left(- f_{5} - f_{6} + f_{7} + f_{8} + \\frac{\\rho u_{1}}{3}\\right)}{4} + \\frac{\\omega \\left(f_{5} - f_{6} - f_{7} + f_{8} + \\rho u_{0} u_{1}\\right)}{4} + \\frac{\\omega \\left(f_{5} - f_{6} + f_{7} - f_{8} + \\frac{\\rho u_{0}}{3}\\right)}{4} + \\frac{\\omega \\left(- f_{5} - f_{6} - f_{7} - f_{8} + \\frac{\\rho u_{0}^{2}}{3} + \\frac{\\rho u_{1}^{2}}{3} + \\frac{\\rho}{9}\\right)}{4}$$</td>  </tr> <tr style=\"border:none\"> <td style=\"border:none\">$$d_{7} \\leftarrow f_{7} + forceTerm_{7} - \\frac{\\omega \\left(- f_{5} - f_{6} + f_{7} + f_{8} + \\frac{\\rho u_{1}}{3}\\right)}{4} + \\frac{\\omega \\left(f_{5} - f_{6} - f_{7} + f_{8} + \\rho u_{0} u_{1}\\right)}{4} - \\frac{\\omega \\left(f_{5} - f_{6} + f_{7} - f_{8} + \\frac{\\rho u_{0}}{3}\\right)}{4} + \\frac{\\omega \\left(- f_{5} - f_{6} - f_{7} - f_{8} + \\frac{\\rho u_{0}^{2}}{3} + \\frac{\\rho u_{1}^{2}}{3} + \\frac{\\rho}{9}\\right)}{4}$$</td>  </tr> <tr style=\"border:none\"> <td style=\"border:none\">$$d_{8} \\leftarrow f_{8} + forceTerm_{8} - \\frac{\\omega \\left(- f_{5} - f_{6} + f_{7} + f_{8} + \\frac{\\rho u_{1}}{3}\\right)}{4} - \\frac{\\omega \\left(f_{5} - f_{6} - f_{7} + f_{8} + \\rho u_{0} u_{1}\\right)}{4} + \\frac{\\omega \\left(f_{5} - f_{6} + f_{7} - f_{8} + \\frac{\\rho u_{0}}{3}\\right)}{4} + \\frac{\\omega \\left(- f_{5} - f_{6} - f_{7} - f_{8} + \\frac{\\rho u_{0}^{2}}{3} + \\frac{\\rho u_{1}^{2}}{3} + \\frac{\\rho}{9}\\right)}{4}$$</td>  </tr> </table>"
284
285
      ],
      "text/plain": [
Markus Holzer's avatar
Markus Holzer committed
286
       "AssignmentCollection: d_1, d_2, d_7, d_4, d_3, d_0, d_8, d_5, d_6 <- f(omega, F_1, F_0, f_5, f_4, f_0, f_7, f_3, f_2, f_1, f_8, f_6)"
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "collision_rule = method.get_collision_rule()\n",
    "collision_rule"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Generic simplification strategy - common subexpresssion elimination"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
Markus Holzer's avatar
Markus Holzer committed
314
       "<table style=\"border:none\"><tr><th>Name</th><th>Runtime</th><th>Adds</th><th>Muls</th><th>Divs</th><th>Total</th></tr><tr><td>OriginalTerm</td><td>-</td> <td>293</td> <td>261</td> <td>0</td>  <td>554</td> </tr><tr><td>sympy_cse</td><td>45.13 ms</td> <td>114</td> <td>67</td> <td>0</td>  <td>181</td> </tr></table>"
315
316
      ],
      "text/plain": [
Markus Holzer's avatar
Markus Holzer committed
317
       "<pystencils.simp.simplificationstrategy.SimplificationStrategy.create_simplification_report.<locals>.Report at 0x7f3aca401190>"
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "generic_strategy = ps.simp.SimplificationStrategy()\n",
    "generic_strategy.add(ps.simp.sympy_cse)\n",
    "generic_strategy.create_simplification_report(collision_rule)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### A custom simplification strategy for moment-based methods"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
Martin Bauer's avatar
Martin Bauer committed
342
   "outputs": [],
343
   "source": [
Martin Bauer's avatar
Martin Bauer committed
344
345
346
    "simplification_strategy = create_simplification_strategy(method)\n",
    "simplification_strategy.create_simplification_report(collision_rule)\n",
    "simplification_strategy.add(ps.simp.sympy_cse)"
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Seeing the simplification in action"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
Martin Bauer's avatar
Martin Bauer committed
364
       "<h5 style=\"padding-bottom:10px\">Initial Version</h5> <div style=\"padding-left:20px;\">$$d_{7} \\leftarrow f_{7} + forceTerm_{7} - \\frac{\\omega \\left(- f_{5} - f_{6} + f_{7} + f_{8} + \\frac{\\rho u_{1}}{3}\\right)}{4} + \\frac{\\omega \\left(f_{5} - f_{6} - f_{7} + f_{8} + \\rho u_{0} u_{1}\\right)}{4} - \\frac{\\omega \\left(f_{5} - f_{6} + f_{7} - f_{8} + \\frac{\\rho u_{0}}{3}\\right)}{4} + \\frac{\\omega \\left(- f_{5} - f_{6} - f_{7} - f_{8} + \\frac{\\rho u_{0}^{2}}{3} + \\frac{\\rho u_{1}^{2}}{3} + \\frac{\\rho}{9}\\right)}{4}$$</div><h5 style=\"padding-bottom:10px\">expand</h5> <div style=\"padding-left:20px;\">$$d_{7} \\leftarrow - f_{7} \\omega + f_{7} + forceTerm_{7} + \\frac{\\omega \\rho u_{0}^{2}}{12} + \\frac{\\omega \\rho u_{0} u_{1}}{4} - \\frac{\\omega \\rho u_{0}}{12} + \\frac{\\omega \\rho u_{1}^{2}}{12} - \\frac{\\omega \\rho u_{1}}{12} + \\frac{\\omega \\rho}{36}$$</div><h5 style=\"padding-bottom:10px\">replace_second_order_velocity_products</h5> <div style=\"padding-left:20px;\">$$d_{7} \\leftarrow - f_{7} \\omega + f_{7} + forceTerm_{7} + \\frac{\\omega \\rho u_{0}^{2}}{12} - \\frac{\\omega \\rho u_{0}}{12} + \\frac{\\omega \\rho u_{1}^{2}}{12} - \\frac{\\omega \\rho u_{1}}{12} + \\frac{\\omega \\rho \\left(u0Pu1^{2} - u_{0}^{2} - u_{1}^{2}\\right)}{8} + \\frac{\\omega \\rho}{36}$$</div><h5 style=\"padding-bottom:10px\">expand</h5> <div style=\"padding-left:20px;\">$$d_{7} \\leftarrow - f_{7} \\omega + f_{7} + forceTerm_{7} + \\frac{\\omega \\rho u0Pu1^{2}}{8} - \\frac{\\omega \\rho u_{0}^{2}}{24} - \\frac{\\omega \\rho u_{0}}{12} - \\frac{\\omega \\rho u_{1}^{2}}{24} - \\frac{\\omega \\rho u_{1}}{12} + \\frac{\\omega \\rho}{36}$$</div><h5 style=\"padding-bottom:10px\">factor_relaxation_rates</h5> <div style=\"padding-left:20px;\">$$d_{7} \\leftarrow f_{7} + forceTerm_{7} + \\omega \\left(- f_{7} + \\frac{\\rho u0Pu1^{2}}{8} - \\frac{\\rho u_{0}^{2}}{24} - \\frac{\\rho u_{0}}{12} - \\frac{\\rho u_{1}^{2}}{24} - \\frac{\\rho u_{1}}{12} + \\frac{\\rho}{36}\\right)$$</div><h5 style=\"padding-bottom:10px\">replace_density_and_velocity</h5> <div style=\"padding-left:20px;\">$$d_{7} \\leftarrow f_{7} + forceTerm_{7} + \\omega \\left(- f_{7} + \\frac{\\rho u0Pu1^{2}}{8} - \\frac{\\rho u_{0}^{2}}{24} - \\frac{\\rho u_{0}}{12} - \\frac{\\rho u_{1}^{2}}{24} - \\frac{\\rho u_{1}}{12} + \\frac{\\rho}{36}\\right)$$</div><h5 style=\"padding-bottom:10px\">replace_common_quadratic_and_constant_term</h5> <div style=\"padding-left:20px;\">$$d_{7} \\leftarrow f_{7} + forceTerm_{7} + \\omega \\left(- f_{7} + \\frac{f_{eq common}}{36} + \\frac{\\rho u0Pu1^{2}}{8} - \\frac{\\rho u_{0}}{12} - \\frac{\\rho u_{1}}{12}\\right)$$</div><h5 style=\"padding-bottom:10px\">factor_density_after_factoring_relaxation_times</h5> <div style=\"padding-left:20px;\">$$d_{7} \\leftarrow f_{7} + forceTerm_{7} + \\omega \\left(- f_{7} + \\frac{f_{eq common}}{36} + \\rho \\left(\\frac{u0Pu1^{2}}{8} - \\frac{u_{0}}{12} - \\frac{u_{1}}{12}\\right)\\right)$$</div><h5 style=\"padding-bottom:10px\">subexpression_substitution_in_main_assignments</h5> <div style=\"padding-left:20px;\">$$d_{7} \\leftarrow f_{7} + forceTerm_{7} + \\omega \\left(- f_{7} + \\frac{f_{eq common}}{36} + \\rho \\left(\\frac{u0Pu1^{2}}{8} - \\frac{u0Pu1}{12}\\right)\\right)$$</div><h5 style=\"padding-bottom:10px\">add_subexpressions_for_divisions</h5> <div style=\"padding-left:20px;\">$$d_{7} \\leftarrow f_{7} + forceTerm_{7} + \\omega \\left(- f_{7} + \\frac{f_{eq common}}{36} + \\rho \\left(\\frac{u0Pu1^{2}}{8} - \\frac{u0Pu1}{12}\\right)\\right)$$</div><h5 style=\"padding-bottom:10px\">sympy_cse</h5> <div style=\"padding-left:20px;\">$$d_{7} \\leftarrow f_{7} + forceTerm_{7} + \\omega \\left(\\rho \\left(- \\xi_{95} + \\xi_{96}\\right) + \\xi_{64} + \\xi_{92}\\right)$$</div>"
365
366
      ],
      "text/plain": [
Markus Holzer's avatar
Markus Holzer committed
367
       "<pystencils.simp.simplificationstrategy.SimplificationStrategy.show_intermediate_results.<locals>.IntermediateResults at 0x7f3aca39b160>"
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "simplification_strategy.show_intermediate_results(collision_rule, symbols=[sp.Symbol(\"d_7\")])"
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "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",
Frederik Hennig's avatar
Frederik Hennig committed
397
   "version": "3.8.2"
398
399
400
401
402
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}