demo_create_method_from_scratch.ipynb 67.7 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": "iVBORw0KGgoAAAANSUhEUgAAAoIAAAAUCAYAAAAZUFxiAAAACXBIWXMAAA7EAAAOxAGVKw4bAAAIhklEQVR4Ae2dgXXcKBCGvang3l0HTgfOXQe+DpKXChx3kLwrIenA6eBergOnhbiDpIMk7sD3/zI4iIUBhIRggfdYSTBoho8ZCWml3cPZ2dlb5NfITFcPDw93j6u/Pg+Hwzm2LlH38VdpPWuw7yWsuYN93/awqnY+ZDIYhT1jMJIZDT4yH9YORoNRmIAsMXxI5sPawUhmFOKD+jfYwzXyb1we8PEe+RaTqM9YHiU0oOAn1P9tV6KObb8j3yM/R/4XckcTSZRnpRg9kPkEJe9KTwZr4BMLt2dGapzor0y8sGGiv8z8tUZGj6Zu/4m+B+O5Zz6xI9AzI/Sd54vm4izW7lgfkOSgq7k4K8lHYmfWwaaqzvmlGMXqieEDGfriLblyhXf7zlwZdV+Qz+06Vf5Sl2ObB4CvLlkts2QZq0fp/7JER04bZd9ufEzbYQvvjP4kC7Ncr7McuTtGqt+8mHnycZTR7x+QZ76vZKthpG2GXeLYarmlS+yfcR6M5175mFxDY9ErI9Xv5uIsxW7TD5asQ1dzcVaSj8kUesVjnrKrimN1KUYpemL4QGaa/2HpnwiijrcPb8zB4boq/+ko505XG5hUPZCn/re2XVttK/t246PGgpM7Xhmx71xygJwTQSXfI6MbFxOUcdLs8+Ma/ChpbJf6ORgwzn0cjuIZsrX4UBE+jcdZSUatxlmS3R3GWRE+jcdZEUY49ibpCR2rVf0lluJE0HmHD214K3F25acGkTN4cSKSEkSpeiDPg97RCS1FZ4osdO3Kx7YV9vCELvJHfXeM0GdO+L46eOmJ8+yObi2MTHthU3BsTfmUdew7KZ5742OzjBmLHhmhz03GWardtj/EbkNPk3FWio/NEXqDxzzI1HI+S/J9u6+x26ljEeKDel7UXz7DhzPhu+PpOSoY6HoBY5pBOhreqzJfvaOJWOTbj1MPbGX5N9juaycqS6mshE+KyZNsp4zow+cYMx40XGlWXhEjl61blPnipfY424LFKvvs0IfIrdU4S7I7w0FajbNSfJLRVhRnpRgl6Ynl450IYkR4d+/oBRLjZPpDGDX9ML4gIldl6KHNr+S9r1K7K5/MHnTFCMHwAvmggsJEd6E2GFx22pWRbcxW263G2VY8Vt5vFz6kmbUaZwvt1t2OWrYcZyX4REH0C+0eZ6UYLdQT5CNNBPmWMB9stdPvqkDfLTDrpcmhKRezvlQPv679M0ZBpszefHLM754RDsy8OucFywcEl8uX92aUM74pbVuNs5Q+7iXbiw95+bYaZxF2e/vsqTipONuAjwdbVHGVcVaKUYSeIB9pIsiTpDSxm32dZg3XH9Z2zmaqHtostcmxxWxbCx/Tptj1wejxodv/MAl854FWCyOPeasXSzHjiufe+CwBPhi1G2d8KF86PizxB7Y5lTjbis8SrrXGWSlGIT1BPtJEkFcw945R4U59SV/1fPcJJJQv1TN9h56gZ6no3nyW2s12XTPCFRQD5zMmgdIjBHszyhnflLatxllKH/eS7cWHnHxbjbNIu519FgpPJs424iOgC1ZVF2elGEXqCfKRJoJO+sbXaHrSZ8rpqx0qzkql9GQZ6Wjcqt2OrmxWtCcjBA7/SYe/J8hfVe8+7TkW3cM/YQCtxtlWdp9KnG3F55RCoRSjNfVIE0FeweiJnT1OfPjQVacnh6xfIy3RQxuyJ6IRxtfAJ8JMp0iXjBA4fMHnuTkJRBnfJubX/HaqgZFt01bbrcbZVjzW2m9PPvTErNU4S7T7qb8JK03HWQE+CShnotXEWSlGiXqCfKSJICdT3IEr8ffXXC9k8G/m+J+/rq+UXfsJlS3Rw5P6Wvol+2rgI9kn1XXHCIFzASB/mZNABYiTQ1eqgZHLri3KWo2zLVisuc+efGji1mqcLbB7iZ80G2eF+CxhyjZVxFkpRgv0BPk8E8jfoe6Fqx4n048o/wGDnk6iWOcdQv4A5JVuwzLkB2TX28dazLuM1WPtgDaLdyRz7VL6dudj9Zubrru0DrFpXLthhPFmIPAgzD8rf6/yDZYsu4afue4g7+5HtNdI3rFFP7qMM4MNVzfjY+kRdVmy3fgQ+w0/XD3Ocn3bGA/v8TrW7lxbWj2fleJjjJVe9ca0FlDL3eOsFKNYPal8KD/9sjSc9Ol/WLmOxJ/XOPo3Bi2HOg7S9FdTWPK5K55UL3S9XqKMry4ze//2TMu6lmyHHNSj20KWk84jO3S9XkIm164q+LA/SHz5gfz56+YsIAOWvdH9NZeqvhtGqr/k4spOH6+IUdTYwt5cf241zorwaTzOijDaKs5yfVuNnfd4nWJ3ri1o31ycleTTapyVYpSiR5/zVRvn+R510/wPS/dEUA0ITy6zv9/SO09ZYh+8c7hoIpioh0HmPLG79pNrF3UhN8OHDJAGI+uix/aNwWh+UTj4yDxsPiPO4nilxBlks88h2EdTx+sUPsrnshi1xmfEWTjOQj6E+mkiKH01DJlJyPc7a6yPTXw26z5WOEPuH7Tl1W9syrWLEFviQy6DUdg7BiOZ0eAj82HtYLQuo9xjNa1p7Xg9fGhdH+Lecv3oJH1InAhi8sZnAX1vVYaHCBJ8rgKLNX5XUNSn9FzC5g+ioKpcw66W+LDbg1HYMwYjmdHgI/Nh7WC0LiPFM/sc0tLxevjQuj5kxGWWH52yD023BtHB2TOCehsAOZG71dupS7R1PqeWup+QPG1Ejv6adi27sJ8m+JAf0mDk8XPtX4OR+zgw+MhcNJ8RZ3GcUuIMsqudQ7CvJo7XKXyUz63CqBU+I87CcRbjQ5CZ5n8HrPBFD/647j3yFQDz7apZwtUJ3wbj3TbeIawuwT6+rcx/inC9/bm5vbXzIYDBKOwGg5HMaPCR+bB2MBqMwgRkieFDMh/WDkYyoxAf1HPe9xqZF0bX/wMgEMnxFJLSjwAAAABJRU5ErkJggg==\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": "iVBORw0KGgoAAAANSUhEUgAAAVoAAAAcCAYAAAAgJ3ezAAAACXBIWXMAAA7EAAAOxAGVKw4bAAAHqUlEQVR4Ae1dgXHcNhD816QAKakgdgd2XEGkDqykA6cDeVxBxulATgWJ3YGdChKrA6sDW9+BsgvxOCCfBA4kiL8XgRk8CeB4t7eAjngQfG3v7+83NYUZ2G63byFxivwE+Rb5NXjb4VjTkTFgtS8t4rKIKWW45cY/R9+2Btpw1zXkXoMnBtgNyu9xeILy8/CVtdUaA1b70iIui5hSxlNu/HP1nfTBQ+Ez5PN+/YrLL3u+/44yOeLstqbjYsBqX1rEZRFTymjLjV+lD3HhaghkJ9AyyELoDWZrn4aEV1znB1VZMvDrVkzN0bnu95ulvrSIyyKmlAGXG79G3w3iKL/1dtJ3ndJm8yfKP/fqVl3ETedpjwDejJj+ezjUz2NhwGpfWsRlEVPKOMuNX6sPcp8QaC+QX+H8nWBu12jRcI3KL2j8QxrrcZ8B8PQFtVyzrTzt03NUNVb70iIui5hSBltu/DF9TftzxAn3rckFWlRySswAciYNKU6sRRY8ud0H4Oi3tfj8WP202pcWcVnElDIuc+PX6IMM12qfSqw4aQAzcLyLBVlc/BL5Dvk0xdHHIAufX8GPUyHuMfi0Vh+s9qVFXBYxpYzb3PgT9H0ATsYMlyTQsuJjU9c5QPEp8ntkzuZ+RV5jkOUujPbuBC7OkWWttsNXLdhmgH1nsS8t4rKIKWV05cafog8TMm4HvcU1brfCFgX1sgEuYkDmWu5qlhjgMwMqHxJyW5ckfgO4jH0DEOF6tMGA1b60iMsippRRlBv/FH24hrGSe+4vGGgZPN+icBZzBBeuMdDegZe9WTz4Inc1HREDGL8m+9IiLouYUoZabvxT9OEartNyu+wZt3fxDSf31lOKI7llAYrBjIGc26k6ux/Qxln3Z+QfAXqH46yUYktzA5oFprnYw8SaF8ivkek3Z9Q/IH8Flmw7HVLsebKL9w385Nt3MhZYDHLhyUaxlerLVB9K4fK4ivKqweTpi3JPg3OTZ4+qguMiN36NvgH/blDHpVc3UePa7Eco2sQy5BgI+eMIfCgUlU+RgU7OqnHYcA3tzr8WZd4ZUJXHJnQVs6XBDDzsCIeJ8kjkmbtArpoyzz9rdGlkoCvJHuSL8WUZm4bbpr+S+NXqnSOXyqvGluVxYQE/+OEkiX/Qzzij/R75oDNaRHzO3P5FZrpE7uO5QF2Wt9VK2qIzyvQGcv4aMGcI5EQ2PHNNuM8JqiYntb0D8GUZm5ZwtQ9ahRnksmKyPC40XBXCv2uwMMa6mRM34EdnqJBdZEYLve0MGedcR3MzOcGEMu8KnTppSz1CTzFbWmw+Jl6DpP6WobXhy6XY82VxvmjfNL63/RPjojQ2n8PQuY8r5kNIT8623Jh8fTg3NS40vJXATxvI/IM+P8GHRF2cHiaBGIcBdxkuGxCczOS4XsfpN1OWGW1JWw+w45+CyZMkD4Pb7TyZyacp9kS2RN/QIbHnOTfKhciWwuZhCp4KLk9o1AdPZtHT3JhEXynuxZ5H0ixORd/C+N1MlpgZaL8htxWsPGDissGNkNDgIKE71HFhOWcqaUuNGx2f9cYSM5xgrzhflrHFeJX2BB/kksWPmTFZHhcaLpfEz0kj063MaKXiofpwn1yX7K9FZluf7blV0lbPdLeIgc837qQPeGPhzK69saDtymvvXjyhNNFeEb4sY9NSPdEHrfpJcgtisjwuNFwtiV8msN9OgIQPoX7SIIKMBIM9cXQktzHcI3Mb1tTUCbLQxaCz9xWhpC2NI3PwND7yZ9VcgMWRN5Y2oZ0DgTP6nVSWttfYXbxvpnCRgk34Cx3ncEu9M3wIwXJb3qB70t/XUpgawJbHRZDTFPwaRQMy7d8uN91LIBt92wsddQ05RmfKMthytsWfCeSWI389lduQmNpfrXko6j5hh7r5FpbsQODTd7e3Fnb6HVrMlgY9sE/CM+Az+XQ+4+huWj7HguUA9hbvmxlcqLEJf6HjVG6pc6oPITzSNhXXwpjU3B87fukH7RH+Mm66N8PcTgMU+NTwHH/Q0Z0HMRno4bu9nSfHsWvG2qGHQEf3j5a0NYbRr8+Jx9c7dl7ano8Dtov1jW9Xcx7DptSRbRxr7GllDtnnGowx7o8dv4YDkYGvnCi53VJcOmD6G7nzldXVTvt4AUO71EsR/fnDNe2yA855p+TMjm9IjaWStsYw+PWT8PgKEs+L2CvZN4n+cwY5ZdxozBThVgOkJ2MG10Tujx1/rzuGi0384oNt941fAi1nJwxqs1Kj/OtEJQRFHJK4bsn/Nju4raukLQEUOs7EE1I92FbYXsm+GfQ3UJmELaCnbSrMbWs3dmIQVxL3x44/1j+99l9Q/iCTTv8/LDCw/YUG/o7ipAQiO/++IUUJruX6LzuOiWuzfImiffLuar2PkrY8s6Onc/CMKg00lLRXsm8CLg82pWIbVNKrLMltz3SwaA1XKvfHjj/YOb1G+MpnNpcSw/xAy6/q/6Ch/hvtHmm1WBmoDFQGtAzwhgJZ/n51u+zZBloqgQAfAFxAgO/W11QZqAxUBioDCQwghnJLF7+Nd555dQIt9TXBlvs2B9dGE2xW0cpAZaAysCoGED/5nInPljobAvYC7apYqc5WBioDlYECDPwPXqbwhQzG1qsAAAAASUVORK5CYII=\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": "iVBORw0KGgoAAAANSUhEUgAAAYkAAAF5CAYAAABqT9akAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3Xl0XmWB+PHvk3RNW8qSUspSoQUXBAb0yCIiKG4BfwPiyDAyqIgL6igj6oz83BiP4CAyqHUZF9zBUWdk9PAjIgdQUFxREXcQKhRa2hRaupckz++PN2+bpHmTd73r93MOp+TNvW+eP27eb+597hJijEiSNJGutAcgScouIyFJqmlaUj8o9A90Ay8ATgDmAaHGohFYB3w/9vXelNDwpJpC/8BM4FTgWGA2brvKidA/0AO8GDgamDXJohFYC9wc+3pvHfMeSc1JhP6By4DTG1zti7Gv94OdGI9Uj5E/bq4Cjmtw1S/Fvt5LOzAkqS4jf9x8GTiywVWXxb7ej1e/SORwU+gfWEzjgQA4J/QP7NHu8UgNOJrGAwHwj6F/YM92D0ZqwIk0HgiA14b+gdnVL5Kak/ibJtfrBg5v50CkBrWy7R7WzoFIDWp2250FPLn6RVKRmJnSulKr3HaVVzNaWHfHtpvYxPWELnzRAdxz1yy6uytf77H3IJ//xX2pjkmayjeX7c4t35zPintmcNwpG7joc6vSHpLUkO9dM4+vX7kXa1dNZ/5eg1zwkVU87aQtEy2abiQAznvfak573fq0hyHVba9Fg5x5wVruuGUO27fWOtNJyqaffLeHr3xwAe/4z4c47NitrHlo0g6kHwkpb577dxsB+POvZ7F2pb9DypdrLu/lZW9ZyxHHbwVg4QGDky2e/sV0V1/ey5kHL+Utz1vML26aPfUKkqSmDA3Cfb+fxfq13bzqqIM4+6lLuPIte7N1c8094nQjce571/CFO+7lq7+9lxeevY5LXr0/D9w9PdUxSVJRrV3VzdAg/Pj6eVx+3f18/Ja/ct/vZ/GlS/eqtUq6kTj8mVuZs1tkxqzIqec+xhOP3MJP+uekOiZJKqqZsytXT5967qMs2G+IPfYe4vTXP8Kvvl/zczf9w01jhIh3pZWkzpi/1zB77D1IqP98i/Qi8dgjXfz4+h62bQkMPg7f/co8/nRHD0c/f1NqY5LqMfg4bNsSGB6C4SF2bMNSHjznpeu57gu7s3ZVN+vXdvGdz+7B05+7sdbi6Z2ZMfh44CuX9fKhN8ykqyuy6KDtvPOzD3Lgof62Kdu+dMle/M8ndh7D/dF1u/HSN63lvIvXpjgqqT6vfPdaHnu0m9cddxDTZ0SO7dvAKy56pNbi6UViz4VDfPIH96f286VmnXexQVB+TZ8BFy5bzYXLVtezeMbmJCRJWZJUJFqZjXYmW2ly21UZ7dh2k4rE5pTWlVrVyvY34b1wpIS0ZdtNKhI/B4abWG8r8Ks2j0VqxE+bXG8bbrtKV7Pb7mPA76tfJBKJ2Nc7ACxrYtUPxb5e/xpTamJf713AN5tY9UOxr9e9YKXpJ0B/g+sMA5fEvt4d93NK7PGlAKF/4BDg2dT/jOvlCQ1NmlToHziCyjOue3DbVU6E/oEAHMXOZ1zX2naHgUepPON6xZj3SDISkqR88RRYqUkhhBtCCJemPQ6pk9yTkJoQQlgC/IHKyRULYozbUx6S1BHuSUjNecPIvwE4Lc2BSJ3knoTUoBDCdGANMH/kpZ/FGI9JcUhSx+R+TyKEsDCEcGMI4eUhhBlpj0el8LeM/d05IoRwcFqDUXmEEOaHEN4aQrh25I+Vjst9JKgcEz4BuAp4OITw/hDCwpTHpGJ7G5XTuKu62Xn4SWq7EMKTQwifA1YCHwIOByZ9NnW75D4SMcb1wJep3NF2d+AdwPIQwv+EEJ6R6uBUOCMT1keNe3k68Br3ZNVOIYSuEMKpIYQfAr8EXgnMpvKH8SUxobmC3EdixBXsrOqskf9OB74fQvidh6LURm9g4t8bJ7DVFtVDSsAK4L+A46nEYfSjHf4rsfEUZeJ6pLbH1/j2RioRWQZ8Isb4cGIDU2FMMGE9nhPYaloI4cnA24GXU7l6v2eCxbYDH4sxviOpcRVlTwLgUioxmMhcPBSl1o2fsB7PCWw1ZNQhpR9RuSFk9ZDSRIGAyu0zPpbU+KBYexJdVHbPFtWx+DCV43r3AcfFGDd0cmwqhhDC7cBxkyzyOLAsxvi2hIakHAshPBW4kcpJEHPrWCUC340xntLRgY1TmD2JGOMwcDn13UO9C5gJ7EGaj3BVbtSYsB7PCWw1IgC7UV8gADYBH+zccCZWmEiM+Dy173I42hDwMHBMjPHRzg5JBVFrwno8J7BVlxjjb4GTqX2YfLyHgR92bkQTK1QkRk6H/SqTnz88OhArJllOAnZMWL8WqGcPYR6VyUdpSjHGnwLPY+pQbCTB015HK1QkRlxB5djwRIapnJ1iINSI/0PticSJPCOEsLRTg1GxjApFrc8tqMxHJHba62iFOx4fY/xTCOGX7Ho67BCVK2P3oTJpLdXrN0z8dLqXj/x7zbjXt1PZW5XqFanMaU1kO/DpGGMqT+kszNlNo4UQTgG+zs4JoeohppOAP4+8tiDGOJD86FQUIYQI3BljPDLtsSi/QghHU3ke9Ubg+VTOeBo9mb0VeGKM8YEUhlfIw00A3wWqp7WOnoO4G5gz8vqaEEJvGoOTJNglELvFGH/C2DmKCNySViCgoJEYdTrsIOMmqWOMmzEUklI2QSAi7DKZvZ0UTnsdrZCRGPF54AtMMEltKCSlqVYgqkaF4hOkcNrraIWck6hXCKGHygUq4ByFGuSchJoxVSCypsh7ElNyj0JSkvIWCCh5JMBQSEpGHgMBRgIwFJI6K6+BACOxg6GQ1Al5DgQYiTEMhaR2ynsgwEjswlBIaociBAKMxIQMhaRWFCUQYCRqMhSSmlGkQICRmJShkNSIogUCjMSUDIWkehQxEGAk6mIoJE2mqIEAI1E3QyFpIkUOBBiJhhgKSaMVPRBgJBpmKCRBOQIBRqIphkIqt7IEAoxE0wyFVE5lCgQYiZYYCqlcyhYIMBItMxRSOZQxEGAk2sJQSMVW1kCAkWgbQyEVU5kDAUairQyFVCxlDwQYibYzFFIxGIgKI9EBhkLKNwOxk5HoEEMh5ZOBGMtIdJChkPLFQOzKSHSYoZDywUBMzEgkwFBI2WYgajMSCTEUUjYZiMkZiQQZCilbDMTUjETCDIWUDQaiPkYiBYZCSpeBqJ+RSImhkNJhIBpjJFJkKKRkGYjGGYmUGQopGQaiOUYiAwyF1FkGonlGIiMMhdQZBqI1RiJDDIXUXgaidUYiYwyF1B4Goj2MRAYZCqk1BqJ9jERGGQqpOQaivYxEhhkKqTEGov2MRMYZCqk+BqIzjEQOGAppcgaic4xEThgKaWIGorOMRI4YCmksA9F5RiJnDIVUYSCSYSRyyFCo7AxEcoxEThkKlZWBSJaRyDFDobIxEMkzEjlnKFQWBiIdRqIADIWKzkCkx0gUhKFQURmIdBmJAjEUKhoDkT4jUTCGQkVhILLBSBSQoVDeGYjsMBIFZSiUVwYiW4xEgRkK5Y2ByB4jUXCGQnlhILLJSJSAoVDWGYjsMhIlYSiUVQYi24xEiRgKZY2ByD4jUTKGQllhIPLBSJSQoVDaDER+GImSMhRKi4HIFyNRYoZCSTMQ+WMkSs5QKCkGIp+MhAyFOs5A5JeREGAo1DkGIt+MhHYwFGo3A5F/RkJjGAq1i4EoBiOhXRgKtcpAFIeR0IQMhZplIIrFSKgmQ6FGGYjiMRKalKFQvQxEMRkJTclQaCoGoriMhOpiKFSLgSg2I6G6GQqNZyCKz0ioIYZCVQaiHIyEGmYoZCDKw0ioKYaivAxEuRgJNc1QlI+BKB8joZYYivIwEOVkJNQyQ1F8BqK8jITawlAUl4EoNyOhtjEUxWMgZCTUVoaiOAyEwEioAwxF/hkIVRkJdYShyC8DodGMhDrGUOSPgdB4RkIdZSjyw0BoIkZCHWcoss9AqBYjoUQYiuwyEJqMkVBiDEX2GAhNxUgoUYYiOwyE6mEklDhDkT4DoXoZCaXCUKTHQKgRRkKpMRTJMxBqlJFQqgxFcgyEmmEklDpD0XkGQs0yEsoEQ9E5BkKtMBLKDEPRfgZCrTISyhRD0T4GQu1gJJQ5hqJ1BkLtYiSUSYaieQZC7WQklFmGonEGQu1mJJRphqJ+BkKdYCSUeYZiagZCnWIklAuGojYDoU4yEsoNQ7ErA6FOMxLKFUOxk4FQEoyEcsdQGAglx0gol8ocCgOhJBkJ5VYZQ2EglDQjoVwrUygMhNJgJJR7ZQiFgVBajIQKocihMBBKk5FQYRQxFAZCaTMSKpQihcJAKAuMhAqnCKEwEMoKI6FCynMoDISyxEiosPIYCgOhrDESKrQ8hcJAKIuMhAovD6EwEMoqI6FSyHIoDISyzEioNLIYCgOhrDMSKpUshcJAKA+MhEonC6EwEMoLI6FSSjMUBkJ5YiRUWmmEwkAob4yESi3JUBgI5ZGRUOklEQoDobwyEhKdDYWBUJ6FpLbX0D8wCzgDOB6YC4Qai0ZgPXAbcG3s6x1MZIASEELoATaNfLkgxjhQc9u96evnMGvOoxz/4utGlh8G1gI3xb7e60fez0AoNaF/YA7wUuA4oIfan7vDwBrge7Gv98Yx75HENhv6BwLwBSoDbcQNsa/3LR0YklTT+FBw/ZrLgWfusuCKe57ItBnb2GfxXyd4m6s4ZcH3MRBKSegf6AauAY5scNWPxr7eT1a/SOpw09/QeCAAXhj6B57Q7sFIkxlz6GnJYWsYGjyh4TfZuul1zJpjIJSm42g8EACvDv0D06tfJBWJw1pY94h2DSKEUGtXSxpjRygWPwlWLl/K0GB33Stv2zKLgZVL2X/pZgyEGtDmz6hmP3fnAQdWv0gqEjNSWpdQ8ZwQwneoTEguauX9VB4xxs286bLzgfpDsW3LLNY8uJgQIlfe8FwDoXqFEJ4FPBpC+EwI4fA2vGVbPnfTPbvpA6/ah3948lLOeMLBnPv0g/j2Z+a3661DCL0hhH8BVgDfAV4MzAZ2a9fPUAnMmb+dfZfcDYwNxfqBLj75r3D+8TM55/Al3PDVeWMCse+Su+melubIlT/7AdOBVwM/CSHcFUJ41cgcWfv99Y/T+dv9DuEDr9pnssXS3YrPuvARFj/xYWbMitz3uxlcdMYBHHLkVg49elszbzeyq3YS8FbgBVRm7GePWmSo5TGrfLq6Kh/6D917CCuXL2XRgX/hoxfuzbTp8JEbt/How6v5wCv3Z4+9A/strSzrkU01ZwjopnIm0mHAMuATIYSrgWUxxrva9pM+8Y6FLHnq1qkWS3dP4uAjtjNjVmV3PIRICPDgvQ3vItXYa5jJ2EBIzauGAmD5H5by8xvncdrrYFYPHHp05IgTAj/5LgZCbTaXSjDau3fxvWvm0bPbEIc/c/NUi6Z/Md1/vHlvTt//EN544kHsvmCQ40/dWM9q4+YaVgAXA/sy+TUYUvOqoXj4fujqgoWLIQ53sebBxRxwcGTVXzcZCHXI+L2LNU3PXWxc18XXrujl/EvX1LN4+gdNL1y2mguuXM1vfjSbO2+bzfSZU070jew1XEBlfmEO9UehCzgphHBQ8wNWqbz7i4dx1Ilj/2rrnjbA7LmVK7KHBqcTQmR+7zq2bOxhy8ady/73x48Np1yxR6LjVZ4dXedyc0f+fTVwdgjhXuCKGOMX61r7qn/r5blnrmefJ9R1oXL6kQDongZHnbiFm76xG9d+anfOvGDdFGtcRmW+odE9oTnAfzY1RpXTT2+AxU8a+9rQIGzZtPPrGAPr1+7BzNmwdtX+O17/5S0fSWiUKo5Gzoar7l08Bfg88MUp1/jTHTO56/YePnnr8np/SDYiUTU8CCuX1zMncTjwT8A5VGIxd/LFd9gAPCPG+KcmR6iSCf0D5wLvHPPi3N1nMzx0AA8/AAsPqLy28r4NPOEpj7P/wQM7lrvyhlfEvt6fJjhc5VgI4e+Bz1K5TqEeG4DtwCeAz9S1xq9u7WHgoem84oilAGzb0sXwMLzhhJl86raJ7hyQ4pzE2lXdfO+aeWzeEBgahB9f38Pt1+/GkSdMOZESY/xtjPF8YAHwJuAuYDPgfZ7UWdu2zGLjugN42knwnc/C0OPbuOdOuOOWeZx8Zl3zaVILtgFbgZuBs4CFMcb3xRgfrGvt0163js/99F4+fstyPn7Lcp739+s48oRNXPLfK2qtkt6eRAjQ/6Xd+fS7FhKHYa9Fg7zq3as58Yy6f9FGror9MvDlEMJhNLd3IdVn9HUQb/3YX7js/IO54AUzmTt/kLPfMY2eeYsZGryH7mnDaQ9VhTNmr6HuKIw3e05k9pydlwLMmjPM9JnD7Lmw5uUB6UViz4VDXHnDA+16uxjjb4HzQwgXAn8HvB1YSuXKwWwdVlP+jL9QLgR402XsuMHf8HAYuY7iYBYdaCjUDtuozFHcDlwB3BBjbO+1XuddvHaqRdI/BbbNYoybY4xfjjEeARwDXEXlUNRGKtdOSI2ZKBDjjb6OYuXygxkaLNzvlhIxk8pew1oqJ+gcHGM8OcZ4fdsDUaekNuRW7l/T9LoTzF18C1jZwlhUNr++9cApA1E1PhR//EW9E5ASwB1ULgZufK5hYm353E0qEq1M6LU8GThq7+IfYoyPtfp+KocQwtFc/6X31hWIqtGh+Mg/3xhC2KvDw1RBxBjviTG+rI17DZumXmTqdZOKxO1NrjdE5aEtUqJ2PFHujz/f3PCtNrq6IosO+hkr7gYYMBRKSbOfuytiX++O02ETiUTs630Q+FijqwEfjH296zswJKmmUY8c3cTAyrmE8NEG32KQ7u73MDxcfWa2oVDiYl/vH6lcZNeIrcD7Rr+Q2DOuAUL/wAFUHgM51bHa9cBtsa93VedHJe00JhAwr/o8iAm33Y++9bPMnvMAr/vA+0deiVQmHG+Nfb2PjLzf6Eeh9sYYpzybRGqn0D9wIJWn1M2ZZLEIrAZ+EPt6xxySTzQSUpbVCsQky0fgzhjjpI+INBTKM0/Tk2g8EI0Y88xsDz0pZ4yESq+TgagyFMorI6FSSyIQVYZCeWQkVFpJBqLKUChvjIRKKY1AVBkK5YmRUOmkGYgqQ6G8MBIqlSwEospQKA+MhEojS4GoMhTKOiOhUshiIKoMhbLMSKjwshyIKkOhrDISKrQ8BKLKUCiLjIQKK0+BqDIUyhojoULKYyCqDIWyxEiocPIciCpDoawwEiqUIgSiylAoC4yECqNIgagyFEqbkVAhFDEQVYZCaTISyr0iB6LKUCgtRkK5VoZAVBkKpcFIKLfKFIgqQ6GkGQnlUhkDUWUolCQjodwpcyCqDIWSYiSUKwZiJ0OhJBgJ5YaB2JWhUKcZCeWCgajNUKiTjIQyz0BMzVCoU4yEMs1A1M9QqBOMhDLLQDTOUKjdjIQyyUA0z1ConYyEMsdAtM5QqF2MhDLFQLSPoVA7GAllhoFoP0OhVhkJZYKB6BxDoVYYCaXOQHSeoVCzjIRSZSCSYyjUDCOh1BiI5BkKNcpIKBUGIj2GQo0wEkqcgUifoVC9jIQSZSCyw1CoHkZCiTEQ2WMoNBUjoUQYiOwyFJqMkVDHGYjsMxSqxUioowxEfhgKTcRIqGMMRP4YCo1nJNQRBiK/DIVGMxJqOwORf4ZCVUZCbWUgisNQCIyE2shAFI+hkJFQWxiI4jIU5WYk1DIDUXyGoryMhFpiIMrDUJSTkVDTDET5GIryMRJqioEoL0NRLkZCDTMQMhTlYSTUEAOhKkNRDkZCdTMQGs9QFJ+RUF0MhGoxFMVmJDQlA6GpGIriMhKalIFQvQxFMRkJ1WQg1ChDUTxGQhMyEGqWoSgWI6FdGAi1ylAUh5HQGAZC7WIoisFIaAcDoXYzFPlnJAQYCHWOocg3IyEDoY4zFPllJErOQCgphiKfjESJGQglzVDkj5EoKQOhtBiKfDESJWQglDZDkR9GomQMhLLCUOSDkSgRA6GsMRTZZyRKwkAoqwxFthmJEjAQyjpDkV1GouAMhPLCUGSTkSgwA6G8MRTZYyQKykAorwxFthiJAjIQyjtDkR1GomAMhIrCUGSDkSgQA6GiMRTpMxIFYSBUVIYiXUaiAAyEis5QpMdI5JyBUFkYinQYiRwzECobQ5E8I5FTBkJlZSiSZSRyyECo7AxFcoxEzhgIqcJQJMNI5IiBkMYyFJ1nJHLCQEgTMxSdZSRywEBIkzMUnWMkMs5ASPUxFJ1hJDLMQEiNMRTtZyQyykBIzTEU7WUkMshASK0xFO1jJDLGQEjtYSjaw0hkiIGQ2stQtM5IZISBkDrDULTGSGSAgZA6y1A0z0ikzEBIyTAUzTESKTIQUrIMReOMREoMhJQOQ9EYI5ECAyGly1DUz0gkzEBI2WAo6mMkEmQgpGwxFFMzEgkxEFI2GYrJGYkEGAgp2wxFbUaiwwyElA+GYmJGooMMhJQvhmJXRqJDDISUT4ZiLCPRAQZCyjdDsZORaDMDIRWDoagwEm1kIKRiMRRGom0MhFRMZQ+FkWgDAyEVW5lDYSRaZCCkcihrKIxECwyEVC5lDIWRaJKBkMqpbKEwEk0wEFK5lSkURqJBBkISlCcURqIBBkLSaGUIhZGok4GQNJGih8JI1MFASJpMkUNhJKZgICTVo6ihMBKTMBCSGlHEUBiJGgyEpGYULRRGYgIGQlIrihQKIzGOgZDUDkUJhZEYxUBIaqcihMJIjDAQkjoh76EwEhgISZ2V51CUPhIGQlIS8hqKUkfCQEhKUh5DUdhIhBCWhhAeCCG8sMb3DYSkxNUbihDC20IIfwghzExudLsqbCSAtwKLgG+ND4WBkJSmqUIRQngb8H5gMfCyhIc3Riji52MIYS7wMNAz8tJm4IwY4w0GQu0SQojAnTHGI9Mei/IphNBD5bMIoDfGuHZUIKqfX3+IMR6aygAp7p7EK4DRH/49VPYo/hUDISkjJtij+DfGBgJgcQjhmMQHN6JwexIhhAD8FThgksW6DITqFUJ4HvB1IIz71h4j/z467vXtwDExxr92emwqhnF7FOMNA9+JMb4kwSHtMC2NH9phJ7Pzl3cim4EXADckMxwVwGNU/rKbVeP747e37cD6jo5IRfMGKp9NPRN8rwt4UQhhUYxxZbLDKubhpv8LzJ3k+9VDTxOe9SRN4OfAQ3UuOwxcG2Nc18HxqEAmmIOYcDHgTcmMaKxCRSKEsAQ4ro5FDYXqNnJo8gpqHw4YbQvwkc6OSEVRZyAAZgL/FEKY0flRjVWoSFA57bW7zmV7gOtCCPt2cDwqjqupb9taTeXkCGlSIYSTgcuZOhBV3aRwOmxhIjFy2uurgel1LL4RWAf8O7Cmk+NSMcQY1wPfonI4qZbNwJWeFKE6/Qb4PJW9z811LD8XeFdHRzSBwkSCXU97HW+YyuGC3wOvBxbGGN8TY3w8icGpED5K5Re6li7gKwmNRTkXY1wTY3wNlYt+3wOsAjZMsVrip8MW4hTYKU573Try7/XAZTHGnyU2MBXKyHZ2D7Bkgm8PA9+MMZ6V7KhUFCGELqCPyt7CkVSOiow/AzUC307ydNii7ElMdNrrRiqnIX4YOCjG+FIDoVZMMYHthLVaEmMcjjH+vxjjM4GnAV9i10NRgZHTYZMaV1H2JG4GnkPlr7mtwHLgUip/2W1PcWgqmBDCfCqHBcZfM3EfsNT5CLXTyPZ2HvAvVK7MngtsAz4cY3x3ImPI+zYdQjgIuBd4HLiOyiElzy5Rx4QQrgbOYuee+GbgnTHGZemNSkUWQuimcijqIuBYKnMXC5KYUy1CJKYDrwH+N42rEVU+IzeJvJmd99zZCizyAjolIYTwFOAZwFeS2HPNfSSkpI2bwHbCWoVWlIlrKTHjJrCdsFahuSchNWFkQnE18CBOWKvAjIQkqSYPN0mSakrseRKhf2AucDZwPJPfyhsqF8HdBlwT+3q3TrGs1FGhf2B3Ktvusew8o6mWR4HvA1+Lfb2DHR6aNKnQP7AI+EcqF+fNnGLxCbfdRCIR+ge6gM8BRzWw2jOp3Pb7tR0ZlFSH0D8wHfgy8KQGVnsWlW39wo4MSqpD6B/YE/galXtD1WuXbTepw01Po7FAVD079A8c3O7BSA14Fo0FouqUkb/ipLScQmOB2LHe6G03qUg8JaV1pVY1u/0FmouL1C5PbnK9MHrdpCJRzzMeakn8SUzSKK1su62sK7WqLdtuYhPXu3jJ4kPGfL19W+D5Z63jnz+6OqURSfXZvjVw5Vv25rc/nsPGx7pZeMB2znnnAMe/uJ7Hm0rpevDeaSx720LuuXM206ZHjnnhBt58xWqmTdyU9CJx7f137/j/zRsDZx96MM8+faoHbkjpGxyE3n0H+fdv38+iJwzyo+vm8OE37cuBh97Hfks8o0nZtuxtC5m/1xBX/+4vbHi0i4vOOIBvfXJ3zrxgwnuPZeM6iVu+OY95ew5y1ImTPfVLyoaeuZHzLl7LfksG6eqGE07bxIJ9t/OnO8bfPlzKnjUrpvPs0zYwc3akd98hjnz2Ju7/U83TY7MRiZu/uRsnnv4YIRvDkRqydmU3q+6fwUGH+uwSZd+pr36UH1w7jy2bAg8/MI1f3zqHp59c81Bp+p/KK5dP44939PDCcx5LeyhSwx7fDh987SKefdpjHPRUI6HsO/KELTxw90xetvQQzn3aEpYctpWTXrKx1uLpR+KGr+7GE4/cwv5LO/7wDKmthofg0vMWMW165IKPPJz2cKQpDQ/Be8/an2NftIFvLb+ba35/DxvXd/GpixbUWiX9SPzg2vk892Xr0x6G1JA4DB86fx/WD0zj4qsfYrpnaisH1q/t5pGHp3HGG9cxY1Zk9wXDPP+sx/jVD2rebibdSNx52yweXT2N557pWU3KlyvevJAVf5nBB76xglk93kpZ+bDH3kP07vc4//vp3Rl8HB57pIubvrEbi5+0rdYq6Z0CC3Dj1+bzjOdtYM5u/pIpPx66bxo3f2M+02ZEzn7qztvGvP6SVbzoHP/gUba966qH+PS79ubbn9mTrq4wIRNOAAABC0lEQVTIU47ewhsvq3l9WrqRePsnPY6r/Nn3oEGuX/PntIchNeVJT9/Gf3z3gXoXT39OQpKUWUlFopXDSR6KUprcdpVXbdl2k4pEK9dAeP2E0tTKHIPbrtLUlm03qUj8iOaq9jjw0zaPRWrEbU2utxH4dTsHIjXoh02utwn4VfWLRCIR+3pXAR8EhhpYbRB4T+zr9WwRpSb29d4DfAwYbmC17cBFsa/XK7CVpluBrze4znbgnaO33RBjcodNQ//AXsAxVJ5xHWosFqk84/rHsa/X3XVlQugfWEBl2+1h8m13HXB77OuteZsDKUmhf+AAKk8HnewGlDW33UQjIUnKF0+BlSTVZCQkSTX9f3sA4pkyxlNTAAAAAElFTkSuQmCC\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//AAAACXBIWXMAAA7EAAAOxAGVKw4bAAAc4ElEQVR4Ae2dUXLUuraGO7d4pqjDC8/NDHJhBNwhwBByzhB2MYJde88AGAEVZsDOCCCZAeeZlxtSTCDn/xuvfZSO7SV3t6Sl9K8qxbaW7fXrU/dqRbblk9VqdYp8iTyWPt3e3r4ZM6hMBERABERgGYGTk5NvOGI9dhRi7cmjxPAe69w5Tf9ON7QuAiIgAiKwF4E/cPSTrTO8xPZrlp0gWw/5OSK0AjCpKImACIhAJQLoNTMYn7OH/D+VfMqNCIiACIiAQ2BxQGY0R/6BvN3tdlztb27p21MfWVvv2iOzlTbv0zVuj8xtXPHd0lL60zHkux6TLThn8P2AzCENDkhXC8YtfaOesymytlnhMEbXHlmftHmfrnF7ZG7jiu+W1tLPMeRb5DXGMFZexn7/HPZ/4u17aHtL315dImvrXXtkttLmx4yxz19kbmN6t8sOqR/n4hgyXNyuFg9Z4EAlERABERCBAgQUkAtA1SlFQAREYBcCCsi7UNMxIiACIlCAgAJyAag6pQiIgAjsQkABeRdqOkYEREAEChBQQC4AVacUAREQgV0IKCDvQk3HiIAIiEABAgrIBaDqlCIgAiKwC4FdAnK1p/RGKtTS94icO0WRtd0ROrIRXXtkfdI28oHKKIrMLUN+maeVsx6dpjo8NvgOi38g/x+3kS5Q9hXLSzxhwqk7i6WWvr1KRdbWu/bIbKXN+3SN2yNzG1d8t7S0fk2/eZe3tkRABESgKgEEeU2/WZW4nImACIhABoFdxpAzTqtdREAEREAElhJQQF5KTPuLgAiIQCECCsiFwOq0IiACIrCUgALyUmLaXwREQAQKEVBALgRWpxUBERCBpQQUkJcS0/4iIAIiUIiAAnIhsDqtCIiACCwloIC8lJj2FwEREIFCBBSQC4HVaUVABERgKYHsuSx4Yjzi9wcW/498g/wc+SPmsbjCUqlTAtHbNLq+TptdsoMSyO4h44txiTp8QQD+c5hM6Hdsn6N8Xatu8PUa+QdydzNFRdQOTc3bdO6zE10ftUdsV2MaWZtpnFoeq/asHjLg/BPg1gjEnwwg1m9Qzu1z5P+18kMv4YPB9wPyv5EZ/LsJxpG1t2xTtKGbIusL3q76vrifrsPvcMjPxCnk3SIz4K7GMmyfkc+3bSjjLEU89sm2rcQ2/PCHoZq/Q9YhmnboCdGmU4yj6zPd0drVdHEZWVuqc2z9mLSjrps4Sg65QxY2BzIbOU03w8aUPd1X67EITLVZlDaNri9Wa0rNgyDgBuShK87KXs/UuNo48owGmTIJRG/T6PoyMWs3EVhMwA3IOCPfEsJkPadfW7/+zgXpdD+txyIQvU2j64vVmlLzYAjkBGSrLC8WTKWnUwaVhyYQvU2j6wvduBLXH4FHGZLnesHWk+G9yQ8uDf86Xyys2BkG56Pfmx29TaPrW/iROI7de/6+RNHOgPx4+LjY8s6nB8GFt7exzIJvarceDG9Je3CJdUelit3S1wpY9DaNrq9Vu0X32/P3pbH2Z9a2HLL4OWzY0mzp8i9sWPBNyy1I067UF4HobRpdX1+tLbWRCXw3cbljyHz444UdlCz5+PTV8OuSFGu1AwLR2zS6vg6aWBJ7I5AVkBFw36Ni1xi64A3MmzSMufBBjbOhqMZirJdew+8hfITSHqhNR9lG15eIDtWuiS6uRta2JfXe5lFqz7moZ6Q4lvoWgXg9FLzE8hW+OFe2Q6klfL7DuTk8Yg8LXKDsK7Yvhy9uKdd7nze49mZtmgk2rL7I7RpZm9fux66dV+tOkTnJzHMEtwd5cc77EMguAiIgAq0I4EeIIw+cmuIka8iilVD5FQEREIFjIqCAfEytrbqKgAiEJqCAHLp5JE4EROCYCCggH1Nrq64iIAKhCSggh24eiRMBETgmAgrIx9TaqqsIiEBoAgzIvNXtN+S5CV1CV0LiREAERKBjAnyWgzF488aQNZZ8m7TNS8FyJREQAREQgToE+CwIY3D2K5zqyJIXERABEThiAhpDPuLGV9VFQARiEVgyl8UKj/ixW83J6G+QOdPbxxpzWcCPUiEC0ds0ur5CzaLTHimB7B4yvhic7+ILAvCfyJz97Xfkc5RzDLpKgq/XyD+Qw80EFVnbVONAc/M2ndLG8uj6Bo36TM414oStx+9LWpVS+rN6yHDOaTbXCMSfTBTW+SYRbnPe2mJv1YAPBt8PyLwbhME/TDCOrA2cZlPLNp0VNhgj64vc7tKW8+nabZ9abHmF7xaZAXc1lmH7jMzZiO7YUcZZinjsk21biW344Q9DNX9L6hBZ21g9oDdEm45pY1l0faY7crtL2914ZW12iOUh2eJcmzhKXblDFjYPMY69kziWzDRl/2XV34gEptosSptG1xexTaWpcwJuQB666azm3IMj1caRO+cdQn70No2uL0QjSsSDJOAGZNTaHhixnlMKYi5Ip/tpPRaB6G0aXV+s1pSaB0MgJyBbZecupj21nbTsikD0No2ur6vGltj4BHIC8lwv2HoyvDdZqR8C0ds0ur5+WlpKuyLAgPx4UGzLOxXAlT8bqrDgm9qtB6N38aVUgq9Hb9Po+oI3r+T1R+CZSWZA/jls2NJs6fIvbFjwTcstSNOu1BeB6G0aXV9frS21kQl8N3E5Qxbclw9/vLCDkiUfn75KejSJSavBCURv0+j6gjev5PVIICsgI+DyUelr3I7EG5g3abg1iQ9qnA1FNRZjvfQafnN8RNZ2T3+gNr2njQXR9SWiI7e7tCUNdeDVImwfLRDJx6PfIhCvh2NeYvkKX5yrBefYaVf4fIcDOTxiDwtcoOwrti+HL+5O5z3EQZG1ZdSvWZtmaOMuYfVFbndpy/x07bBbabYn0HSKzElmniO46eLcDo2kQ0RABERgVwII8hx54NQUJ1lDFrs60nEiIAIiIAL5BBSQ81lpTxEQAREoSkABuShenVwEREAE8gkoIOez0p4iIAIiUJSAAnJRvDq5CIiACOQTUEDOZ6U9RUAERKAoAQZk3ur2G/LchC5FRejkIiACInDEBPgsB2Pw5o0hayz5Nmmbl4LlSiIgAiIgAnUI8FkQxuDsVzjVkSUvIiACInDEBDSGfMSNr6qLgAjEIrBkLosVHvFjt5qT0d8gc6a3jzXmsoCfpr7pfy615DKnSzYREIG+CGT3kBF0ON/FFwTgP5E5+9vvyOco5xh00dTSt1exyNo87bl21PE18g/kIjNceTpa+/f0zdlbam/pe45Jji269lL6snrIcM5pNtcIxJ8MJtZvUM5tzlvLWbmKpJa+vQpF1uZp9+yoG4PvB2TehcMf3arBuLV/1Hfn1FJ7S987AxsOjK69lj5e4btFZsDlPLT3MmyfkTkb0R0byjhLEY99sm071DbO3cy3V4fI2jztS+yoJ3+Qi7bznJ7W/ue0ebaW2lv69rh49ujaD6kP59rEUTLJHbKweYhx7J3EsWSmKfsv635/p85dw7enPLI2T7vsIiACwQi4AXnoplP23IMjRcaRW/r22imyNk+77CIgAjEJuAEZsu2BEeuRpjWZC9Lpfruut/TtaY6szdMuuwiIQEACOQHZZM9d1HlqOxVatvTtVSmyNk+77CIgAoEIPMrQMtcLtl4i700ukVr69uoTWRvv2+YPxYVXiS37GS4sXG2VHd2m2O3W5OK2G7f0KAbkx0OBLVM776rg7W0ss+Cb2q13WORdfC19p5UcW4+sjXqpD4tityOOMXkoZWK3W0uK227ccNQzO5JDFj+HDVuaLV3+hQ0Lvmm5BWnaS6WWvr06RdbmaZddBEQgBoHvJiN3DJkPf7ywg5IlH5++Gn4Zk+KDrrb07VUksjZPu+wiIALBCGQFZARcPip9jaEL3sC8ScN4ER8YOBuKiixa+vYqFFmbp32hfey/o4Wn2Gv31v73Ed9Se0vf+zDjsdG1F9GXc1HPwHI88i0C8XooeInlKwSlGheBWvq2+k8tI2ub0pxVjrZ+hx05LGUPwFyg7Cu2L4cfo6zz7LpTa/+76uZxLbW39L0Ps9bccrSXZsurdafInDjoOb5kRS7O5VRU+4iACIjAMRJAkOfIA6emOMkasjhGSKqzCIiACNQmoIBcm7j8iYAIiMAEAQXkCTAqFgEREIHaBBSQaxOXPxEQARGYIKCAPAFGxSIgAiJQm4ACcm3i8icCIiACEwQYkHmr22/I1xP7qFgEREAERKAcAT7LwRi8eWPIGku+TdrmpWC5kgiIgAiIQB0CfBaEMTj7FU51ZMmLCIiACBwxAY0hH3Hjq+oiIAKxCCyZy4LP57Nbzcnob5A509vHSnNZNPWNes6mllxmhWUYo2uPri8DsXYRgUUEOH7BV7yvEVw5sflohp3zXbw2O9Y529E35NnjbP99lvDRzLenO7K23rX3zNbYow6cp+AH8hMrq7Vs6durY2RtnnbPvrRuw/447e0qq4eMXgqn2WTg/YTlJmGdbxLhNucELvZmipa+f9V0+m9kbdOqf1mia4+ub44vtLOz8gGZdzDxonmRqRpx3nuppe97YrYKImvbkrp481B1yx1DfgOFY28F+YLy00HM4kpkHtDStycxsrbetXfLlp0V5DfIvJXps9cQh7S39O3VI7I2T7tnP1TdcgOyzYe7rYtjyUxT9l/W/f5OnbuGb095ZG29a++ZrcdedhEYJeAG5KT3O/fgCP8tO3hq6durTGRtvWvvma3HXnYRmCPgBmQcbA+MWI80Pd9ckE7323W9pW9Pc2RtvWvvma3HXnYRmCSQE5Dt4LkLE09tp0LLlr69KkXW1rv2ntl67GUXgXsEHt0ruV8w1wu2ngzvTS6RWvr26hNZW+/am7IdhkwuPIhb9jNc2LnaKtNmJwSitDkD8uOBmS3vIOTVQ4hlmQXf1G49mCLv4mvpO63k2HpkbWN607Lo2lvro3/wKnYrZ9oWWo9BoHGbPzMKHLL4OWzY0mzpkre8WfBNyy1Ij90Sl+63z3pL357uyNp6194zW4+97CKQEvhuG7ljyHz444UdlCz5+PTV8OuSFB90taVvryKRtfWuvWe2HnvZRWCUQFZARsB9j6OvMXTBx0A3aRhz4RN8Z0NRkUVL316FImvrXXvPbLfYj/1nubVLsc2Wvr1KRdbmaffsO9ct56KeOeeY2lsE4vVQ8BLLV/ji1LiQ0dK31X9qGVnblGYrj649uj7jeG+J78k7FHJIzx5wuUDZV2xfDj829445VEFL314dImvztHv2Q9SNV+tOkTl5z3N8UIpcnPMqIrsIiIAIHCsBBHKOPJwj/p5kDVkcKyjVWwREQARqElBArklbvkRABERghoAC8gwcmURABESgJgEF5Jq05UsEREAEZggoIM/AkUkEREAEahJQQK5JW75EQAREYIYAAzJvdeObDa5n9pNJBERABESgDAE+y8EYvGJAXiPzbdI2LwVWlURABERABCoR4LMgjMGbgFzJp9yIgAiIgAjMEdAY8hwd2URABESgIoElc1ms8Igfu9WcjP4GmTO9faw0l0VT36jnbGrJZVZYhjG69uj6MhBrFxFYRIDjF7fIawTX1VSGnfNdvDY71jmj0Tfk2eNs/32W8NHMt6c7srbetffM1tijDpyn4AfyEyurtWzpe986Rtd+SH3DuYDsdpXVQ0YvhdNsMvB+wnKTsM43iXCb89ZyVq4iqaVvr0KRtfWuvXO27Kx8QOYdTLxovvN0jDh2UQK3Zr4XCR3ZObr2Gvpyx5DfgN/YW0G+oPx0EDqC+CBFLX17FYisrXft3bJlZwX5DTJvZfrsNcQh7S1971uP6Npr6MsNyDan6zbzm6Fgyr69/y7bU+eu4dvTG1lb79p7Zuuxl10ERgm4ATnp/c49OMJ/yw6eWvr2KhNZW+/ae2brsZddBOYIuAEZB9sDI9YjTc83F6TT/XZdb+nb0xxZW+/ae2brsZddBCYJ5ARkO5gXC6bS0ynDgcpb+vaqEFlb79p7Zuuxl10E7hF4dK/kfsFcL9h6Mrw3uURq6durT2RtvWtvynYYMrnwIG7Zz3DR52qr7Kg2xW3/5mZAfjycxpZ3zsoriwDNMgu+qd16MEXexdfSd1rJsfXI2sb0pmXRtbfWR//gVexWzrQtHtK6uO3cms/sSA5Z/Bw2bGm2dMlb3iz4puUWpMduiUv322e9pW9Pd2RtvWvvma3HXnYRSAl8t43cMWQ+/PHCDkqWfHz6avhlTIoPutrSt1eRyNp6194zW4+97CIwSiArICPgvsfR1xi64GOgmzSMF/EJvrOhqMiipW+vQpG19a69Z7Zb7Mf+s9zapdhmS9/7Viq69iL6ci7qGViOqb1FIF4PBS+xfIUvTo0LGS19W/2nlpG1TWm28ujao+szjveW+J68QyGH9OwBlwuUfcX25fBjc++YQxW09L1vHaJrL62PV+tOkTl5z3N8UIpcnNu3kXS8CIiACDxUAgjyHHk4R/w9yRqyeKggVC8REAERiERAATlSa0iLCIjAURNQQD7q5lflRUAEIhFQQI7UGtIiAiJw1AQUkI+6+VV5ERCBSAQUkCO1hrSIgAgcNQEGZJvDwpZHDUSVFwEREIHKBBbPZVFZn9yJgAiIwNEQWDyXxdGQUUVFQAREoBUBjSG3Ii+/IiACIrBFYMlcFis84vcHjudk9DfInOntY6W5LJr6Rj1nU0sus8JkFAER6IpAdg8ZQYfzXXxBAP4TmbO//Y58jvJ16Rq39O3VLbI2T7vZUYfXyD+Qi8xgZX52XUbWJ227tWpkbjk1KqU/q4cM55xmc41A/MnEYp1vEuE2563lrFxFUkvfXoUia8vQzuD7AZkTSvFHNVQwBtuw+qQNn5YdUmRuOdWppZ+zvd0iM+CuxjJsn5E5G9EdO8o4SxGPfbJtO9Q2zt3Mt1eHyNo87akd9eAPbtF2TP0tXY+sT9ruxoTcto3MLacOh9SPc23iKP3mDlnYnK449k66Gbam7Hd23nFj6tw1fHuSI2vztMsuAiIQjIAbkIduOmVfz2gvMo7c0vdMXTemyNo87bKLgAjEJOAGZMi2F5lajzStyVyQTvfbdb2lb09zZG2edtlFQAQCEsgJyCabF1mm0tMpw4HKW/r2qhBZm6dddhEQgUAEGJB5lf035Kne7lQ5q2G9RN6bXCK19O3VJ7I2T7vsIiACcQhcQQpj8OaiHsd/+cCHBVeW/51w5c+GKsbs1jss8i6+lr7/BjCxElnbhGQVi4AIxCRwClmMwdl3WfyFfS348jhLFqRpL5Va+vbqFFmbp112ERCBYARyx5D58MeLEe18fPoq6S2O7LJ3UUvfnvjI2jztsouACAQjkBWQEXD5qPQ1bvXiDcybNNz2xQcKzoaiIouWvr0KRdbmad+yj/33s7VL083I+qRtt49GZG45NSqi/1GO52EfPh79FoGYY85ML5FfIShxQLp0aunbq1tkbbPa0ZbvsAOHnewBlwuUfcX25fBjM3t8aWNkfdK2W+tH5pZTo9L6TyCCA8qcOOg5voRFLs7lVFT7iIAIiMAxEkCQ58gDp6Y4yRqyOEZIqrMIiIAI1CaggFybuPyJgAiIwAQBBeQJMCoWAREQgdoEFJBrE5c/ERABEZggoIA8AUbFIiACIlCbgAJybeLyJwIiIAITBBiQHw82W07sqmIREAEREIECBJ7ZORmQfw4btjSbliIgAiIgAuUJfDcXGrIwElqKgAiIQGMCCsiNG0DuRUAERMAILJnLYoVH/DhnJyejv0HmTG8fK81l0dQ36jmbWnKZFZZhjK49ur4MxNpFBBYR4FwWfAX8GsF1NZVh53wXr82Odc529A159jjbf58lfDTz7emOrK137T2whUbOQ/AD+YnHu7Y9sjaPxTFpH+oKJLerrB4yeimcZpOB9xOWm4T1G5Rzm3MCc8azIqmlb69CkbX1rj0yW2hjZ+QDMifjWiMXmYoR512cImvzKiPt+W8MeQOYY28F+YLy0wGkx3tXe0vfnubI2nrXHpYtOyPIb5D5HrTPHuia9sjaPA7Snh+Qbb7cbaYcS2aasv+y7vd36tw1fHvKI2vrXXvPbD32sovAKAH3Louk93s9eoZfhfy37eCppW+vMpG19a69Z7Yee9lFYI6AG5Bx8D+GE1iPND3fXJBO99t1vaVvT3Nkbb1r75mtx152EZgkkBOQ7eC5CxdPbadCy5a+vSpF1ta79p7ZeuxlF4F7BB6hhFeKeXFiqrc7Vc6TWU+G9yaXSC19e/WJrK137T2z9dg/WPsw1HSxsIJnuJh3tfCYg+/eWDvrzxi8ue2N47//QuYtbDfIdxKvfEIsyyz4pnbrwRR5F19L32klx9YjaxvTm5ZF1x5dX8pS6/8lwHbDVrFbYP/r6fBrjbWfokaMwX/akAWD8lziLW8WfNP9LEiP3RKX7rfPekvfnu7I2nrX3jNbj73sIrBNYBODLSBvG7e3+fDHi+1CbPPx6avh12XEfJCilr69CkTW1rv2ntl67GUXgVECWQEZAfc9jr7G0AUfE92kYcyFT/CdDUVFFi19exWKrK137R2xHfvP0cNfyx5Zm8fgKLXzol5u4tjQWwRiG954ie1X+OLUGJBv6dvjE1lb79rDssX34B3gcsjOHmC5QNlXbF8OPyYe+2L2yNq8Sh+7dl6t44AyJ+95jg9SkYtzXiPILgIiIALHSgA/Qhx5OEf8PckasjhWUKq3CIiACNQkoIBck7Z8iYAIiMAMAQXkGTgyiYAIiEBNAgrINWnLlwiIgAjMEFBAnoEjkwiIgAjUJKCAXJO2fImACIjADAEF5Bk4MomACIhATQIKyDVpy5cIiIAIzBBQQJ6BI5MIiIAI1CSggFyTtnyJgAiIwAyBJXNZrPCI3x84Fyejv0HmTG8fK81l0dQ36jmbWnKZFZZhjK49ur4MxNpFBBYR4FwWt8hrBNfVVIad8128NjvWORvTN+TZ42z/fZbw0cy3pzuytt6198AWGjkPwQ/kJx7v2nZpm45n+7bFIdkO54Kk280bQ7A9n9BL4TSbDLx8q8gmYZ1vEuE2563lrFxFUkvfXoUia+tde2S20MbOyAdkTsa1Rg4zVaS0oTUKpRpsc8eQ36COY28F+YLy00FoIQyrlr69OkXW1rv2sGzZGUF+g8z3oH32QNe0S1s52jXY5gZkm/N1u7YcS2aasv+y7vd36tw1fHvKI2vrXXvPbD32sovAKAE3ICe937k3AfPftoOnlr69ykTW1rv2ntl67GUXgTkCbkDGwfYiU+uRpuebC9Lpfruut/TtaY6srXftPbP12MsuApMEcgKyHTx34eKp7VRo2dK3V6XI2nrX3jNbj73sInCPQE5AnusFW0+G9yaXSC19e/WJrK137T2z9djLLgKTBNyAzCuLw9EWfNOTWQ+myLv4WvpOKzm2HlnbmN60LLr26PpSlloXgUMScAPy4Iy3vFnwTf1bkB67JS7db5/1lr493ZG19a69Z7Yee9lFYJRAbkDmwx8vRs7Ax6evkh7NyC57F7X07YmPrK137T2z9djLLgKjBLICMgLuexx9jduR+JjoJg23JvEJvrOhqMiipW+vQpG19a69I7Zj/zl6+GvZpa0c6SJsHy3Qy8ej3yIQr4djXmL5Cl+cqwXn2HXXlr49zZG19a49LFt8D94BLofs7AGWC5R9xfbl8GPisS9ml7ZiaDnJWdF2P4H0U2RO3vMcH6QiF+fK4dGZRUAERKBvAgjyHHk4R/w9yRqy6Lu6Ui8CIiACfRBQQO6jnaRSBETgCAgoIB9BI6uKIiACfRBQQO6jnaRSBETgCAgoIB9BI6uKIiACfRBIA/I3XO273cq8OV9JBERABETgAAQQX+/FWZz27zjL+5B5qxvffDCWatxjPOZXZSIgAiLwEAnwRdGTD5X8B3IxDIy8iCo+AAAAAElFTkSuQmCC\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": "iVBORw0KGgoAAAANSUhEUgAABbsAAAA1CAYAAABofTaAAAAACXBIWXMAAA7EAAAOxAGVKw4bAAAgAElEQVR4Ae1di9XcNBMlOSkgQAWEDgJ/BUAHPCqAdBAOFXBCByEVQOggoQIeHSRUQPJ1wH+vP81G67VlydbL3qtzvH5IlmbuXGvssay9899//72nJASEgBAQAkJACAgBISAEhIAQEAJCQAgIASEgBISAEEhF4M6dO09wzn0sD7C8xvI94o03WB8+XbPuvRr3joLdvZpGcgkBISAEhIAQEAJCQAgIASEgBISAEBACQkAICIF+EXDB3qeILzLI/R72n2P1APuf9Ct1HsmuWfc8CJap5W6ZalWrEBACQkAICAEhIASEgBAQAkJACAgBISAEhIAQEAIHR+DLkX4/Yv8hAsEc5X30dM26d2vbuyDfWyz/jRZ+fqAkBISAEBACQkAICAEhIASEgBAQAkJACAgBISAEdo8A4l5PsHy+e0X6VMAPbNv0Jf6xPqXOI5Wv57XpvhlBXJP3sfBrgOTkrulxTPvtPVfTV1j/7dX6xtvWphAQAkJACAgBISAEhIAQEAJCQAgIASEgBISAENglAgyKUXBMrfFylwp0LDQw/Xgk3kO3/+fo+OF2r1n3XMYEhje4Pp9i+QvbqVPf8CuCp54s5N4zC3a/RoXD3DpeAW0KASEgBISAEBACQkAICAEhIASEgBAQAkJACAiB3SKAINpjCM85pDnQU6k8AnyxcDV/UDmC85p1H0ERv4tr8+WagDfO40h6G03P+eLvs9W78U2rpBAQAkJACAgBISAEhIAQEAJCQAgIASEgBISAENgHAgh+cU7lHxTormMv4M1g70vg/VOdFvtp5Zp1z2EFcOZn1PPG4bipSgW7N8Gnk4WAEBACQkAICAEhIASEgBAQAkJACAgBISAEekMAQTOO8nyGRSO6KxgHeH+HZu4jaPmoQnNdNXHNumc2BK/V74Dnprn1FezObBVVJwSEgBAQAkJACAgBISAEhIAQEAJCQAgIASHQHAH+6d2vCL5qnu7CpnDByY8t0M19LDZ3d+HW21Z/zbrnRh784ZQk32JZ9YeVJo+C3YaE1kJg5wigg+XbL75JVRIC0QiAM6v/+Ti6ERUUAkJACOwEAfnSnRgqQcy9+DlxL8GoKnqGgDh+Bod2KiLQO/cgH6cv4ejQ7yvCcpVNAWsGtTl9yR/E3WFP3A//34DXrHspsiPg/Rvq3jSdyR1U8BbLZ6js71KCql4hIATKIoAOln+48Q2u49R/ri0rmGo/IQAb8UbrCyz/Yvkflh976XedbE/EH1hFSQgIgatF4Jp8acgnIY+ffPPl+cdYXsE3nObcdHn/4PhHOM6RN7tITt9u/Rzk033cDpjkeKR7uRW2uiaOh3iCPPWvK/iz9hTg/Qrncu7oq5tSIwaznFxFXYwrkt9nCdgz5thdulbdQ3r3ZiTIypdVHN3NrwWiX5rgPL54+V3B7t4sKnmEQCICrsNiJ7CrB89ENXddHDbiW+7T3GVmM3Ta7/eiGGTig/b/IJPms+vFKJJDCAiBaghYv4wGD+9Ll3wS8p/CFzwyTHxfhWP0FQwanz28uuMfIo8vdD9Efnej6JyM3fk5w/kauAcdd5tgJ93LrbTeNXF8iSfIT+5fCTvOG0bNom/ly5buEuTr7jkCMvGl7VMsSYGy7sAtJBDwCfZpyF/F1ULiZq32WnVf0jsryJkqg8zJL6xwzhDs1jQmmYygaoRACwRwIfPtKQPd/DpjNyOsWmDVqk3YiG8kv4R9/BEFf+IYpw/Z9KcLOXWCfBy59wAyaSqcnMCqLiEgBLpH4Jp86ZJPQv4DGOwvZzS+/ByPpGGg5WzeU+c3hgC38yX8hHnTPIuu/ayrHv0ccNJ9XFYrl6kMdtK93Epor4njSzxB/pr+9SHOY8CWzxE8v8vUY/8KoPjSlaO6x36sSwxrClWCqzXl39LWteq+pPcWTAufy5cynK734quBpXYV7F5CSPlCoG8E+M/SdOKbpyFCB/IEy1MsL9w6uUNpCVXH8tNG7KT99Knb6Q1jBjbIg97k8rHTthAQAkIgNwLZfGmKYI381pJPeoN7ip+dHl9j/ctIJ/qvF6NjDCicjuF8zrPI+Tp79CW9+bms3GvEqREd1u92LP/SdbNe6fxnHpbjHfPDrLjEk+T+lc94WBjo7u4FointrbvhHrjCkZ18OTB+BvPEXbdZgocl6lzQLjtX2R70eOx0GdYLMrTKLqK7058vp160Umyh3SW9F05vlv2razl5QN69ZiKrYSEgBDYhgI6UTpwjTTin5qZEp4QK+KnS8OYb+7yh+h3LLuYA71V+yMWR23zYt04am0Oi7Zg2v6S4rSbPL+0PmTlij86QN6xKQkAICIFDI4A+L5svTQGqhd+K8UnwAzfUwytrgW/7jJ4+7TSyG+W4z4DCGyx+Yj30gQx8d5N68nPALiv3UJ/u5QowzbsWdC+XiG9OjvfO7xiepPaviXA3L95T/wow+ILgBjKd/FUOgErwsESdIV1LcRX1Mhh5msYM+3zp/Rw26OaZsqDu9Of2FXd3X2DE6B3iTMs89puQnzGTb7DwS/TopJHd0VCpoBDoDgEGJH/jjUUGyRg099OP2OGbye46a19Ib7tX+encX7OT9mTlJjtrHs9hu1HVm3dp+15H5G1WThUIASEgBEYI5PSlo6qDuy38VopP4nQlY//F4DWDB/6LWrtPGPs5Br8/CCLQLrMXP5ebey04ldOKvcqfct3kxGNLXUfkeK/8MDul8CS2f7W697TuhXv8MilroNsZoQQPS9QZ4kwpru7hK68iuvO+CMsjgN7rFxgpeoe40yqPXxkmx6YU7G5lLrUrBDYg4ILQfINIp5Ir2QMr67OHVv9YrnZK1ePL2ov8DAyc3WjBdrQbFzqd7pILYjCQkfypUHfKSCAhIASEQACBQr400OJFVm2/leKT6Kf8oDaFZ4CG/znBUd7j+RMZ3B6n++MDPez34OcKcq82p3KbtEf5U66b3Hisqu/AHO+RH2ajFJ6k9q/WRvfrTrhHfOl/XhQCrAQPS9Q5p352rsKnEW/qML4X4DM52+slZde9F8UW5EjRe6GqJtn2lWASl+42EVWNCgEhsBUBBrn5BjHLyGDUw3+p9oOyvElgGh5qbzf7/e1Rfs/pn25e3DG+8X0EmcdBhJ4A5h/h/NCTQJJFCAgBIVAAgay+NEW+2n5rhU86u7/A+XzA+BSLHed9Ax9iLX1gG2493h9lN99t7eeyc682p3JbsEf5V1w3uWHZUt+hON4jP8w4K3hi/ehQRUT/ak3tZd2aexYQ859ts2BXgocl6pxTtiBX7XnXvy+gGAx+d3E/UFD3Obi7OL5C7y7k9oXANcI+k9zioIvodC+6JAoCKJL4czR2mr8v5fwjlwU2/PwkW/AxFSvZJhWx7eUb25yfZvEzsVKJ8z5+j2t97LBKtZe73h7ktxutR+DKY6cg51f/Crj2HOimqLw5vA+52d9nv1FkA6Gk/iyETrm8ln2abF7Orj3X3JJzDpfSvjQF/tJ+K9UnMRj7DDYy/8UHDf6PB//EmMf4SSmTBW3u3+6efrnfs69r6ueATQ3ulebUydiFNnqQP/W6KQTFqmqPzvEe+GGGSeVJbP9q9e9t3Zp7DIhxyi3zTyXxK8HDEnUaBqW5yuD2OI3vD8b5tfZL615Lj9R2UvVOrb9WeQ7CtAGZUW3eQam3WD5bCr7gxpYk5QTzSdH0KCkOUggYccQmA4Q1OtYTarLNCYrqGy1sjjb50ukVlk+Wrts1gKB+Otj7qJvzTu0u9SI/5OCohk+B4y7+5HNsaMhP3/Az5OcNebWEduVrqqF92RDwr+7HZPNLO1zTkRacI75ot6gvTbEhZCnud9FGMZ+EunlPwi+WTi9Hcew/HHsfx25SsKhZFjK28nPFuQfdinOqpK16kR9yFLtuSuJndR+V473ww8O5OE+gMwNWT9GncuBM96kV9wiMa/tPYFU0blWChyXq9MmC+otwFfUyCPkXFn71dYqF4Tj97I84lvTHgr7MubZL6e7Lhza6u05r6O1jUGobevC+5jG4xBh2MKEs+fj73WCp88zfsTsb/EKF/EOzt1h6eXNzLn2dvW/RTItJ6YO2qaP61bbSwubD2zlc6NlHTOH65TzNew509yQ/R239suMrg0ELexNcUw31ZzXRvmyrRZ8mm1/a4ZqOtOAc8S3mS1OMV9HvlvRJfHg+/Q8FdOLXjvwD7W4D3c5GrfxcUe5V5FQK1aPLdiZ/yesmGpMNBQ/H8c74YaapwZO9xViacA/8IE5cTgFXM1LOdQkelqhzQudSXDW8xzzlfvaYxYReMYdK6R7TdssyR9H7FUHEdRI9ujsq2O0uPL4dMxIPxsJxfub+HAuj7N9gGZN7KHctP+6m/iXweFxLZ7TF4N6FbWq1f+3ttLA5MOdI4bNrMYcdwCU+fPFt7PBSi/tYojuTHDJsqSOn/KiLL+9WX8c4l6O22B/yRm+v6Q8IXtX+wE39WWO21O7TZPPGBu+g+dqc81Qu4ku9+hc3wf8qfre0T4INOWLrFdrhFCd8JvgCx07B70Ug2hWo7uecqsW4V4tTpUyWU37UpXu59947FMdz8iMXhyFT0Xt+1o+FLxTZt3KbsZfVzyi59I6opxX3aA+mITB2u5n3F/hn990l6hxrjTaKcdXdyzE+8cG4Xexv+g8wyL2pL6c8JXWf0Dfboa2671XvGQCNR8OAgZkyZ4fvne3N7/Az9i/G2Y7Uw80sgGSQgiM5rj1xHuV/sNT6VGPSNtduhMr617Y5L/CswW5cvwxq8ibqR3aqDj8GvffwsGpv+HqSnzcTr9FH9vIm25k0aTVwjE4SemTlW0AK9WcBcCpm1ezTZPOKhu24qZqcMxiy+1KrOGZd2e8W90nwE7Xue2PgjS3Tws9RtiLcq8ypWIyjy3Uof/HrJhqc9QUPw/EO+WFWKcoTdw/OZ7K9pVbcoz2Yijy7lOBhiTpvIbj4LcpVtGZfeQ2DvaBXT195ldbdwO5t8G8tvU3/kmu7pqOncloMdoOkQ4dRMdhREqDidQOnG2D2GkvxP3aTbYqbM6qBmjZ3AvGazD1imNMIsHN+7toYVtTN3+94uyv5gRvtE90Rd4qrORS+CLHtYqKqPysGbXLFtfo02TzZNIc9oRbnRgCW8KWjJoK71fzWQXxSEMyVmebbqvg5T8ZS3KvGKU+XnJtdyX+Q6+ZIHO+KH0b8g/DE1Mm5bsU9G1ls7efUiXWV4GGJOi/0Ls1V1P8T7u0fY+EANCZOjdrFwLnSurtnGg7g4cvs4QsMrP8gJlg3S6X1rqkYdGGclU3aNb7Y/GKwGzXwjUzuwNqiYEsFoCiJxNHm/2L5HxZOfJ9lFCXq5h9xzl6Yjswv0ObZn/Fg3xLx4vmlcZu1DWRk4JKj7RlwewV9Theak59/IPARSYN11oT6Z22zIBdl5qj4InKFlIRce7G5qfHKNubWC1jzQevEAfDg/bl6Sh3PyZMW8pfCZWu9GbhsIthNYrRDsRNXrmf7M9bn8Zm77PN5Q0EeM0jxIZZ//X4O+9lSTq5mE2qhogw8qOHHstnc48fhfd6Crs38KCm5kXc1ODe+clJ8Kc8N9j0Ltmnud8fK73l/I9dM9dp+ztrlepF7LORxirtB/rW4F4J82e75W8hPUHtMmfhN1brmeG5+e/VVvxfokUdrZDoA93gfxBQd4/B4w/Oy9LNenYtcPFLfV+pZjIbpOUFv9rV7/AKjGKwZ+5KxjHaNj49P7vMfUh/CQO9NLchjUPe7qTz/GMtg4T+v8w3OZF25jqMNvi3ivxEP7WCbN1pvbX/LGvVQD47KntUB+QSYuvLfQC/K4Tjr+GsqL+cxtDFrG+Q9YVtIF9jgGOf5Qval7FuPod6gbZA/2A3rJLlQngGtF1vlmzof9e7J5sa9mGuyCQemMB4fA+ZFeDJuZ+0+5GMQbvL6XltnjfNycNmXE/XN9nN+uRzbaCvUn5H3A5/ZFhKvWQYKBhu57SJ9LurOzlXK7eod1jnw8+tA3bvo0yBnNpsTT8eNJN/i47Zm29lx9n4E+ck+D+cE+bGmTodNMT/q6t/EO+jF84tcx2Pbop0oX+rKRfc9KN+Eh2P9jr4PnDdxzccHdVXzc2wXKYp7Xtlo/vl61diGLtn9Y065IZ/u5W451yXH3bWQld/kpLt2qt4L5ORty7qAX7a+1dmhKvdcm+yX2G5UTIrlsGTlocnh1uLixpgT7LPLvpz237rsVXfInbUvMRxRL2PXi88KKMNnnrf38LOUOBLlzVKhWvl4Q0CyfwmF+ZbMEicr559l5pg65BPU/bNVPLVGPofQc6S0vS0fFyNe7DhLp0nbQDYe559CMH2FZSwnR8RnH3W+ZBsnB0cTM0XJhTpJVHtLRr1KpD3Z/IMYAFpxIFK24DWMOsjXJJ7EtHslZXJweQwVR03XSJP9mWv4B6w5p68l9v8sb301+4hxP2dlV68L9Wl0/h/Cj3BkOkft8U9Xgl+WrFAgBw9q+LEsNm/V37Xgh9M1qX/EOTX8KGm6lXc1OGeXU5QvReHovqcVD02hK1tv5doYrlp+ju3Gco9lo/nHwjVTif6vpvydt5Wb31S3R45n5bf64Cys3jv3CEJqDCYrDymAuEgUlK4cgRJ9iUEafY3HBLt5U3ZjNXewfgYZhiCBJ8unbjtace/c0yY6Jr55O5uzmJk4zrcHn5wK3m5wVOHfo2O2y6BLqcCstcH1nG3eQN7fXMGvsfaDRDxMPccY8vjWtGQbymXBqSm5aMczWVGeGD9ytqHcWdMObW4cZ1AglFpxICST5WXniVUMez7ENt/OX/yhrpXx1yjPt/9TvCLOHyD/G7+82+YfT/JlzWzCecPwrdkCXgbqGiaf8g6t2kSbufqvcfvGufHx3Ptz/Rnb4TRVN16DtPNLO4Z19pd3rq0SXGXfy+D8kCD7b7DdcywcgeLraEWS1hl5UMOP5bJ5q/6uBT+686MkaCbe1eCcXU/Wry350pS+pxUPTafZNewT7ZNmK9mQgb4ti5+jCJm4NtbG+DA+XmLf2lriHttO4V8JWUN1luj/aF9+2cLA7L9cgzuLzys4R/dyIUvd5hnvlktuL2FtLXE8N7+b9MFH6V8L9a1kk/FhO7PiauC9JUfU3sQVL9LPNuHikr69c7VUX05crlX3FL1z3auhzVIxCZqSfmW4xrmzlO4tFegp3wHHDvPXkVwMfDD9fbta/cvpXE5zW7MWtMm6OXJ8nD5GWQvcjvOa7lvn7uF1ktPpQ/myBoe8tmZtEyEXbZtVrghDHMLmYz0jsOYptbG2h9TgNRwh+wVPHK8tgBj9ogltTT5EoT6OPn+A/LP+YIzz3D7Oy/ZgP9fGxPFDcpl6Gic8nelEJ23nldm0WaJPQ53k7tRo5hscp072knKL7IfgQYrNraxns136vBh+ROh60T9uIVPCuYfg3Vhfw9s7Ptv3WNmaPPTkCm5CthY+KSjThsxDcm0KD+OUlzfLP69M8U2P41nv+VHvqi+fgNPk/QDq071ccTasbyA3v60+j5/F7wWoPdo9Sv96NX2rz1rjjXdscz9rddbmoqfD5GbvXIV8RfpygnGtujfSu5u+5O7klXB+kNFzPjz1kIapL2A0Bgb89A12ONpy62fs/tQoVj87vBe2wzU6rqVgGt82bJXFb3Jue8k2xOvvEV7U5wbHtr4YGMuUYhuOuqW9fDuWkmss53h/bzY3zGLfaNXkwBjbqf0iPCGfsTDYffFlxpQQBz2Wi8sDPOjnrN83zpWGbak/M7ns5WbplzUluGq+Y4wpdY+9ppfskIsHNfxYbpvX7O9a86MnP0pO5uBdDc7Z9WPXYPR1hz7xoTt5qe+pyUPT55rWObg24NXAz7HdZO7xpAT+sXjpVKL/o8wMdJyeuXBfxxfAnOrrPjOvJGXjN/HysDPe1YDR2irRv8bIrz44BqXLMkfgHrXivaXPfe5GpQL9rLgYhbwKHQyBrH3JBDbDNT5x/OJQTLCbQdtoZ3XRQt4DDIiePWS4TokPIOxMtqYHXidnN5b8NG48JQIDak8DjTGgYY4+UGxz1pJtKMc46E5dzjDcLMVtBSm2ob3+HrVLuYYR9LDBd1hq3djuzeZ2ccfiU5MDI5NO7u6VJ5PKdHYwF5dNLev3OWVTjTTbn6E/8B92ySG+oT/1Ich/XKDPKMlVu459XGOvaf+cqe1cPKjhx3LbvGZ/15ofPflR8jAH72pwzq4ZuwaD193KvqcmD02fa1rn4JrhVdvPsd0o7rHgSv7x1NIpe//nfDivHcPHdLjBxuD37cDB1zn5Tai65XhBfqsPXneRHIF71Jx9BpNx/3Zv5rcgD9miuDiDuw4fGoHcfYkPFq9ru8b945PbdyePnh/8G7vj+arPS9zuzT4w8AYGy39Y7E+Vps4PHmMdKMAOg8uQ3DEGnh+NAh9r26Ouv6PeJ1gYzP4dC4PoDL4+dcf5D6BP0d44iIzDp0S8ggFl1LVWxlMj2FiyzZmMaJM3i1xOoyZY2VZZeD6qibIN20OakovzddtxThETTeKN8lex+UYZB9D44+HC+QxjkmE6lIUckxyIqWirDjwf7TTjSYyOrctsxDgXlw0G2ovpze1q/nej3FbxZH+GuslZ9vNcM529fEQ+OcWvVU59xlZ5eD7qLMnV8Q342f5G+XPxIOjHNspIOzJls/ltdScfMuxCxsn+bqvsPB8NNONHQNdWfpQi5eBdkHNsZKvtWAeT11/M+lLHn6S+57b2OB66ssFVLn2DjVTOzKBTDq6Z1rX9XBT3KNwG/pluk+ut+PN8VFyi/2OdTDe3q9PvG2yd+chTTqcbGzHOyW8i1CXHS/HbUSLq2Wejnbpk30adds89ZxTO989k3L/dm/gtzEO2GMXFCdEuDm207UV9OiAEQghs5FvuvmQs6tl1Nc709+/5OzPbDIzOjmIGEMzjTYgFIhgs5ghd/qnjMFcWHypwjEIxwLv2T7isfv5ZIf+8hIlD5L9C/QT0lDa0xz9HfIiF9bNO/ovoa7THz+p+wMLEP9M4zQF2e+jilw+c314c9Q5skNGrZQhaz9oGBSn3sxFePP8sEJ9BlmjbsHGksVzkBh9y+ZKB2P+CJTptlL+KzTfKOMaCDwL2UDDOG++PsbbPSs44MD5paj+DDk15MqVTb8c2YpyLywaLPVyOHzwt/7TeKLfVM+dr6E/4KTPfErN/OL2AxPbwAhXtn/XJGeQpxVVzzuMbcO6f/NhG+XPxIOjHNsoIdYeUzeauvqj+LoPsTfkxo2tLP0qRcvAuyDk2ksF2rMbSki9N7ntcxVE8NCFC68z6hpqqlpdBpxxcM31r+zlrd4l7LLeWf9bG5DoD/qX6P5P3jW1467HP9LL629yIcU5+E5xeOV6E344NUX3wRjv1RzxItFGnI3GP9jHuh2xVkodsN4qLIQEtb6NtrRqthUAUAhv5lrsv8WXm/cDUfYJf5mybI5U5iTg7yMkF+fyMnX/WNpkfexx1fImFwe7kenAeg7oMoEefu6W9lHb8stQPyyv/WGh7q4xsC0uUbVAuiOFaWZbqDekfm4c2eHMdxHWt/LEyzJVDu9E2zyEj6mCAKOlaMNlxXpADVi60XqtDjrZDcjEPaZEnS3W4ethXPY4pGyqDOvhSglMh2ULbfRk6x2t/VV+5VHdsPuTkH0UR1Gg5qFtK+bEsODe6PxufO7W/Vh6ct/k6mZKHx5yOn/v5ODaJM45vwtNvI2Ub7Vbr0xweUT4sRQeWRQraEfmr8F2qN1VOvzzq5jUQxQ//vJht1sv6Q2WRvwqTUJ0xeWg3mnOsL4ecqGO1L43RycqgnSAPrVxonUPfUP1TeWhzlf+aqmvqWAudxnJAhup+jjIgVeHeWF9/fy3+OG8zn305bBv1PsRCcM78Afb5nLrqfgznsT9bda7JxTXSqmvBtR99D+W3mWsbMlwtxw1DYBDkbAs7oc1VnDKdltYtdBrLBBlacY/XPS/cs3upsXwt9iFTkIsxMqGOqvdpaK8oVxN03tyXx7RlZXrQm7I4e2/WHfXwPpv8s5gEv1w887emu79Gmap889seb0MW6kBQnozzxvsow3uKt3fxE5MICt9MbU3/gyA3Kyv5Guf9knjulvYSmzoV/wFbJFJs2irjpG0wAvI5ltO0MdgmOeh0QnZcK8sa28TiY+Uo/1JaK/9SvUv5KTbPISNHgNLxBNNKDgTrdJlrdajBkxj5Y8pw5PvZaOGYk/wywJ82+h593mnBPvsGXpu8FkNpLcahOlPz+BXATWKfvVXuyf4sVXCv/Fp5SnKVHDj9xwS4wJuI32ZwXiu/B8GqzZp9Whabr+zv1uLbCz9SjXsUP0q919rOxyzKl/onLG2v5OFStczPoW9MO0MZ6LHFf8W2U1WnGaFa+DmKkp17M/qFDq/Fv1T/F/XlU0ihiTzdy91+BV37Xo6maMLxlX3w2mthgnLLh9S/BjHKYQvrSxaflYOSbMxcycWYVnNgFNMOp9OqcS8QI8vmvjymESvTkd4UabPu0If3//9g4WDJIS6BbcYE/3J52JxN1fg2K8G7DLumbaqidzkzW3dnjp8dBigM+vATcmvgLD9mxwEZLZhfp2uXRqKxo9KW9qIamCjk2uRbxJ8msi8O5ZAxYBu+zWBQxRLf3pDckxiulQXnkRNJtjGBYtasHwv1YECE2wwUPh6fi2OUYRW/xnWl7Lt2o2yeUcY/IOMwJdCCrEkcWKhryF6rA84rypMY2VPK4DpJfTCYqp5BTM73z/WQUO9vbnP2pdNajF29OVccBTrZX0w1kkPuQH821WTw2Fp5cF5RrkJH+odXaIfTNrFf+wLHTsFvU2qt/Hb+2rVrt1qfltHmSf3dWnxxXhf8SLEvZcZyCD9KvdfabgKzWF86cersoSQeztbiZWTU16t1cXOV/1qs1RVopNOUeNX9nBOiBPem9Js8thZ/nFes/4MvuIGwDFJ9MLta5I4AAAmPSURBVCE0pxpITqzT1Zt8rnfCqmthLcZeu7k2r43jSX1wIzut4lQsIRrpNCVeE+7hmueLFiZOkdoyJXExRtAGti3K1RidWSZTXx7bHMt1oTcFyaT7M1fXaSAf6qW/5XM+n0UnUwO+TcrhHeQ9CJO90LrdW/hdnMYEYLAKBhNfcHvNgnO/W3Oea5udZfDT23HdW9ob1xW7T3ywLH4OYPXlkhH1XNgGx4gZg8Jc+JC7NFXNKvug3mTbmP4515BjlfxbZUC70TbPJSPqob15UQZ1Rn4SB2KwWGpzrg4nS9I1PFdX6DjaoXMq3k5IBstzOrN/PfuMDvu03ew0NMgL2tXqL7mGDFEc82XIJbdre7WvMZnWyoPz1KdF+rG1GJuNbJ3D5s5uV+PzDLuc61z2TJUJ7Ub7UdadS07Uk9zPLemGOrvxu0uyhvKdHsn+K1Snn5fLhn6dqdtr7J9L7jVtp+oXKr9WD5xX1D+i/uGZxWTHPu/pntt+i7XTOflawHm6l4t4VsltU2ev4vcCW+Rey6nYNq+de8QJiVPCzT5rxWK5pRza3/39gNMhuf/bglsP5x5Nb+gzGXvAcQa6Afl0fBd5zf2YL5vJi/VivBVl+LLp7R3+YPkMFdlbMOxOJ0T3GU1n4Ob0VmC65PUdBTacluAlsEl605ALKdkmF5Lx9bS0OdpmcIAjVi5GhMZrcJySjv8cLc0bC/ZTHEH9B/CJ+soCZaskyDk8uKGxRz33o56cH7fo05w95WuqsPJdI8C9mR+Tzd/Z4Zq2WnKOOKN9+dJIwgGrXfivSHVoe9OnlZ8T9yaMBbswUPmhy7qPe5BHE8WaHvK4o3u5gCWAkzgewMfP2gunfJlD254+rfpXDvRjoI6xLqWMCHi27br/y6jyUNVe9YbcHNjBeC/jlF/4uCCP/pYB7ybXqS9LzDbk5UwV/O+zxesaZRns/j0p2B0jhMoIASFQHgFcwAxKcXL+98u3phZyIACb0dlwHn3Ozzw7jUmOtrbWAVl5k/gp5Gz9CeBWVXS+EBACQmAWAfnSWWjOMvbkv84ED+y09nPiXsA4HWft6VoQxzsmkifanjjliR3c7IR7fJZ5H88yN0FhlRmNwBG5GqP83vWG/Ax2vx4/1+M4A90MeHNKzeipS2MwK1HG6fEnZD0L2k+1hbJDsPvuVKaOCQEh0D0Cv0JCztvNkUlKHSMAG3G+XDoTBrrpSH7sWFwT7eudyGnyai0EhIAQWIOAfGkAtZ36r4BGZ1mt/Zy4d2aOvnd2ei2I4x3TaqecikW0Nfdsnn9+8au0EYGDc3UWnQPpzVk5GPwdJ7s+bC7scX43+7QFhOHAQX4xFJ0U7I6GSgWFQD8IuLfUHB38Qz9SSZIpBGArvknln8N+jHzOIfcPOmxzLlOnND0G2R5TAMhrf6bZVB41LgSEgBAohQD6OY74ki+dARj47Mp/zahxcbgHPyfuXZil6wN7uxbE8a7pNAi3N07FItoJ9zg9L/37N7Fyq9w8Akfl6rzGtzlH0Rt68D73pV2b1A7bHDBpo7nf8FjnyQZ4JsUnNI1J51aVeEIghAA6qug590P1KK8eAs5mbPAjOJ/uPq1z8jE4r/9mqEcLtSQEhEBDBORL48Dv3X/FaTE85PHeqQs/J+7FWq2vcr1fC04+cbwv2gSl6Z1TQeG9zF64Bzk0b7dnl5ybR+FqKiZ71xvyM2D8ARaOkGagmy+DOMjtEzz3L/5/I8o1S5CdI7r5x5QcPLiYUF7TmCyipAJCoH8E+KbuWf9iXp+E6GQ5fcnUZ0H8tI5OprvR3ZCXDo8j+RTovj7KSmMhcM0IyJd61t+j//LED2526OfEvaDF2mbu8VoQx9tyZqn1PXJqSSfL74x7/DM7jmDt7nnL8Op9fWSuhrA/qt54vuf/hv2M5ScsDG5bnOJ1CI9O8ngd8wVWUtI0JklwqbAQ6AsBdliQ6I27uehLOEnDKUtewTYMbHef6NghJKfF+ap7YSWgEBACQiAjAvKlF2Duyn9dSD9zoEc/J+7NGKufw7u6FsTxfogTkGRXnArocZbVG/fQt3LkKoN4eq45s1TSziG5GoHAtejNZ38GwLv70ty3EfqWYYpVHEsejKdgt4+ktoXAPhGgE3+EjmDqjwf2qdExpKbjeDnhQOhYmGyerNu9hr/gDgPy/DzoW8i7h7e7DdFS00JACBwUAfnSd4bdjf96J3J4q3M/J+6FzdcydzfXgjjekiZJbe+GU7Fadcy9J9DhOydfrDoq9w6Bw3H1nWrBrUPpDf7zGnjrXwfYZtyIC78u6z1xMN6qoPy93jWTfEJACIQRYDAVHdYXKPWCawUrw3hVzL1wHrAP58pisPvRRBC8omgXTfFTvyeQKelPHy5q0QEhIASEwE4RkC89M9ye/NeZ4IGdbv2cuBewWvusPV0L4nh7vsRIsCdOxejDMl1yD33rz3j2YsCbwbIL3Cm4UhCBC8w6fpYNKpKYeTS970P/02A22JD7nAaX8YjT8USMqhR3fKO8FzaJEUB/UBmDksoIgR0ggM6AQdT76LS6/oOBHUCZTUTYhPNL+Z/P0UYMKvc2qvvTnmTKZgBVJASEgBBIREC+9BawPfivWNNCFz4ode/nxL1Yi9Ytt4drQRyvy4mtre2BU7E69s49yMeBRgzGv49nnZtYvVTuFoEjcTXFpkfTG/rwpQ8T74eYnu/h2R9y/wVZ/4SsjwapI39w3vAHlQp2RwKmYkJACAgBISAEhIAQEAJCQAgIASEgBISAEBAC+0BgbcBsH9pJSiFwTARw3fJFFUegf5T6osqC3XePCY20EgJCQAgIASEgBISAEBACQkAICAEhIASEgBC4YgT4le3XLgB2xTBIdSGwDwRwrXIEOgPdX6UGun0NFez20dC2EBACQkAICAEhIASEgBAQAkJACAgBISAEhMDuEUCwjPMSf4uF05koCQEh0D8CDHT/imt309SvCnb3b2hJKASEgBAQAkJACAgBISAEhIAQEAJCQAgIASGQiACCZr/hlJcYMaqAdyJ2Ki4EaiKAa/Q7tMf/oUuap3tKRgt2P0Cl/mITl0+do2NCQAgIASEgBISAEBACQkAICAEhIASEgBAQAkKgewQseIa41+PuhZWAQuAKEcC1+TnUfoRr9YsU9XHefT+ejXMf8HwLdvMN1ytv+YGZSkJACAgBISAEhIAQEAJCQAgIASEgBISAEBACQmDPCCCIxvm7P3ZBtT2rItmFwKEQYMAaCnE092crFGP82o9nD19w/B/NOHrruCSKBAAAAABJRU5ErkJggg==\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
177
178
179
180
181
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
      ],
      "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": [
    "from lbmpy.maxwellian_equilibrium import get_moments_of_continuous_maxwellian_equilibrium\n",
    "\n",
    "eq_moments = get_moments_of_continuous_maxwellian_equilibrium(moments, order=2, dim=2, \n",
    "                                                              c_s_sq=sp.Rational(1, 3))\n",
    "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": [
251
       "<lbmpy.methods.momentbased.MomentBasedLbMethod at 0x7f4335401128>"
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
      ]
     },
     "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",
    "method = create_generic_mrt(d2q9, relaxation_info, compressible=False, force_model=force_model, cumulant=False)\n",
    "method"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Example of a update equation without simplifications"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
282
       "<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>"
283
284
285
286
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
      ],
      "text/plain": [
       "Equation Collection for d_0,d_1,d_2,d_3,d_4,d_5,d_6,d_7,d_8"
      ]
     },
     "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": [
313
       "<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>40.47 ms</td> <td>114</td> <td>67</td> <td>0</td>  <td>181</td> </tr></table>"
314
315
      ],
      "text/plain": [
316
       "<pystencils.simp.simplificationstrategy.SimplificationStrategy.create_simplification_report.<locals>.Report at 0x7f43353c50f0>"
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
      ]
     },
     "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": {},
   "outputs": [
    {
     "data": {
      "text/html": [
345
       "<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>expand</td><td>20.96 ms</td> <td>116</td> <td>227</td> <td>0</td>  <td>343</td> </tr><tr><td>replace_second_order_velocity_products</td><td>9.31 ms</td> <td>126</td> <td>235</td> <td>0</td>  <td>361</td> </tr><tr><td>expand</td><td>10.97 ms</td> <td>118</td> <td>227</td> <td>0</td>  <td>345</td> </tr><tr><td>factor_relaxation_rates</td><td>19.82 ms</td> <td>118</td> <td>184</td> <td>0</td>  <td>302</td> </tr><tr><td>replace_density_and_velocity</td><td>4.98 ms</td> <td>118</td> <td>184</td> <td>0</td>  <td>302</td> </tr><tr><td>replace_common_quadratic_and_constant_term</td><td>15.69 ms</td> <td>106</td> <td>148</td> <td>0</td>  <td>254</td> </tr><tr><td>factor_density_after_factoring_relaxation_times</td><td>15.23 ms</td> <td>106</td> <td>136</td> <td>0</td>  <td>242</td> </tr><tr><td>subexpression_substitution_in_main_assignments</td><td>14.23 ms</td> <td>102</td> <td>132</td> <td>0</td>  <td>234</td> </tr><tr><td>add_subexpressions_for_divisions</td><td>2.65 ms</td> <td>102</td> <td>132</td> <td>0</td>  <td>234</td> </tr><tr><td>cse_in_opposing_directions</td><td>9.09 ms</td> <td>102</td> <td>116</td> <td>0</td>  <td>218</td> </tr><tr><td>sympy_cse</td><td>27.07 ms</td> <td>86</td> <td>73</td> <td>0</td>  <td>159</td> </tr></table>"
346
347
      ],
      "text/plain": [
348
       "<pystencils.simp.simplificationstrategy.SimplificationStrategy.create_simplification_report.<locals>.Report at 0x7f43353c5668>"
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "simplification_strategy = create_simplification_strategy(method, cse_pdfs=True, cse_global=True)\n",
    "simplification_strategy.create_simplification_report(collision_rule)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Seeing the simplification in action"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<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\">cse_in_opposing_directions</h5> <div style=\"padding-left:20px;\">$$d_{7} \\leftarrow f_{7} + forceTerm_{7} + \\omega \\left(- f_{7} + \\rho \\left(- \\xi_{72} + \\xi_{73}\\right) + \\xi_{71}\\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_{72} + \\xi_{73}\\right) + \\xi_{71} + \\xi_{76}\\right)$$</div>"
      ],
      "text/plain": [
379
       "<pystencils.simp.simplificationstrategy.SimplificationStrategy.show_intermediate_results.<locals>.IntermediateResults at 0x7f43356a2c18>"
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
      ]
     },
     "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",
409
   "version": "3.6.8"
410
411
412
413
414
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}