demo_create_method_from_scratch.ipynb 65.4 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": {
Martin Bauer's avatar
Martin Bauer committed
33
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAoIAAAAVCAYAAADSDI/HAAAACXBIWXMAAA7EAAAOxAGVKw4bAAAI5UlEQVR4Ae2d0XHcNhCGaU0KUJwO5A6kuILIHUjqIFEHyuTp9JaxO7BTQcbpwE4FGasDuYM46kD5/xPBIXHgAiBAEBCAGYogAGIXH3bJJU8ku91ud4Plvl9OHx8fO9OC+hNTeQ5lOeiWgw7SXOSgXw46NEZm/5a4qLoc5i8HHRQP0zoH/XLQwcRGleWgXw46KB76OgfdctBB5zLezkG/HHQYMxnnbbqhfhL3HXVd9wOW69vb21dY7pA/SCi/QeHpQUU+BSe9jptoVAAfcmmM7NbRGMmMGh+ZD2sbo8bITkBu0WxI5sPaxkhmJPJBzPIOyyt08ReWly8QGb5F5hMKP5v6RfkFyl9j/atejzLuy/QvFnb6FmVfWRAzuchBGwarX7HmwJIlyNucjxosdGGw/hHLGfIPqlytUVYtI4z9BByUDf+I/Dduo3xy8ZMjI84f9BLnlm1CE2RY/blmPoqvbS5qZoSxF+lnrnorGwhZQ1ZxfpaSj2ILmeIxD/VZnc9SMXKV48IHbfbxnxgIotExJuVvrM/U5Kg1yr4g/zvW+8Crb8uyN8hHCwZ95PRtf8L6IAhSesdcQ04OfKjDH1gY2DDAofN8P8egUkY8OfEi5RLrfUKeDsADCe11chGUCyPo4TW3TyNb9rcfs5M/V8rHay4qZVSqn3npvczDnvYq1M9S8inVz5Iwgv14ybEdh1C/DwSPLEbNRu/1Ntj5F5QdYz3cfUP+AWXcPmiv7++6vUAOZVPnVGlTPhwkuWO5xHKNzT8dBl4dIzDhPP08ZgNevDtIm+UdVD3lwsh3bvVxOG2Dha8/V8WHEMHIdy6qYwRMRfrZAr1pEt6pVD/DQH3n1ZuN2qFUP0vIyHcunI5DtkDwChPzQU3SaM07K5Of1Pq6f7A+xz6M6mMkLzm9rtQ5lnzbGLbmY9PvoL5SRucA8Z/BLngnkBc0vMoaUkaMBp1WzpTqZytjWd59hTZEWKX6mZfey62iK9XPUvHxRpuRn6Vi5CXHlc9sIIgOKHDuJ17W8adIPfEOCxPrY6QlcqjzVQzhUh+Z8JFUlOpqY8SAj/8/quxTZ2O6cMiBka7nWtul+tlaPGL1W5MNkVmpfrZE7yU2UqqfpeKzhCn3ycHPUjFaIsfK5zuBPK9eKHSScDI1nTQnbbDxUi/w3Q6QQ53fYDHdyfRVQ2q/KR9JMYe6qhjBljhXpsT/p+xQb7q7vSkjk7JrlJXqZ2uwWKHPKmxIcSvVzxbqrYbttC7Zz1LwcYI432hzP0vFaKEcK5/ZO4JgzgcP7g3sVZBnurui7hK6BIuGridFS+VQ58lPfZNe421szSdkJNUzgkMxCKSdqCeJdZ5bM9L1WWu7VD9bi0fMfmuxoVlmpfqZg96zY56peFZ+tgKfGWxOxVn6WSpGDnKsfKRAkMGcCuycZmPUiO8mTJFMcqhzjEDUpn8JfObG0Bg9PSTyAU70bgZSCYxmVI9eXKqfRQfh2WGzoXL9jA+RSccHT1Nwal6Sn23BZw5irn6WipFNjpWPFAjyCubBQJ6dziV11cP3CoampXL4e3iKO4Jb8wnhWzUjBH98kuoz1nzSei5tzWhOr9jlpfpZbA5r9FeLDRnZlepnjnobxywUPhs/W4mPgM5alZ2fpWLkKMfKRwoEjfQhWAWHprtuqoyCg1KAnLkALUgf150D9HYVEaNdtYwwP3xVSoe1FASyyaaMqECKFGCvVfAJnINqGZXqZx56e5nGc/Gztfh4wTxsnJWfpWLkIcfKRwoEGcypwE5Hz38+NN11o0Am1sdIS+RQZ+nqK4Ze7CMHPkvHUiUjOM4FgPFTikMQiDw/xWOy5RwYLZ1f3/1K9TPfcaZuX5MNDWxL9TNPvYfxemSK9rMEfDxQTppm42epGHnKsfKRAkEGU6YTJGeAv0nzYQk9naHgDko+qArkqcTS5CxnJIDBqPWOZKBeFJcDn9GwvbLVMcJ8n4KQ6VOJDA5NFw45MHKe1EB7LtXPUvFxlqM1rMmG9kNfw88CbVtNiXS87lz1DtSlWD9LxEfNle86Cz9LxchVzgiilY8UCN6ho9ejzoYsFOGrWb5hzZPoPiHPgI/v7xu+4NCX8UW+/PScd8J+TnK0jhmMftLKJpuhevWdbc5nMii/jaoYYb55QcODMF8e/X60sOwa28OFywjj5oxGuohZ6E/fq87PRCijylA+o658s9XYEMGAc3Q/izh30vHaSe9QXbB/qeezJHx8nWvUfnM/c7X9CDbkNBcjNsxa+RxpO4w3GUydjgu0PDvnHZYbLsjze7f8zi8dbp+Q5wl2/xMq8jxZLUlWOVqnfGnn8Ok7rW6/GUmvXPh0GM9HLNTnt368X/oyzosp1caIAR8diP8fOF54IWMKAlG8fyn65nbkMrdoU62fJeTTucii4YxSFn7mojfaxLCh6H4WSS9OiXS8dtI7ki4lns9S8unAucTzWSpGTnJGxyBmrcehbrfbvcVy/vj42OkLyu+xnOrlvtvo4wLLse9+vu0h4wTLvet+oXpRFpZi+JAL9G2MDLY+tpnG6PBY0PjITMZ8mp+5sfLxM7QNPoegj6KO1z58epsLYlQan+Zndj+z2RDq9/GfdEeQkSQ/cDz8Yz0LFibeOZy787KwS+NufDkwdXZNoXqVxodcGiO7dTRGMqPGR+bD2sYoLqPQYzW1Ke143Wworg2xt1A7epY2JAaCCN74Pw1zT1Xapwgt0McxVjHeKyjKgxz+9EddqbM1xdCrl1UEHwKBvo2RxTIaIxlQ4yPzYW1jFJcReEY5h6Cfdj4TpqYkPs3PhInsq3yOQ2Ig2PfH77TyBbxL0xUUmvt6w9I+TftRx7lvyprax9KrFD5k0BiZLGFa1hhNeehbjY9O5HC7MTpkopf4MIp1rKYOpRyvffhwXLEYlcKHY26MSGE+OfN5wd+I0Y96+vcSQdvwsIfqv48sz7F2utum9ku1hl58AIBfirC+NmYNnXLnwzE3RvaZb4xkRo2PzIe1jVFjZCcgt2g2JPNhbWMkM7LxQT0fJOW//fHVMpf/A58tUf0oNZcSAAAAAElFTkSuQmCC\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": {
Martin Bauer's avatar
Martin Bauer committed
59
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAVoAAAAaCAYAAAD2fpSuAAAACXBIWXMAAA7EAAAOxAGVKw4bAAAHQUlEQVR4Ae2d73UUNxTF1z4uwFCC00GADkwHhlSA00E4fLK/5UAHMRUQ0wGkAjAdQAcx24Fzf8No0G5mxk8zWvmtVzpH1t95774r8UYraZe9m5ubRQ3DDJyfnx+q9VXb46hNX6h+OfxUbfHIgNex9IjLI6aUOZUb/1x5B+vgJfBEdUulH9fbdrT8Wlz8HmxX/i/lrxR/CXU13RoGvI6lR1weMaVMtNz4TfLkH+j3ch3oflyhDjjZp0qrk/1JzKn4OP5ZXLxW/kh1v0Z1NbsdDHgdS4+4PGJKmWW58VvlvZNv+LAOtHO0amw+IivtVm/rnXe0DB+fd9T2+2a217H0iMsjppT5mBu/SZ785xeBvFT6Rwx2L+zRqgEvTIeLuEPNrzIgfljRniitWwer1GxdyetYesTlEVPKhMuN/zZ5av8qfI+UNmc5zR6tChzyPFb6dAy82vm4fBkLGOt/39pa+9leeXTfbNs1e7yOpUdcHjGlzNfc+I3yOMt5q/gMrOEwjM3b3pWshLKlwAPXio8Vw8m7srsTxAN2s5rt3lK7Y/39stTrWHrE5RFTymzMjT9BHv6Ug7FDxWVwtKeq7F3N0kltjVdWnn2HnTsEkt2Nk1XacKQUDriZ8U1pDVvEgNex9IjLI6aUqZYbf4o89cU/sF/7XPHiQIXgOOuBT88otuTyMeBlxBUb4/+7wtHzeK1yxIDXsfSIyyOmlKmUG/9EefhUFqkXrGiPFb9JULNpq3wNqwxwZ5btE9IuiC+cbQ3bxYDXsfSIyyOmlNmWG/8UeRyIsaJd7J2dnXG4xb3QWw941IetA/YpH+R2zJKHM2MLg9P8ryq/UdqEtu2t0mYLI9RPTUvqsmKMMPHIE8UXimxZ/KZI+KQ+739k5/9N0Rf13fjYYFmkj+IoF1HfItgAZAkRLrqP2mCRl6NPbkyRvCLcR/qycBrJ2wh+yWcRy22uB/v6g4PjoOuuwysBw7kCLHzlNWDircBpf65QUpcVMxvnb1oOPukhDiCPVWaLgjHiBZczpOgrzZdnbNYxSLHBKnNuv9yYPM8LC1ebxh/86hGO9qHi0oJqU33kTFi54VwIHDitHzJRl+XbaiV1YYwlCBMr+XVHyoslrGAZo/V2i+jePin6SvPlGVsvmT2VKTb0PL6RqtyYPM8LC4GF8Ae/+pA9Wg8r2msZHpwKq9c/18hiCb5et9bFXCypywrqs+yPXy68WL6EOqVZtkwiMCn6SvPlGVtE4Wg2xYZRQRkbc2PyPC8stJXAH1a0C1a0eF1WTHcW5Egaz68Uh4rj7+70qo5bEdTlWtEW0yXMpiAbuQYSB3h4F1fkzKfoU9+ifHnGZh2DFBusMuf2y43J87ywcFUIf/Cr1zhaT4GVGyu55h93CwynE+6k5cRaUpcZt2zHXkJY4S9Ux6VnXjbZQ4K+4nx5xmYdiAQbrCJn98uMyfO8sHC1Sfzh3+wSR8tH1lBhATbYRwM4Vw57tfFHaHT17s+W1DVocNQwFQ/PKX5QDA62uTamcswDm/bLSN1C5Ulc85xisj7p3vjYlMAWcziUn8ot8mbYMASnq5+Ka5OYBM7zvOi4G8mY8Y/IGGpaWdFyCMVXa2eFdhJ8V8p9s6khdi4LycL5ELmJ0IWSujqlI5mZeIKN7BnhPNnX6Zxqy0E4KGxQlNbXml5ibJK5SMHW9h1NZnKL7Kk2bBLXRjC1gD3Pi1FOU/BbBPX0wYnzaXzJPVoGobnrRUVP54XqL1WPE8Ahk0Iu+4rc7Yzvu3JBlzDp9wAkC9lca0L+v4rcb2vu1qptfUCL6RKGW4PwTcLT2syNguZ5lbnixY0DVvK8tOJNexV/hDvQt/GxmcGFGVvgbyydyi0yp9owhie0TcW1YUxm7rcdfxgHayp7+UbpQ6XPFvxMopztdxwu+blRck4UD+fKaXG9lqyrIVkldQ1hiOtz4onlDuVL64txSHexsYn1WvK3YTPKyDaPLfqsfe5yzC0Yb+N+2/FbOAh9ZOuV4inl/dY7/62UFVSO8EQefJkqCO+v2G07KM+bktXs2LWukrosJk3CYxE80KeIvpJjM2DnYPVEbIPyooYi3Eb6rFk3uCZyv+34TePU+i9uTOFbO0fLEhenNiu0wvnIPyWwhRFfaWIL4UIyu9P3WGhJXbHeofxMPENiB+sL6ys5NoM2DzQkYRuQsVJdmNsV3WMFh7iSuN92/GNj09PG9wHey+Zm0XlABxW4UvVRkf85oNex9Qjqq3qu57s9274OI3V81ZTfXOD3FNibZYU7dne2pK4R2F3THDydkIRMSX0lxyaBgqZrKjaL/JLcWvCEPt5wpXK/7fjDOFhSuOl2CeL/yoaP6v/Iud364zIWLbVPZaAyUBnYRQbkQ9kd4Bplt+jsHC2EqKE57Vba3OXcRZKqzZWBykBlYCoD8p1c6eLTeLeaRVY4DGvkqpFtg/gye1Nf/1QGKgOVgcqAiQEWqXzbbCX8B9yIT2mRLUjhAAAAAElFTkSuQmCC\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": {
Martin Bauer's avatar
Martin Bauer committed
85
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXMAAAF5CAYAAABtIcr0AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nO3deXRlVYGo8W8nNc9AioJiLkAUlEFaRZBB1MaIrx1BbURQVJwQRX1Ka7e+9xwaFWgFRcSxFXyNc2sThxYUROQpKs4IMoNVRYoaoaYk+/1xc6uSW7nJHc/4/dZiQW7OSbauU1+du8++54QYI5KkfOtJewCSpPYZc0kqgGm1L4SBwQA8EzgeWAiEOvtGYB1wPfDD2N830q1BSo0IA4PTgX7gaGAuHrvKiTAwOAs4GXgKMJvJj901wHWxv+/asd/YIebAPwOnNTGOFwNXj+4npekS4OlNbP9i4KvAe7ozHGlqYWCwF7gCeHITu50aBgY/F/v7Lqi+MG6aJQwMLgZe1sJ4TgkDg7u3sJ/UEWFg8FCaC3nVKWFgcGmnxyM14ak0F/KqM8LA4KLqF7Vz5k+Y4LVGBOCwFvaTOuXwNvb12FWaDm1xv17g8dUvasM9s+XhtLev1K4Zbezrsas0tXP8bTvuJ5oz39F5z96LO343i97eytc77TrE5355VxsDkLrvq5cs4rqvLuT+O2bw1Oes5/zPLE97SFJTfnDVfP7j4l1YtXw6C3cZ4tx/W84TT9g40aaNxRzgrPeu5HmvXduxQUrdtsvuQ5x67ipuuW4uWzbVWx0gZdPPvzeHL31oMe/41IM8/qhNPPTgpL1uPOZS3pz44g0A/OU3s1j1N4915ctVH+njlDev4tBjNgGwZK+hyTZv/GLnlR/p49QD9ufNz9ybX/5odnujlCTVNTwEd/1xFmtX9XLmEftx2iHLuPjNu7Lp0brvMBuL+Sv/5SE+f8udfPn3d3LSaWv4wKv25L7bp3ds4JKk7VYt72V4CG66Zj4f+e69XHrdPdz1x1l88YO71NulsZg/4ehNzF0QmTErcvIr1/GYwzfy84G5HRu4JGm7mbMrd0A8+ZWrWbzHMDvtOszzz36YX/+4bndbvDdLiHi3RUnqjoW7jLDTrkOExq/bTx3zdQ/3cNM1c9i8MTC0Fb73pfncdsscnvysR9oZq9R1Q1th88bAyDCMDLPtGJby4OkvWst3P7+IVct7Wbuqh/+8YieOPHFDvc2nvsI/tDXwpQv6+PDrZ9LTE9l9vy2864oH2Pdg/1Qo2774gV34+ie2zzHe+N0FvOiNqzjrfatSHJXUmDPes4p1q3t57VP3Y/qMyFH963nF+Q/X23zqmO+8ZJhP/uTejg5SSsJZ7zPcyq/pM+C8S1Zy3iUrG9nc+5lLUgHUxrydq5peEVVeeewqTR05/mpjPuFn/hv0aDsDkdrUzgV5j12lqSPdrY35r4BWLmwOA79sY0BSu25ucb9h4BedHIjUpFaP3Y3ArdUvxsU89vetBz7awg+9KPb3rWlxQFLbYn/fncAXW9j1Yo9dpSn29/0G+HqzuwEXxP6+bWf1IU7w4Z8wMLiMJp4BGvv77mhyIFJXhIHBQ2jiGaAeu8qKMDB4GHAUUz8DdDXw49jfd8+4/SeKuSQpX1yaqMILIXw9hHBJ2uOQuskzcxVaCGF34C4qFzp3jTF6GwoVkmfmKrpXU5lnHAFOTXksUtd4Zq7CCiH0AMuBxaMv/THGeEiKQ5K6JrEz8xDCziGE74cQzgghzErq96rUngWMPdb2DSEcltZgVB4hhHkhhDeGEL6TVO+SnGbZCBwDXA6sDCFcEELYI8Hfr/J5KzB/zNczgXNSGotKIISwfwjhk8AK4GLgicCWJH53YjGPMW4EPkVl/eR84FzgjtG/uY4JoYm7sEtTGL3weULNy73Ay0IIPiVLHRMq/j6EcB3weyrXaeZQ+TT9BTHGkSTGkfQF0I9TuRAFlbOkWcDJwPeB25yCUQdVL3zW8kKoOqI6lQLcS+UTnCdQaVr1+cgB+EJi40n6AmgIYQA4iYk/4bSByh/Ay4CPxxgfSHJsKoYJLnzW8kKoWhZC2B94G3AGlV5N9E5vCPhcjPHspMaVxtLED1H/DnfzGD8F812nYNSC2guftbwQqqZMMpVSb8puK3BhUuODdM7MA3A7sH8Dm0cqt3h8AHhajPGhbo5NxRBC+B6Vd3/1DANfiDG+OqEhKcdCCAcA1wGLqJxwNuKnMcZjuzeqHSV+Zh4rf3t8kMqUylQClZvO7IwfcFID6lz4rOWFUDUjULnpYKMh30ClcYlKK5BfaXC7EWAV8JQY44oujkfFUe/CZy0vhKohMcbbgWOp3GmzEeuoLOpIVCoxH12meDmTr7+shvyoGOOdiQxMuTZ64fMcJp8vr5oHvL27I1JRxBhvBY5j6qA/SoLLEcdK7eP8IYS9gduY+A/eCLAGeJIhV6NCCH8P/BcwrcFdRoAjYoy/7d6oVCSjF85/Tv0Tho3AbjHGRs/iO6bRg77jYoz3jl4ZfjbjlymOUHnHsDOwOY2xKbf+BFw9wev/OPrvq2peH6ayRlhq1BD1Qz4EfCmNkEPKN9oKIRxH5UyqemGhOrVyHJU/mAB7ut5c7QghRGBFjHG3tMei/AohHEJlWSJUPqb/Y2DBmE02AYfFGP+S8NCA9FeI3EDlHgYwfo78z2z/2+9+7+EiKU01Ie+NMf6aHefQf5lWyCHlmI9ZpriVmoudMcbNGHRJKZsg5COww0XRzaSwHHGstM/MobJM8fNMsGrFoEtKU72QV40J+qdJYTniWLl4OEUIYSaV+ShwDl1Ncs5crZgq5FmThTPzKXmGLilJeQs55CTmYNAlJSOPIYccxRwMuqTuymvIIWcxB4MuqTvyHHLIYczBoEvqrLyHHHIaczDokjqjCCGHHMccDLqk9hQl5JDzmINBl9SaIoUcChBzMOiSmlO0kENBYg4GXVJjihhyKFDMwaBLmlxRQw4FizkYdEkTK3LIoYAxB4MuabyihxwKGnMw6JIqyhByKHDMwaBLZVeWkEPBYw4GXSqrMoUcShBzMOhS2ZQt5FCSmINBl8qijCGHEsUcDLpUdGUNOZQs5mDQpaIqc8ihhDEHgy4VTdlDDiWNORh0qSgMeUVpYw4GXco7Q75dqWMOBl3KK0M+XuljDgZdyhtDviNjPsqgS/lgyCdmzMcw6FK2GfL6jHkNgy5lkyGfnDGfgEGXssWQT82Y12HQpWww5I0x5pMw6FK6DHnjjPkUDLqUDkPeHGPeAIMuJcuQN8+YN8igS8kw5K0x5k0w6FJ3GfLWGfMmGXSpOwx5e4x5Cwy61FmGvH3GvEUGXeoMQ94ZxrwNBl1qjyHvHGPeJoMutcaQd5Yx7wCDLjXHkHeeMe8Qgy41xpB3hzHvIIMuTc6Qd48x7zCDLk3MkHeXMe8Cgy6NZ8i7z5h3iUGXKgx5Mox5Fxl0lZ0hT44x7zKDrrIy5Mky5gkw6CobQ548Y54Qg66yMOTpMOYJMugqOkOeHmOeMIOuojLk6TLmKTDoKhpDnj5jnhKDrqIw5NlgzFNk0JV3hjw7jHnKDLryypBnizHPAIOuvDHk2WPMM8KgKy8MeTYZ8wwx6Mo6Q55dxjxjDLqyypBnmzHPIIOurDHk2WfMM8qgKysMeT4Y8wwz6EqbIc8PY55xBl1pMeT5YsxzwKAraYY8f4x5Thh0JcWQ55MxzxGDrm4z5PllzHPGoKtbDHm+GfMcMujqNEOef8Y8pwy6OsWQF4MxzzGDrnYZ8uIw5jln0NUqQ14sxrwADLqaZciLx5gXhEFXowx5MRnzAjHomoohLy5jXjAGXfUY8mIz5gVk0FXLkBefMS8og64qQ14OxrzADLoMeXkY84Iz6OVlyMvFmJeAQS8fQ14+xrwkDHp5GPJyMuYlYtCLz5CXlzEvGYNeXIa83Ix5CRn04jHkMuYlZdCLw5ALjHmpGfT8M+SqMuYlZ9Dzy5BrLGMug55Dhly1jLkAg54nhlwTMebaxqBnnyFXPcZc4xj07DLkmowx1w4MevYYck3FmGtCBj07DLkaYcxVl0FPnyFXo4y5JmXQ02PI1QxjrikZ9OQZcjXLmKshBj05hlytMOZqmEHvPkOuVhlzNcWgd48hVzuMuZpm0DvPkKtdxlwtMeidY8jVCcZcLTPo7TPk6hRjrrYY9NYZcnWSMVfbDHrzDLk6zZirIwx64wy5usGYq2MM+tQMubrFmKujDHp9hlzdZMzVcQZ9R4Zc3WbM1RUGfTtDriQYc3WNQTfkSo4xV1eVOeiGXEky5uq6MgbdkCtpxlyJKFPQDbnSYMyVmDIE3ZArLcZciSpy0A250mTMlbgiBt2QK23GXKkoUtANubLAmCs1RQi6IVdWGHOlKs9BN+TKEmOu1OUx6IZcWWPMlQl5CrohVxYZc2VGHoJuyJVVxlyZkuWgG3JlmTFX5mQx6IZcWWfMlUlZCrohVx4Yc2VWFoJuyJUXxlyZlmbQDbnyxJgr89IIuiFX3hhz5UKSQTfkyiNjrtxIIuiGXHllzJUr3Qy6IVeeTat9IQwMTgdeCDwNWACEOvtGYD1wI/CN2N+3uVuDlMaKMW4OIcwCNlEJ+p4xxgfqHrtvuxSmzVgYBgb/vfojgNXAdcB3Yn/fiCFXmsLA4Gwqx+7RwDzqd3cEWAX8d+zvGxj7jR1iDnwcOLGJcTxrdPvXNLGP1JaJgs41D70beMYOGx90JPT0TgeeUvOdfuCJIYSrMeRKSRgYDMCngKOa2O25YWDwkNjf99HqC+OmWcLA4IE0F/Kq48LA4MEt7Ce1bNyUy9Jl9zO09aSmf8iWTS9n3iJDrjQdQXMhr3pFGBicU/2ids788W0M6Alt7DtOCKHeWwxpnG1B3+exsPyeZQxtnejd5sS2bJ7ByvuXsfdBYMjVhA43qtXuzgQOrH5RG/MZLQ+nvX0JFU8LIXwNWBVC2K+dn6fyiDFu5m2feDnQeNC3bJ7Byvv2BeBD33imIVejQghHAqtDCF8IITyxAz+yI91tbDXL+8/cjZc9dn9euM8BvPLI/fj2pxe28cvHCSHsHEJ4C3A3MEDlIsB0YOdO/Q6VwJx5Q+yx7HZgfNDXDvbwyXfCG47r5fQnLOP7X54/LuR77P8Xps+IaQ1bubQb0Au8HLghhPDnEMJrQgjzu/Lb7vnzdP5hjwN5/5m7TbZZY29JX3rew+z9mBXMmBW56w8zOP+Fe3Hg4Zs4+MktrWAZfYtyDPAW4GQqV2jnjNlkuJWfq5ILPZE9lt3OA3ceyPJ7lrHbPnfysfMWM206XPS9YdYO/o33n7EnO+0aWLqsEnJn9NSaYSpBnwMcBFwEfGz0YvrHY4y/6thv+sQ7lrDskE1TbdbYmfkBh25hxqzK2UsIkRDggTubfmtQ5yx8FuNDLrWuGnSAe/68jF/8cD7Pey3MmgOHHDXMoccGfv49Q65OmwfMZsez9Xlt/dQfXDWfOQuGecLRj061aeMfGrronF15/p4H8obj92PR4iGOOXlDI7vVzIU/AHwA2JvJ11JKrasGfcW90NMDS/YGIqy8b1/2OgCW3/OIIVeX1J6trxydWz+i6Z+0YU0PX7mwj9d98KFGNm/8yv95l6zk3ItX8tsbZ3PrDbOZPnPKecbRs/C3Upn/nkvj8e4BjgshLG54fCq3d1x2GEc9e/w7vN5pDzF7XuUYGhnpBWBh3xo2bpjNxg3btx349yeF57y3rQv4KpWjaKxl1bPylwOnhBDuAy6MMV7R0G/57P/q48RT17LbPkONbN54zAF6p8ERx2/kR1cv4JuXLeLUc9dMscfFVObDm71twHwqf6tJjbn5+3Dg4eNfGx6CjY+Mf23tqkXMnA2rlu+57bVf/PBfExihyqt6tn4gcDkwdcxvu2Umv/vZHD55/d2N/pLmYl41MgR/u7uRM5nHAm8Czhz9utH5o7XAM2KMt7QwOpVQGBh8CfC/x704f6eZjAzvw4r7YMleldf+dtcG9nncFvY8YHDbdh/65utjf9+1CQ5XORZCOBm4Emh0Vd96KhdML6MS86n9+vo5DD44nVccuj8Amzf2MDICrz92JpfdcM9Eu0x9xrxqeS8/uGo+j64PDA/BTdfM4WfXLODwY6eckI8x3hZjPAdYDLwW+BWwEdja0P8gqVVbNs9g/ep9OOIE+M8rYOuWYe64FW65bh7POPWRKfeX2rOFyq0mrgdOBxbHGP8pxjhhiHfwvNeu4TM338ml193NpdfdzTNfsobDj32ED3zt/nq7TH1mHgIMfHERl797CXEEdtl9iDPfs5LjX9jQBVCAGOMm4CvAV0IIB9Ha2brUmLHryM/7+B1c8LoDeOtJvcxbOMRp75jGnPl7MbT1TqZNb2guUmrCuLPwhuNda/bcyOy525doz5o7wvSZI+y8pO6y7aljvvOSYS7+/n0tDWgCMcbbgHNCCO8AXgC8HXjc6Fimd+r3qKRqPxAUArzxAujpHWbpfncSR8K4degGXe3bQuXa4C+AC4H/ijF29rg6632rptoktfuZxxg3xRi/EmM8ksqNZi4HNoz+M2vSnaWJTBTyWmPXoTd7Lxdpu1lUzsLXUFms8dgY43Exxm93POQNqo15Ox9rbnnfCebWv03lg0VSY3574x5ThryqNuh/uLkvgRGqOH4HfIftc+HntzydUtGR7tbGvJ0LQw3Podcz5mz9JTHGKd9WSDD6hKCvXXoh0PgnO8cG/Yp//o8kHhKtYogx3htjPKWDZ+HtdHfbvrUxv5nK3E+zInBTGwOSWrLtCUG33QJ77P+npj7ZGXoiS5f9hjv/AF1+SLQ0iZ+1uN8gcFv1i3Exj/19g8CHW/ihF8X+vhUtDkhqybhHva1f3UsIH2nyRwzT0/NehrZ09SHR0mRif9+9wKVN7rYV+JfY37ft5DvEuON0TRgYXMr25yhOZh1wY+zve6DJgUhtqffMztFj9xjGfqDjY2+9gmnT1/HGD79tzI94GLh+9ASGEMJMKuuCAfaMMXpMK1FhYHBvtj8DtJ5I5Yz8+tjft3rc/hPFXMqyZh++HEKIwIoY46T3gzboyrPUliZKrWg25M0Y90xRp1yUM8ZcudHNkFcZdOWVMVcuJBHyKoOuPDLmyrwkQ15l0JU3xlyZlkbIqwy68sSYK7PSDHmVQVdeGHNlUhZCXmXQlQfGXJmTpZBXGXRlnTFXpmQx5FUGXVlmzJUZWQ55lUFXVhlzZUIeQl5l0JVFxlypy1PIqwy6ssaYK1V5DHmVQVeWGHOlJs8hrzLoygpjrlQUIeRVBl1ZYMyVuCKFvMqgK23GXIkqYsirDLrSZMyVmCKHvMqgKy3GXIkoQ8irDLrSYMzVdWUKeZVBV9KMubqqjCGvMuhKkjFX15Q55FUGXUkx5uoKQ76dQVcSjLk6zpDvyKCr24y5OsqQ12fQ1U3GXB1jyKdm0NUtxlwdYcgbZ9DVDcZcbTPkzTPo6jRjrrYY8tYZdHWSMVfLDHn7DLo6xZirJYa8cwy6OsGYq2mGvPMMutplzNUUQ949Bl3tMOZqmCHvPoOuVhlzNcSQJ8egqxXGXFMy5Mkz6GqWMdekDHl6DLqaYcxVlyFPn0FXo4y5JmTIs8OgqxHGXDsw5Nlj0DUVY65xDHl2GXRNxphrG0OefQZd9RhzAYY8Twy6JmLMZchzyKCrljEvOUOeXwZdYxnzEjPk+WfQVWXMS8qQF4dBFxjzUjLkxWPQZcxLxpAXl0EvN2NeIoa8+Ax6eRnzkjDk5WHQy8mYl4AhLx+DXj7GvOAMeXkZ9HIx5gVmyGXQy8OYF5QhV5VBLwdjXkCGXLUMevEZ84Ix5KrHoBebMS8QQ66pGPTiMuYFYcjVKINeTMa8AAy5mmXQi8eY55whV6sMerEY8xwz5GqXQS8OY55ThlydYtCLwZjnkCFXpxn0/DPmOWPI1S0GPd+MeY4YcnWbQc8vY54ThlxJMej5ZMxzwJAraQY9f4x5xhlypcWg54sxzzBDrrQZ9Pww5hllyJUVBj0fjHkGGXJljUHPPmOeMYZcWWXQs82YZ4ghV9YZ9Owy5hlhyJUXBj2bjHkGGHLljUHPHmOeMkOuvDLo2WLMU2TIlXcGPTuMeUoMuYrCoGeDMU+BIVfRGPT0GfOEGXIVlUFPlzFPkCFX0Rn09BjzhBhylYVBT4cxT4AhV9kY9OQZ8y4z5Corg54sY95FhlxlZ9CTY8y7xJBLFQY9Gca8Cwy5NJ5B7z5j3mGGXJqYQe8uY95BhlyanEHvHmPeIYZcaoxB7w5j3gGGXGqOQe88Y94mQy61xqB3ljFvgyGX2mPQO8eYt8iQS51h0DvDmLfAkEudZdDbZ8ybZMil7jDo7THmTTDkUncZ9NYZ8wYZcikZBr01xrwBhlxKlkFvnjGfgiGX0mHQm2PMJ2HIpXQZ9MYZ8zoMuZQNBr0xxnwChlzKFoM+NWNew5BL2WTQJ2fMxzDkUrYZ9PqM+ShDLuWDQZ+YMceQS3lj0HdU+pgbcimfDPp4pY65IZfyzaBvV9qYG3KpGAx6RSljbsilYjHoJYy5IZeKqexBL1XMDblUbGUOemlibsilcihr0EsRc0MulUsZg174mBtyqZzKFvRCx9yQS+VWpqAXNuaGXBKUJ+iFjLkhlzRWGYJeuJgbckkTKXrQCxVzQy5pMkUOemFibsglNaKoQS9EzA25pGYUMei5j7khl9SKogU91zE35JLaUaSg5zbmhlxSJxQl6LmMuSGX1ElFCHruYm7IJXVD3oOeq5gbckndlOeg5ybmhlxSEvIa9FzE3JBLSlIeg556zEMIu4cQ7gshvKDO9w25pMQ1GvQQwtkhhL+GEOYlN7odpR5z4E3AEuDKEMLzx37DkEtK01RBDyGcDVwELAVOT3h444QYY3q/PISZwEpgwehLG4F/jDF+y5CrU0IIEVgRY9wt7bEon0ZbtWn0yz1jjA+MCfmc0dfvBfaNKUU17TPzU2rGMBu4KoTwPzHkkjJigjP0dzM+5AA7AycmPbaq1M7MQwgB+BNw0CSbGXI1LIRwNPAdINR8a6fRf6+ueX0IODbGeFu3x6ZiqDlDn8i1McZnJDWesdI8M38KsOck398I/ENCY1ExrAXmUYn32H+qal9fBKxJeIzKtzOBRyf5/tEhhP0SGss4acb8XYx/i1KrOuXy/Em2kbaJMf4B+HMTu/x3jHFFt8ajYplgjnwivcBbkxnReKnEPISwFDiJHd8O1zLoatZHgQ0NbLeByh9MaUoNhhxgOvCqEMLc7o9qvLTOzN/I1CGvmg18PYSwfxfHo+L4Go0dW48A13Z5LCqAEMJRwGVMHfKxEl+mmHjMRy8gvAmY2cDmG4D1VM62HuzmuFQMMcaNwJeoXNysZyPwcS+uq0G3AZ+iMlf+SAPbzwXOH13kkZg0zsxrlyPWilT+D7sdOAfYNcb4ztE/pFIjLgW2TvL9HuCzCY1FORdjXB1jfAOwG/BO4AGmnspLfJlioksTp1iOuJlKyH8EfAj4WVqL75V/IYRbgUPrfHsgxvicJMej4ggh9ADPBP4JeDIwjcpcea1ElykmfWY+0XLE6lTKx4ADY4zPjTHeaMjVpnoXQr3wqbbEGEdijD+IMZ5A5YThM0w8BZPoMsWkz8y/xfa1449SmQf/IPB/Y4yTLcSXmhJCmA08RGX+cqwVwFLny9VJIYT5wBlUllwvpPJ5h63Ap2KMb05iDImdmYcQdgP+BzAMXENlaeJBMcYvGHJ1Wp0LoV74VFfEGNfHGC8F9gZeBPyEyprzV4cQmlkF07LEzsxDCL3Aa4DvxhjvT+SXqtRGb9b2CyrLW6FyXWYfPyikJIQQDgSOBT6fxLRxqndNlLqt5kKoFz5VWGnfNVHqtuqFUC98qtA8M1ehjV4IHaSyYsoLnyosYy5JBeA0iyQVwLTaF8LA4GzgNOBpbH+cWz3rgZ8CV8b+vsnu8St1XRgYnE/l2D2ayjrfyawBbgCuiv19m7s9NmkyYWBwMfBy4O/YvvqqngmP3R1iTuXuYE9tYhxHAcePDkRKRRgYrN5v5bAmdjuGyqeSX9eVQUkNCAOD84ArgX2a2O0YKu09u/rCuGmWMDB4MM2FvOpJYWCw3n0wpCT8Hc2FvOrpYWAwlSfDSKOeRXMhrzohDAwuq35RO2f+2DYG9Lg29pXa1c7x57GrNHXk2K2N+UR3/mrUjDb2ldrlsau8aufY3bbvRHPm471g7wPHfb1lc+BZL13DWz62so0BSN23ZVPg4jfvyu9vmsuGdb0s2WsLp79rkGOe28gDBqR0PXDnNC552xLuuHU206ZHnnLSes65cCXTJm7/1DH/5r23b/vvRzcETjv4AI57/vqODVjqlqEh6Fs6xL9++15232eIG787l4++cSn7HnwXeyyb7ElEUvouedsSFu4yzJV/+CvrV/dw/gv34hufXMSp566ZaPPm1plf99X5zN95iCOO96k/yr458yJnvW8VeywboqcXjn3eIyxeuoXbbpmV9tCkKT10/3SOe956Zs6O9C0d5vDjHuHe2+o+brO5mF/71QUc//x1BD9rpBxa9bdelt87g/0O3pL2UKQpnfyq1fzkm/PZ+EhgxX3T+M31cznyGXWnCBuv8t/unsafb5nDSaev68hApSRt3QIfes3uHPe8dex3iDFX9h1+7Ebuu30mp+x/IK984jKWPX4TJ7yg7rNHG4/597+8gMccvpE995/sQblS9owMwwfP2p1p0yPn/pv3Mlf2jQzDv7x0T4569nq+cfftXPXHO9iwtofLzl9cb5fGY/6Tby7kxFPWdmSgUlLiCHz4dbuxdnAa77vyQaa7ClE5sHZVLw+vmMYL37CGGbMiixaP8KyXruPXP6l9DOI2jcX81htmsXrlNE481VUsypcLz1nC/X+dwfuvvp9Zc7xFqPJhp12H6dtjK9+6fBFDW2Hdwz386OoF7H1Q3fsITb00EeCHX1nIk565nrkL/MOg/Hjwrmlce/VCps2InHbIAdteP/sDy3n26Z6YKNve/dkHufzdu/LtT+9MT0/kcU/eyBsuqPv5nsZi/vZPOs+o/Fm63xDXPPSXtIchteSgIzdz0ffua3Rz17OX18EAAADhSURBVBhKUgHUxrydaRSnYJRXHrtKU0e6WxvzduYRXX+uNLVz/Dl/rjS1c/xt27c25jcBwy38wBHgZ20MSGrXT1vcbwtwcycHIjXphhb32wL8v+oX42Ie+/vWAO+luaAPA/8n9vcNtjggqW2xv+9B4AKaO3aHgHfH/j7voqjUxP6+XwJfbHK3IeA9sb9v2ydCQ4w7TteEgcFFVB5JtAAI9cZA5RT/57G/b3WTA5G6IgwM7kLl2J3L5MfuWuCm2N/n9KAyIQwMLqXyxKxZtHDsThhzSVK+uDRRkgrAmEtSAfx/3Sy8hI22nQYAAAAASUVORK5CYII=\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": {
Martin Bauer's avatar
Martin Bauer committed
109
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWQAAADhCAYAAAD/Ec//AAAACXBIWXMAAA7EAAAOxAGVKw4bAAAbXElEQVR4Ae2d0W0cOZeFpYWfDWMMOAApA2sdwVgZ2M7AOxkM4Keet4Gcge0IBusQZiIQrAxmAzCwguEI9pze5v+XytVFqrtYvBx9BEqs4q3iPfzYfcVmV7NON5vN85OTky/aptLn33777fWUgTIIQAACELgfAcXTv3XF2dRVsp0+Ghjea98nD9P/DA/YhwAEIACBowhcTVx9qbJXLh8G5A+K0ATgCVoUQQACEFiCgGLsx3E9KnPRDwF5fN7ksS72FMd/a7vQ/rfJkyoVtvSda1Jkbb1rj8wWbblX17Q9MrdpxXdLa+kfjpDvehwcyfkTHX7SdqvtP7VNzoGofPHU0neuMZG19a49Mlu05V5d0/bI3KYV3y1dQ39pQPZIePvlnkT9qn2PkldJ8tfMd66BkbX1rj0yW7TlXl3T9sjcphXfLV1D/3/cdckRBCAAAQi0IkBAbkUevxCAAARGBAjIIyAcQgACEGhFgIDcijx+IQABCIwIEJBHQDiEAAQg0IoAAbkVefxCAAIQGBEgII+AcAgBCECgFQECcivy+IUABCAwIkBAHgHhEAIQgEArAgTkVuTxCwEIQGBE4HSwHvK5fhq4d7U32bygkNe08FoWzn3ujbZr2bx0Z7XU0neuUZG19a49Mlu05V5d0/bI3KYV3y2toV91/pe8eLXN0+KAfFcWRxCAAAQgsASBYUBmymIJotQBAQhAYAECBOQFIFIFBCAAgSUIEJCXoEgdEIAABBYgQEBeACJVQAACEFiCAAF5CYrUAQEIQGABAgTkBSBSBQQgAIElCBCQl6BIHRCAAAQWIEBAXgAiVUAAAhBYggABeQmK1AEBCEBgAQIE5AUgUgUEIACBJQg8uk8l+onf1e78/1V+ru1KZXvXv7hP3SXnytdznec1NS60/63kmijnRNUuXU37NNc/HegL+5qM+prL9bntD1V7cUAWoC/i9LvyzztgXmDoi44vtVULyqrbfj5pu9XmhY3OtHWRomuXviZ9Wtp5UfVF7tfI2nL9jvaTk6KALFBejeiJ8m0wNljtf9sdf9DhpctqJPtRva9dt/Z/VeYRSRcpsnZpa9anJZ0XWV/wfuX9UvICW/icpV4TpXPIDoheanOcrlXwUmI8iiX1RSB6n0bX11dvo7YLAqUB+aVa4ymDcfJ/YyfbSX0RiN6n0fX11duo7YJANiAXjn5/6qK1iNwSiN6n0fXxMoJALQLZgCzHKdim0fBQSxo1M2UxpBJ/P3qfRtcXv4dR2CWBkoBc0rCnJSdxTlcEovdpdH1ddTZiYxB4JBmPd1JSPlaWRsHjch+nkYzvS/7Hpd1H57/UsPt8Anit626Cw4jep9H1Be/eNvJ6fr801v4s9ZgD8mySUN/e5nOmglIqq3Yf8qy4yka3XS4uKrtZvfrofRpd3+od1onDnt8vUbR7yuL7rr9TPtX9f6rwbMKQRsi2k/oiEL1Po+vrq7dRG5nA1ySudA7ZP1f2r+TGyaPHm91/l7GN49gEovdpdH2xexd1XRIoCsgKuB/Vulvlr1Irte/pijfa3qYy8n4IRO/T6Pr66WmU9kTgdLPZPJdgr2lwrjfB3rngXQB+p/PSF3gvtO+1Lap/gSUfHi35H4BH6c6t036vZXuvPGyKrF3azLJJn5Z0WGR9wfuV90vJC2zhcw59Teg6L2PwQflpcUBeWDvVQQACEICACAwDctGUBdQgAAEIQKA+AQJyfcZ4gAAEIFBEgIBchImTIAABCNQnQECuzxgPEIAABIoIEJCLMHESBCAAgfoECMj1GeMBAhCAQBEBB2Tf0/uLtrkFXYoq4yQIQAACELg3AS8T4Bh84oB8ps3PxUvrUmiXBAEIQAACKxHw03Ecg7cBeSWfuIEABCAAgTkCzCHP0cEGAQhAYEUC2fWQh1r0E7+r3bHXszjXdqWyvetfDK9dYl++vO6Gf6d/of1vS9S5VB2Rtc21Ubqb9umcNts60MdrMteJE/Ze3y+pKbX0FwdkCfACRF5M6LNFKffCNF+UX2qrFpR3fj7Jl7909OJCnvMOkSJrKwEk/U36tESbz4mqL3K/o6301XX/89ZgWxSQJcSrET1Rvg3Gbor2/SQRH3sy+tJlNZL9qN7Xrlv7vyrziCREiqwtB0jam/VpTpvtkfVF7ne0lby6DjtnDbalc8gOiFPLbF6r/KWEerRM6otA9D6Nrq+v3kZtFwRKA7Jvy5i6TznN49pO6otA9D6Nrq+v3kZtFwSyAblw9Ms9zF109/+LjN6n0fV11NVI7YxANiCrPSnYptHwsIlp1MyUxZBK/P3ofRpdX/weRmGXBEoCcknDnpacxDldEYjep9H1ddXZiI1BwAH58U5KysfK0ih4XO7jNJJJz9mbOoeyeASi92l0ffF6FEU9E3iWxDsgf98dpDzZtrnm89JUxdS0RCqrdh/yHTEcLEIgep9G17dIJ1AJBP5N4GvaLZ2y8GpEZ+miQZ5GyLaT+iIQvU+j6+urt1HbBYHSgOyfK/tXcuN0oYKbwYhmbOc4LoHofRpdX9yeRVm3BIoCsgLuR7XwVvmr1FLte7rijba3qYy8HwLR+zS6vn56GqU9ETjdbDbPJdhrGpzrTbB3LngXgN/pvPQF3gvte22LG+VVk3x4tOR/AB6lO7dO+72W7b3yZimythwUaTfLJn2a02Z7ZH2R+x1tJa+uw86pwVZ1ehmDD8pPiwPyYfK5CgIQgAAE5ggMA3LRlMVcZdggAAEIQGAZAgTkZThSCwQgAIGjCRCQj0ZIBRCAAASWIUBAXoYjtUAAAhA4mgAB+WiEVAABCEBgGQIE5GU4UgsEIACBowk4IPue3l+0zS3ocrQjKoAABCAAgUkCXibAMfjEAflMm5+Ll9al0C4JAhCAAARWIuCn4zgGbwPySj5xAwEIQAACcwSYQ56jgw0CEIDAigQe3ceXfuJ3tTvf61mca7tS2d71L+5Td+7clr571pbTXmIXe6934vVELrT/reSaJc9p7f+YtrTU3tL3Mcx8bXTttfQVB2QJ8AJEXkzo8w6YF6b5ouNLbVWDckvfbutciqxtTnfOpna5fz9p85e9XtTJ3zWsllr7P6ahLbW39H0MM18bXfsa+ooCsoR4NaInyrfBeAfv2+7Yk9GXLquRWvrOtSeytpz2nF1t80j4tc/T/q/KPEpeLbX2f0xDW2pv6fsYZr42uvY19JXOIfuNObXM5rXKX0qoR1O1UkvfuTZF1pbTjh0CEAhGoDQg+7aMqfuUPYpysr1Wauk716bI2nLasUMAAsEIZANy4ei3yj3MLX3n+imytpx27BCAQEwC2YAs2SnYptHwsCVp1FxryqKl72E7p/Yja5vSSxkEIBCcQElALmnC05KTKp3T0neuSZG15bRjhwAEVibwSP4e73ymfCwhjYLH5T5Oo8T0nL2pc44pa+k7pzuytpPdlMpfasR9Pr281nU3uYb/0+2wO6yH4XYYN131LF3pgPx9d5DyZNvmguzb27w/9cZOZVXuQ27pe9v4mT+RtVm29Sm7mGkCpj0EYLcHTKYYbhlA+81fk6l0ysKrEZ2liwZ5GiHbXiu19J1rU2RtOe3YIQCBYARKA7J/Nutfa42TR2A3u/+MY9tSxy1959oQWVtOO3YIQCAYgaKArID7Ubpvlb9K+rXv6Yo32t6mshp5S9+59kTWltOOHQIQiEfgdLPZPJcsr1NxrgCzdy54F4Df6bz0Bd4L7Xtti+pfArX0rTbOpsjaZoUXGNU2fwLwP15/OnLu14f7+1q298qrptb+j2lcS+0tfR/DzNdG115Dn+r00hQflJ8WB+RjQXM9BCAAAQj8SGAYkIumLH6sghIIQAACEFiaAAF5aaLUBwEIQOBAAgTkA8FxGQQgAIGlCRCQlyZKfRCAAAQOJEBAPhAcl0EAAhBYmgABeWmi1AcBCEDgQAIOyL639Bdtc4vlHFg9l0EAAhCAQIaAl2BwDD5xQD7T5ufipXUptEuCAAQgAIGVCPjJQ47B24C8kk/cQAACEIDAHAHmkOfoYIMABCCwIgGvh1yc9BO/q93JXs/iXNuVyvauf1FcccGJLX3n5EXW1rv2ntmavfR7rRivCXKhfa9RvVpq6TvXyMjactpz9mPaVhyQ5cQLEHkxoc8WpNyLzXxRfqmtalBW/c18u61zKbK2Od22RdceXd8+vtLt98Ynbf6i3Asz+XuaVVJL37kGRtaW056zL9W2oikLOfNqRE+Ub4OxxWnf/+19vJ2MdlmN1NJ3rj2RtfWuvXO2fsqOH4flb87/yPXFknb5bOY7147I2nLac/al2lYUkCXmtbapZTavVf5SYjwiqJVa+s61KbK23rX3zDbHHjsEJgmUBmTfljF1n7JHyU6210otfefaFFlb79p7Zptjjx0CkwSyAblw9FvlHuaWvidpDQojaxvInNyNrj26vkmoFEJgAQLZgCwfKdim0fDQbRo115qyaOl72M6p/cjapvQOy6Jrj65vyJJ9CCxGoCQglzh7WnJSpXNa+s41KbK23rX3zDbHHvsDJfBI7X68a3vKxyjSKHhc7uM0kknP2Zs655iylr5zuiNr6117U7a7KZO/BPE+n/x8V8VNDjz2mAQa9/mzRMUB+fvuIOXJts0l1LfReH/qxZnKqtyH3NL3tvEzfyJrm5G9NUXX3lqf/QvURY4j9n8OgcZ9/jWRLJ2y8GpEZ+miQZ5GyLbXSi1959oUWVvv2ntmm2OPHQKTBEoDsn/66V8cjZNHETe7/y5j21LHLX3n2hBZW+/ae2abY48dApMEigKyAu5HXX2r/FWqRfuernij7W0qq5G39J1rT2RtvWvvmW2OPXYI7CPgOeTS5NHwO71R0tTFCx3/rOM1vsho6TvHJ7K23rV3y1bvC4/wPWhJnyy97ovfK9fK3+c65hh7S9853ZG15bTn7Eu07XSz2TyXIy/ec64Kq3w5l2sIdghAAAIPlYDirtcK+qD8tGjK4qGCot0QgAAE1iRAQF6TNr4gAAEIzBAgIM/AwQQBCEBgTQIE5DVp4wsCEIDADAEC8gwcTBCAAATWJEBAXpM2viAAAQjMEHBA9q1uftTM3IIuM1VgggAEIACBIwh4mQDH4BMH5DNtfi5eWpdCuyQIQAACEFiJwEv52T6blCmLlYjjBgIQgECOAAE5Rwg7BCAAgZUI3GctixP9tO9qp8sL0p9ru1LZKj+3buk71xeRtfWuvWe2Zi/9XprA61pcaN/rLK+WWvo+tpHRtdfSVxyQJcDrXfyu/LNhK/fCKV4w5VJb1aCs+pv5dlvnUmRtc7pti649ur59fKXb741P2vxFuRcX8vc0q6SWvo9tYHTta+grmrKQEC9+8UT5NhgbvPb9397H28lol9VILX3n2hNZW+/aO2frp+z4kU7+5vyPXF8saZfPZr6PbUd07WvoKwrIAv1a29Qym9cqfymhHhHUSi1959oUWVvv2ntmm2OPHQKTBEoDsm/LmLpP2aNkJ9trpZa+c22KrK137T2zzbHHDoFJAtmAXDj6rXIPc0vfk7QGhZG1DWRO7kbXHl3fJFQKIbAAgWxAlo8UbNNoeOg2jZprTVm09D1s59R+ZG1Teodl0bVH1zdkyT4EFiNQEpBLnD0tOanSOS1955oUWVvv2ntmm2OP/YESeKR2P961PeVjFGkUPC73cRrJ+L7kGqml71x7ImvrXXtTtrspk78E8T6f/HxXxU0O/D/ZDreDe/dZutIB+fvuIOXJts0F2bfReH/qxZnKqtyH3NL3tvEzfyJrm5G9NUXX3lqf/QvURY4j9rsE4HaXxz2OvqZzS6csvBrRWbpokKcRsu21UkvfuTZF1ta79p7Z5thjh8AkgdKA7J9+pseZDyvyKOJm959xWL7kfkvfuXZE1ta79p7Z5thjh8AkgaKArID7UVffKn+VatG+pyveaHubymrkLX3n2hNZW+/ae2abY48dAvsIeA65NHk0/E5vlDR18ULHP+t4jS8yWvrO8YmsrXft3bLV+8IjfA9a0idLr/vi98q18ve5jjnG3tL3Mbp9bXTttfWdbjab5+LgxXvO5azKl3MGTYIABCAAgR8JKO56raAPyk+Lpix+rIISCEAAAhBYmgABeWmi1AcBCEDgQAIE5APBcRkEIACBpQkQkJcmSn0QgAAEDiRAQD4QHJdBAAIQWJqAA/LjXaUpX9oH9UEAAhCAwH4C/1rLwgE5rWGR8v2XYYEABCAAgaUJ3Hsti6UFUB8EIAABCIwIMIc8AsIhBCAAgVYECMityOMXAhCAwIjAfdayONFP+65213tB+nNtVypb5efWLX2PmP1wGFnbD2InCqTfP5/32gsX2vdawKFSZH1oO+ylEplbSYtq6S8OyBLg9S5+V/7ZgpV74RQvmHKprWpQVv3NfLutcymytoxu998nbbfavADOmbYwSVzD6kPbYS+TyNxKWrSG/qKALCFe/OKJ8m0wtnjt+0kiPv6g7dJlNZJ8NPOda09kbQXaPRJ+7fPUjl+VeZQcJklTWH1oO+xlEplbSYvW0F86h+w3rpcOHKdrFbyUUI9maqWWvnNtiqwtpx07BCAQjEBpQH4p3f5oO04exTjZXiu19J1rU2RtOe3YIQCBYASyAblw9PtTjXa19J1rT2RtOe3YIQCBmASyAVmyU7BNo+FhS9KoudaURUvfw3ZO7UfWNqWXMghAIDiBkoBc0oSnJSdVOqel71yTImvLaccOAQisTMAB2bes/aItjXbHEvaV+7w0SvR9yTVSS9+59kTWltOOHQIQiEPgT0lxDD5xQD7T5lvXUnDV7r+T5krTVMXUtEQqq3Ifckvf/yYwvRdZ27RiSiEAgaAEfHOAY/A2IJdodAR34B6nFMRtr5Va+s61KbK2nHbsEIBAMAIeIZck/6w2Pc58eL4f034zGC0ObUvtt/Sda0NkbTnt2CEAgWAEigKyAu5H6b5V/irp176nK95oe5vKauQtfefaE1lbTjt2CEAgHoHTzWbzXLK8VsS5AszeueBdAH6n89IXeC+077UtbpRXTS195xoWWVuBdo/w/Y/Vn36cu//dn9dq13vlTZM0hNWHtsNeGpG5lbSohn7V6eUhPig/LQ7IJWI5BwIQgAAE7kdgGJCLpizuVz1nQwACEIDAIQQIyIdQ4xoIQAACFQgQkCtApUoIQAAChxAgIB9CjWsgAAEIVCBAQK4AlSohAAEIHELAAfnx7sKUH1IP10AAAhCAwGEEnqXLHJC/7w5SnmzkEIAABCBQn8DX5IIpi0SCHAIQgEBjAgTkxh2AewhAAAKJAAE5kSCHAAQg0JjAo/v410/8rnbnez2Lc21XKtu7/sV96s6d29J3z9p61x65381W+rwWjNfcuND+N5dFSZG15Rg9VO3FAVmAvACRFxP6bJjKvRjNF+WX2qoGZdXfzLfbOpcia5vTbVt07VH1SZdf+5+03Wrzwkxn2kKkyNpygNBeuEC9QHk1oifKt8HYYLXv0YCPtyvdu6xGauk7157I2nrXHpmtX/vaXmvzY3f+yLFe0x5ZW44D2gsDskC+1uZlGcfpWgUvBdIjhlqppe9cmyJr6117z2xz7LFDYJJA6Zd6fuaTP56Nk0fJTrbXSi1959oUWVvv2ntmm2OPHQKTBLIBuXD0+9Nk7UcWtvSdkx5ZW+/ae2abY48dAnMEsgFZF6dgm0bDw/rSqLnWlEVL38N2Tu1H1jald1gWXXt0fUOW7ENgMQIlAbnE2dOSkyqd09J3rkmRtfWuvWe2OfbYHyiBR/p46C/rTmfan0bBU6ekkUx6zt7UOceUtfSd0x1ZW+/ae2abY/+Pte+mmv5SA+/zidl3q9y0htJSu3z7IdLeThyQnyv3fb6TDzmV3bf4yDwJOYGvch9yS99u8FyKrG1Ot23RtUfXl+P7UO3uN7X9osf2t9Qu3/96yGnplMWfgnw2ATqNkG2vlVr6zrUpsrbetffMNsceOwQmCZQGZP801L9IGif/N7zZ/XcZ25Y6buk714bI2nrX3jPbHHvsEJgkUBSQFXA9v3Gr/FWqRfuernij7W0qq5G39J1rT2RtvWvvmW2OPXYI7CPwaJ9hotyj4Xd6o6Spixc6/lnHa0zIt/Q9geJOUWRtd4ROHETXHlafXvcewXtQkj45el0Xvxeulb+fYL1aUWRtOQgPXfvpZrOZ/VIvBxA7BCAAAQgcTkD/hO79pd7h3rgSAhCAAASKCBTNIRfVxEkQgAAEIHAUAQLyUfi4GAIQgMByBAjIy7GkJghAAAJHESAgH4WPiyEAAQgsR4CAvBxLaoIABCBwFAEC8lH4uBgCEIDAcgQIyMuxpCYIQAACRxEgIB+Fj4shAAEILEeAgLwcS2qCAAQgcBSB+6xlcaKf+F3tvHlB+nNtVyqrshbyuFUtfY+1jI8jaxtrHR9H196BPi894HUtLqTV6wGHSdKDtkq9UYttcUCWAC9i/7vyz26jci+s4gVVLrVVDcqqv5lvt3UuRdY2p9u26Nqj6pMuv/Y/afOTTby40Jm2EAlt9bphDbZFUxYS4sUvnijfBmM3WfseDfj4g49rpZa+c22KrK137ZHZ+rWvzY8e+kWc/8ixXtOOtnq012BbFJDVxNfappbZvFb5Swn1iKFWauk716bI2nrX3jPbHHvsEJgkUBqQX+pqfzwbJ4+SnWyvlVr6zrUpsrbetffMNsceOwQmCWQDcuHoNz1bb9LJoYUtfec0R9bWu/ae2ebYY4fAHIFsQNbFKdim0fCwvjRqrjVl0dL3sJ1T+5G1TekdlkXXHl3fkCX7EFiMQElALnH2tOSkSue09J1rUmRtvWvvmW2OPfYHSqAkIKdR8BSiNJLxfck1UkvfufZE1ta79p7Z5thjh8BeAtmArPm8NFUxNS2Ryqrch9zS915iO0Nkbb1r75ltjj12CMwRyAbk3cV/Kj+bqCiNkG2vlVr6zrUpsrbetffMNsceOwQmCZQGZP80ND3ufFjRhQ5uBiOaoW2p/Za+c22IrK137T2zzbHHDoFJAkUBWQH3o66+Vf4q1aJ9T1e80fY2ldXIW/rOtSeytt6198w2xx47BPYReLTPMFHu0fA7vVHS1MULHf+s45uJc5cuauk715bI2nrXHpatXvcewXtQkj45el0Xvxeulb/Pga9pR1s9urXZnm42m+eS78V7zuWsypdz9fBQMwQgAIG+CSjueq2gD8pPi6Ys+m4u6iEAAQj0QYCA3Ec/oRICEHgABAjID6CTaSIEINAHAQJyH/2ESghA4AEQICA/gE6miRCAQB8Ehre9/a1v+caqP6vMC4WTIAABCEDgSAKKp3+rirN91Tgg+1Y3P4pmKnEb3BQVyiAAAQgcRiA9KHry6v8DrMXoISD5ct0AAAAASUVORK5CYII=\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": {
Martin Bauer's avatar
Martin Bauer committed
150
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABbsAAAA0CAYAAACjIeUlAAAACXBIWXMAAA7EAAAOxAGVKw4bAAAgAElEQVR4Ae2d6/XctNbG/2SlgAAVvNBBCBUkdJBABcnpAFY+5XxjQQeQCrh0AFTApQM4FQDpgPf5eaTB4/HYkqybZ7bW8sgXeWvvZz/WlmVb89Y///xzZ8kQMAQMAUPAEDAEDAFDwBAwBAwBQ8AQMAQMAUPAEDAEDIEYBP773/8+UPmX7pz3XP5c+9/EyNlj2Vu2vWd/3XeOOdPxFkh5ZrTtMAQMAUPAEDAEDAFDwBAwBAwBQ8AQMAQMAUPAEDAEDIFQBL7QGOJ/fGGtf6X1X7W87/ddcX7Ltjd366Ux7XvSDAL+PVm+aK6xKWAIGAKGgCFgCBgChoAhYAgYAoaAIWAIGAKGgCFgCGRAQANj32l5kkGUiThF4MUEV8YU39O+h6fFrnLrlm3P4lDx5IGWH8gTBH6nc6Zj2t/dd4I+ktAfE4TaKYaAIWAIGAKGgCFgCBgChoAhYAgYAoaAIWAIGAKGQLcIaMyLQTEG1GzsK7+XeKv7l/xidyHxlm3P4iBdk2+0fCZhPyl/zHaoYJX9aFxW2y+0/cwPdo+P2bohYAgYAoaAIWAIGAKGgCFgCBgChoAhYAgYAoaAIbB7BDQAxkD3H8q/3r0xHRowgysDwOD9W4fqZlXplm3PCSRc0fK5ZGaZ/oZpTCwZAoaAIWAIGAKGgCFgCBgChoAhYAgYAoaAIWAIGAJXhYAG0HjT86Fy3hy1VBgB4czUJU+1fFC4qu7E37LtOZwh/L6XnB+V83BqU7LB7k3w2cmGgCFgCBgChoAhYAgYAoaAIWAIGAKGgCFgCBgCvSGgQbP3pBN/lvisN92uUR+HN/N1f6D14KkorgGLW7Y9p/+EI18F8HCKh1TJyaYxSYbOTjQEDAFDwBAwBAwBQ8AQMAQMAUPAEDAEDAFDwBDoFIEfpNeXGji7+uk0WuPvB3uVD3MoK+cNb+Zi/qO1bqXrv2XbC2HLVxivheu3WpIemthgdyHPmFhDoDYCagR4Yv2rcpuHrDb4O65PfOEfj3/SEvVHEDs22VQ3BAwBQ+AiAhZLL0Kz2wN7iXPGvd1SrLnixvHmLrhZBXrnnvTjzdB3lNv0JYVZKoz9G/SfaZ1BbhJv6F499rds++DlAj/C9HstLyWarwTgUXR669WrV7/rrP9IkP0jbTR8doIh0AcCun4Z6CaQ2+dZfbjkTAv55ol24h/a3A+10BHo4im39KBD8lqLDXgLBEuGgCFwmwioLbyZWLoUk3SMh6AMELyv5XdtfzlmhLaZR/G58qQ3bcayaq1L167jnPS7Ge7V8nmJeuQn68slAntLHF/iiY5Z+5rIoZTThPffOu9r5Vc/4JqIz8U2LZarDmv4fZK0/62THZ1sSK+btH3J7k5cc1TD6cqXGe9rPXjcRGXpwz67f5RkK4aAIbBLBNzF/Ej5zf0BxF4cJt/wRPKB8uGppHKCK/8y/HYPNkgf/vmYG20GMIbPznrQy3QwBAwBQ6AWAmoD6RjfRCyVrWsx6QuV4UUY/lyKB6HHwW6H01Plu3q4Ln27jXMO05vgXq3ruUQ98tPadVOi2mCZxvFgqIoWDOCJta9FPfCvcPniU20x+Pr5v3ttzSOQm6uS18V9rbdvKb9V2wPsXoKt+jHpyx9VMv0QD6ui3+62P6is7jKr0BDIh4Auft5UovP9OJ9Uk5QTAfmIwQIGBsYNNE8mGfxm0LuLJF2G6W+U0zG0ZAgYAobAzSCgdu9mYqlsXYxJOs5nyDyMJfHwc/omDfvO5j3VeQzgsHxKzsm9JenVXZyTTjfDvd74EKOP/LR43cTIKlnWOF4S3XXZazzR8dT2lXuGF1p827yuTOUS0q279lUQcO/FVAi7+QqpltuEyWKbpuNJXK2l/5Z6btX2Nbu3YFr4XB5W0f6dfTWwVq8Ndq8hZMcNgb4R4I0rPs2yIN6vn/DR9MafDgQputE+nFbsl04hgxVev2IVmWBDwBAwBDpC4JZi6VpM+ksxwP/3x8fy0TcTP/GQ9mSfyhPj/lTO9Fy8Bf6D8h8m5/Wy2VucuyXu9cKBFD3WrpsUmaXOMY6XQnZd7hpPUtpXHojx5RH3DL3dN0wR6YZ7ikHgxv3M51MlbXtAIDtXd4Trrdq+ZnevLvTTbdMORiWbxiQKLitsCPSDgII4T2QJ5Jvf6pYsOk8vnXV+oHM383H2qr/0YlAAbL912PoMv5HO3o477G7zK33/0MIgB1Oa2HQmbdxgtRoChkBFBNTmZYulMWqr3upxV3WuxiSVeYMdo7J+4Jt99A/Q2994UJTEDcixL6JyfHbKgPd7WqZvhg8ntPpBHy1dxDnpkZV7kledUzn92Kv+0mv1usmJw1ZZ0vcqOd4rP7y/QniiMtHtq87hXoFpmGgvuk7SsRvuCSjua984/LLhJnnZ29kSMpcMVn2rbZrKRHOVOnWef8HrT22+q+2u5kqXPiVthxu8JMA0cF1NLxtit/TuMkl3rmP6nTxMO06rF6Ls/ZBCVsYQMAS6RIBg8iUNQAbthvnjvBzJZLCTT+X4c6o9pF71Z05TOqhTH32i/XQIuxoEcI6GV/whWXeDFE4/ywwBQ8AQyIlAzlgao1eLuBUTkyhLnBrHL24S77Tv+KCWWKFd3OCNy1GMbR7sWpwDjfmUm3stODVvWdreXvWPuW7SLM9/Vi99uZwc75Uf3nsxPAlqX73gneW9cI94NX0wmwPKEjwsIXPJ1iJcVX8A3/OV1zAgqfyJFh589/QCVSnb6e8MfSTl9Il6SzF296Y7+vC/Yl+JS1HjE/d6tMR0MgQMgWUEdKHToHKDmevTrBeS6RtoKidY0Zj4N5DZ13PqVX8w/WUMnMMUXAk63SXpx8CEf3ranX6mkCFgCBgCuRBw7XHOWBqjWou4FROTwOU4qO0MI24NgwfCDv25ofM3dX+5Mj5j+x2/0VPeQ5wrxL0WnMrp2l71j7lucuKRLOtKOd4rP7yfYngS2r562bvJO+IesemHAsCV4GEJmUuml+IqX3kdHzCIC6wz4A3fe0lFbJeNvNzGIH+PD/jBPsbuXnw11sN/JY8dwelecElXsDOyxqpfrHwPuPSgQzGAOxTcGG8+zeIz4TeZoOGzkJNB2Uxya4npTn/5hk4Wwf0Y4N0+nkzyedN0EKEWViH1oGP0vFghgkPLCJ8jbqHnWLltCLTGvHX929Czs1MQ6MDnuWNpDAxV45awjo1JJzdsOp8bjEda/P73tW+pD8JAN3X2mlrHuRLcq8qpAo7tTv+E66YALMkir43j3fHDeyaBJ74dHUQktK++6l7z1tzjC1rSceD1sJnltwQPS8icNbYUVyWX+zZi/rRfwHYXL8+Vsn0W6I52JtjdkfYHVWQDPGKJ+krgfowlquRTladxPmmgY2RccVnegn2qJWoemVx4mG9yIRklp6XPmbeNwJgliT/HOTmdQGTz+XLPA7JH2zvV3z95ZGCbtpPEtDDPdoArnUP++Z027XsUr5kcXhZraoJ+qKtZm2Y+r+/sTmpsxjlnf9ZYGoOpOF877sbGJObZfC09h8+StU6bzByUfG5NTPtGC8nf2DK47dfZ/0BLz/cLTeOcsMnOvQacws/ZUqf6x1432fDIIOiqON4pP7ybYnkS2r56+XvLW3OPwVXm+c0eg0rwsITMBcKU4ioxn9TzV16lbD9Y3u9vrN29WsKLmVEPTu6HWqKLkE4Zk8zPDuZqPxXzFO8DrY87u6FV7LqcbOYt24daqg8OUafAu+ibXQPbsfKtfA7PHCx0JLInJx9OdfXHCqGGdqQ/Tx75pImO1my7GWpT7XLo7PTGhqqD3arX2rPaDnf1Cfsmccx83sjhHVTbinOYrrqLxtIYeJ0upeNuVEySTvTl56bbOtmncn5eb3+jOza92wfmTm/ic4s4V5x7lTg19nXW9Y70j7pusoKwUdg1c7wjfngvRfFE+ge1r1743vKW3HNYPVJe/IvlEjwsIXPCn9pcfUf1z/UPJmpV2axtexWjAiqJsjtAXqsi9Cn9wH2QDvdCSumig6AvlfMU8pjYr+U7LfyZ3WstN/3ZuXBgQAucql3Qrq4z3xydZCtFEWjhcxlEAL9T3SWeVnMN8xbXLh9aCZOe9KcxLvJAQnJrJHQfuFajMuqQ/2ZjTa36rZ7BB1XjmPncWCcOVOXcCPFisXRUx+qq7K8Vt0rGJN5SP96AyCbW/cPeVQwaFqge55ytRblXkVNFXNeZ/iWvmyL4TYReHcc744eHe+888XbkzJtwT/zgXoIl+33yGJwSPCwhc6yzWy/F1TdOPoPb41TcF+PKVtZL2b5SbfPD12L37yCp6wR7gtL9oFKHwS8GtE+SKoLUwxseWueTRv+mwkm5G9sAJwYLs00xsYIfdZ35ZuUcO5wXgdo+H57O5TVhaDiGG25dy8i/U871XOQTsNy6O32z6S/baUT5vH76mXmQ6jqPwI4+Jf4YJUiHDIUIKLXn7bb2LIPjMoio2aaZzzM47ApE1OSch6tILPXCQ3LFimxxa6m+0jFJ8j/TwvQmXM9/amHKrsdLOnVyrEWcw/Ri3JMPqnCqlP9y6i9Z1pe7u7sqjufkRy4OS6eifX5nM+MKns/Ey9+1v/evRltxjzaQRP1FUgkelpA5NV51FOOqZBf7ykuyPfeT7svBoaTtU5xzbm+1fa92X8DQf63BGBUP01bT/dUShwIfC6hag7eBKvVZTDh9reVvLXT8eRhQOplvSiO8Ir+Bz2nwgy7wFdWPh2UDHQM6T/DWP7Timj/5muN4QmcrBfSnM8CSmvxbW1n9lKpM4nnDGxHwQUutT9GtPUt0Vs7T5O+accx8ntN5O5VVmXMepeyx1AsOyWVzzbhbPCbJnl30Fya+aRHnUKEI9ypzagLl9s0C+ltfzr3dKmxr9uWKcLwAP7aT9iChaPsqu2mnrH0N95Yf7C7yZncJHpaQeQGuolxVnf4rr+G+UXYR63J85bW1LQeO0rZTR4m01fa92j2Hpb+meaEiKN1fK+VI6gWvFbfjBwTA62MtyU+fQoA034SgVK1MFZ87a2j0cl+Tv0omcsmPSRzby0Ou3vTnyePwhvwRzP2teI7RaSw+2G3tWXcEKd6mmc+783lrhYpzbmJgiVg6qWJxs2bcuoaYtAhm4sGqcW6kYynu1eTUyJxsq73pfw3XzTVxvDd+eOJfA0+8LTnzVtx7xxnh689pE7JK8LCEzDm7i3JV/fqev/IqbTv3yz1+gVHU7jmSldonfjHjAOL9Nb5a1epgtyQwTUl3bye6m2R04xOVD7VwcWVp1CSHeciH6Vkk9yzpGGRmeoJnWp8bBAIvBrqKDnZTv5ZZ30gvOtJMQcCTj5NPndyx18ov2qhzkpPk8hQP2We+WdKLCnWcPzl9rvwN27US9Wq5iIeO9eJzD8nqp1nSOZgDKvu2F1wrV53ZeNJC/yWcpA/8nb02l87LcUx1b+WyV8O3p8EBxZ+YmF9sz5Anuzyf2fxQy3MtXJefaCH9rDLfH1bz/kpuNq7m1eyyNOm8lQc14lg2n4/4cfUxb8lWGIHvlVWPo75u1Z8aS2twDjXHKSaWct5i27PkG3fs2PfSdrW4q7qaxaQx2DnXZdPWNg51ase5MQSr3KOw7AyOfSpbjVPeENWZLT620N/bMZdLn2bXjerOwW/M6prjsjMrv0fyqvUFWvJkjrdb910B9+AUies3KI14Q/m1OB/Uzo5krnJRZYNkBhmzUEj1FG/TVEeXXyGUtl3yaWu7s7203Qt0u1PdueLYtBp/jU/3n23fO9tzvoNX34M6ZOenltkj4JgTkIHm/2hhviqmX+CJ2OYkeU8l5JsAQQyy0MGbS+DF8dJpyTf8aSXYMCj/cqIIb51jZ/akOtd8w5yOs3ppP4PzT5UHB6ccBqi+3fhcunpe/RVgexMOBOh1Jzt2x5MQu1qXycTlwQzJ8tdhcEDZaP9Se4booe2QXrQfP2vhT5GfaJuOBTrCqexJ8nfHVem8lzYtp8+btHeN+NFdHOXCy8C7Wn0ndI2JpTFtTxMeZm/4OheYgWuDhZJTO87Fcg89Y/hX1XPCb3fxsSpAiZXl4jfVS1bvHM/Nb2uDE3nn+JKj/zho0IJ7zvR3XR5yr+yKFmlnjYseXctvDoGccWwCHjEt+EW8+5OT5zYZRIhpLOZkZNvngGNAlKdkPvEk5YH2MfCx9W3KjyRjceoGHWcCfgZc/NNyr4fPwQvcSqdZ30g3buIYDCLxhvlUT/ZtxQnZJ0n1EiAv+kbH0MM/lLik128nQrUhudjJAD0PNz6YHs+wvTefr5osnJpwYFUxFZBuRXgSUvcNlMnB5SlMvtM43Z97e7Y9oxJxhgdh08FseOSfoBP0pse1a1sqxVXJ9br+KQ3f1ba3Y5vC/56dgwc14lgWnwu/Ju1dC344W3uMo7BvK+9qcM5fJXBvNQnv4LanFQ9XjbjOAlu5NkWlVpyj3iDuUTCGf5SvmaSb9eXKAZ6b32jaHcdz81vymvQFytGgieS9cw/QgttYCufmoZNpXAQIS7eMQIm2BDyj7hVCBrsZRHjTkad4m286MECDQopq3A6n/Purxo43tfn89yRp/w9aGJwdJ95AOhuYdQUY1PU6jc/JvX7JN39JX/85P4PEn08qxs7pvkmRpM0136CXn9olSC+VfyhN0Je0yb8HEae/kr9Hn2PE2jXZigOnAM9vZefJfDXre+V/vgrx/BqfwLV1p+NzD774o42Ln+m78/4ZC1tal6y3lo6HHpOcXFyeVpn9uptW4LYvtWcc/kX2jR/a0R7jh2Gf8kV/OPkpWXauSlcGuv9UzgPTO+VPtMzFmBR9B3k6cS9xLJfPW7V3LfiBrV3FUc/jDLyr1XdC5aGNV74WS2PanlY8xJ7FJM4Ex6RFQYkHVX+WOEf1klUi1tWKc5gQyj3KxvCP8jVT9vYvVXlxwvpy6+D1yPHc/G7SBot/V9G+FmpbYWZN7lHf0MbKnrX4TllSbh4iswkXqXgp9c5V6VekLQeTW7U9xm6VzdJXk5wS/bQxtX0/arxvdv3+7N5OdzrgaDC/naj40G1fGnyeFL+4yT9VDwMQvoS2kT0eZPGH3tcxf8Pp9/kcB4Q2sP6cbLn0GupWDtHA66ins4d9Wd/sHtV10TcrevFw4EwvnYNPf1POGyQl0t58ztMsElhdTCtYw+kzrC8Ky3hAenlOZuUJKko2NkV9AaBz5gazkQXf3lN+0h5QT0jSeVmCRUhdozK5uDwS2ceq8Jy27fBo+tAzq7KqsxRXeVP0sVdW9fyohcFu+DYXa3zR0DwXD1rHsWCfC7c3gKPc+2y3MU9mLPJjxdZWcRT4c/CuJudCY2nXPAT4kCTetIhJIaqllMnBtZR6c50TxD0qk9+C+ZdLuRA50su3tSX6ctFfPkkf68uFOK5emSCO5+a35FXvCwCp6r2W9nXvbWsSw3Pz0HGiCRfXAOidq9KvSFvufNL1dVrK9kY+76Ytubd2Ueg4N98MIvWQeHuPgc+hARkp9InWmVpk60DBeGoUL54O3Q9+g1z1cEO5lMDLB/qlcluPrflmDi/s4Z9Mpx3orbrM1YXMOd9QFn+N/YheYJtbL8Qupb35fMmWuWNzfinFgbn6p/vm9KHMJp6INwzgM0jEtcdyiykXlwfshKnHsdZ/Nqy1Z16voa3Qhv96hXaDaay8vrl8n52r0tEPRr6ZKMk2HM6RcvEAPEvHsdw+n/NZqfZuri78l9yWRfKD+nuJo9idg3c1OIeuSUn+CW175rhRiodJtuz8pBxcGyAYxY1acS4Z+gj+JdcRceIcxzk9uf3jZNnov3z6TOu8bMCD4JN7MMpdecrGb3ASfrSrpK45Lj1D29eDNcu/c/y0NngZM45eC/eGvuuI++uWuxKZeYhU42Iw+lbwihDI2pbM4BJ8fxoy2I2wtcHdGR1Od6U0OKcShi0C1S/j/ZLLAAELjckxJdbHm3XHAQe3TsfrOIWJk8sn4p8fKztfeUe7VgfeE3Uc17bmG/w21QNbzt7qzqBLsG9UP3r9pmWc8N+gl3R5EatPbPlRxdV8vkHHkbp3DIqR4FhICuZAiLAMNhThifTiIRg3RlO+h5jVVZkNGOfissfDc8xzzu+fzTfo7eXNtmfI1cINL9whDU/9tT32NX8Cc6Kntv0N3uGs+N8SXPU6Yes4se3xHvZv0D8XD1bj2AYdve1ZfS6hwe1dBt2b8sPZ+psH0uWt4ijV5+DdKueoKIPvEOPbi5PrjgM+UY+W6LZH5wfz0Ne1lGeyd6mK6sc22pSDa95m73/PB79/Nt+ot5fp6/J1+/0nOXVpSeHfiZy5jQx2lGj/UJWXFo73KNKT9SfKuaZ2lTZgnJPfYOZ55nm3iOMGvcdyfV2+7vGxYZ16tBThtyoIboMz2XtmX8sdG2y6Bu4B/Sr/vH8K85BqgrnodVrKN/h2SawdMwRmEdjAt9xtyVg/4oq/xsf7Z9fvze493cnN1Ienu+K2HFB/K/d/qhQnQKWdDBoMliG5fQw888eFx5s+tz+lPmT8pPP5Z2jk/qSFm0cGX7/TwjxC/9PyldaXQOZPFBffRND5D1QmRUeddkxrvhkPBt2pTjqnLCe6bdXFnR/kG6f5nF6PdMzvZ4qYJXydmEO2Uf8qPt+o49jev9wG/AlJHtOhrPSY5UCIoK02uPOb8STExtZlNmKci8seBs8xzzm//yzfqLeXh/5zscZzlvnv0Al9ju2D9nH8Zy3HtFUfd35NrhK4Pd53G/XPxYPFOLZRR++rbD53AoPau626u/Ob8WPB1lZxFJVy8G6Rc1Sy1XfIcIl2hHS87g6bJ7/RbY87O4iHJzVd2Mho74Ua6u/OYFMOrnnDvf89H/z+szyD3l6mr8vX7fdP81T+TeWcbG+1w52fvf2TXGSCyZsThQ/bDyf7ut7ciHFOfoOT55nn3UXsNuo9luvr8nWPj/n1Ivx2woPa4Iz2epua5xttugbu4QP+/J20xL9DicO4CFwMvsfwJwbmQVwMkbXRtyFVWBlD4IjARr7lbkuOermVk+tqenC8fW+8cWGdgdGLnQwBwSAwZV668391+z718rRNxwWleIob0vD4U8c5DRGJge1PWbTOW9fPtH6cn5MC2k6tj7e1eXMc2XS6PpAsPpX/TAv1Mx/w59P6tG+aKHv8xH56kO0NOo7FLfpGBdGbuhi8xyYG7knHtybY0LFUvDidFOybQ/EzveAEN7nvOD2/ceWCso36V/H5Rh2PODg5+Gvu85BjudFKEAdG5S+uZrChKU8uGtbRgY0Y5+KyR4QBWBJ8W0wb9fayL7VntFe08fDnherize7nynnoOExdo/ykvdV2r22ax9JjK1OG9EC/x8C9Uf9cPFiMYxt1dGYPD17n+hfRPncCg9q7DLqXasuC+HHBVjjUKo6iUg7eLXKOSjL4DjFeDngvxdKiPBwUWfnJZe9KNVUPZ7ApB9e8zb4t9tee33+WZ9B7kOnkrHGPsqn8O9N9vCODHaXaP9owkh8oPWwdtr2f/L6u840Y5+Q3OHnseuN4EX47YtTqC3THw1vnnnPIKtdHjivJQ6oJ4uJIn4urG317Ua4dMATmENjIt9xxbKwifYVpP2F8/GT9rVevXjGHFwPIXOyzSccow6Ayo/TJSefzp2/8GVdMIzTUp3N4q/qRcm7mgtKW+oIqmCmkOhkk57OspRuo45lbddT5wb5RWR4OPFE+i2GqLjov2jdHAAJXnG4M2l/ENVX/QBUuFlO9wT7PoaNk8IUET6CP0+tcVG5yQOcscmBSfHYz1Qad1wVPZo2a7HQ28glO0h9UenE6H278x23TOLMNjy+2t5TV8eS2kvNzJOnAQDI+e1vrQW32Vr11fnB7FmJjqj46rxhXJftv6f5Y+TGeav0f7eOLluOAN/ZpuwkPVG+1Nk11ZfU5uPkk2YvtXSq+Oq8Lfng7Q3Nn7+7jKPam+m6MlWQkx9KxnLV11bPIw7XzOZ7D3pB6xmVUZ1L8GstYWm9h01Qf6VA9zqGD6q3Cvam94+1U/HVekfZPcnnwCS4nfQ7tJ2byolF0f0znEEOtL1e5LyfM74R9c46jB0m6LLbBOl69r6U6rX09uOfkN4cvnAy+0P9I64v3WyeVV9iQPotcDFEhB0Yh9fgyqq8oV309S7mzeXNbvlTH9FgPdqNTTtsli/jtEw9CeXls8R7f1Z80fusrypVLF8ZS6BN8qfXhQdIl2TpO/+7ZvUsFJvu5MP2AzeRQ1OaHqngR0AVpvEkQ22BtqW9BlcVDAA9eoWmrjrO+Ec689UhHY0hahxw4nSctl1KqLim+uaTDlv2p+m+pk3NjfJ5DR74+eLSmdCIH1sRyPNWGXngSYiPtVGpbNch31xx/suQX2lA6XzwMo2O9lFIxXpIZe4yHYtM/wFuTsVXv2fZsrdKF46n6lOSqf0t9UFtcoC7mnD8Z6HY2peq/AEnQoZptWhafC79riXkx/AhyZmShPXAOk3LoGRRLY/BL5GFIFTnsDalnKCM76DP62EUeE79C66lq0wWlWsQ5VMnOvQv2Le1Oxb9UfPR9Lv8mstcdLs7FR398KUeml7tU7uKxjddCKsYX9Uk4cFMcT2yDq/ppI6dCKVDVpgtKteKeby8YpG2WErkYom8131biaojNm9vykEp8mY7sRqXNtmOPFl4u+l05LznTp2NM8H9aX7tOqvFN+qwlr+ufawX98aDBboHAzRdPU3wF/vzgXOc+UOFgxcaC3bnUzWfuQWlLfUEVzBRy+IATeK2mHDq6uuZ8Q2f0m5ESr7X+tcqffO7vj6fq4s6L8o2vMySXfGxjMOSlFtYZ0GBKlpOkfcn8OhEUuaF6sT3I5xl15CFGyJRAURwIMT3VBndeMZ6E6B5TRvryBDPoOl6Qy8Mlpt0YD2x/68rD59nksEpqK2cFpu8cBmFDT8+ht8N8rj0LVeNYLlUfd14xrko+A8nvKh9PL/X4qHJrGRUAAApUSURBVLhbSdV/Kid2W/VWbdNU3+b+hbMxqr1Lxded15wfMX7Fp1quIo5it/NBjjYyNJbGwB3FwxDBGe0Nqc6XSYpf/uS1vJFNc2pVj3NOiRLcm7Nvdl8q/u68Iu2fZDNAxU09/flpOn4JNT2wtC2Z1pc7TAEXjF8qN2b80IrjUW1wRntnILi4y9rXGWhy+UJyPN8ZbG+ZorgYomgujELqcmWKcjVUD9mdoy0PrY5yXdiNIplsZxyQ6YKPX0hpneuEB+/cl84mlWkyvjarzGEn/Q+Sf6B12Fr4vb9wbHqI+Z6HT0KmBwK3Px4DHHiOLza8xarzY97s3lKfrzc259MAPy92yLm5dJzzDcTlBpeBYab+YJB4Cb9UXVJ8E4LNUEY6Q+aLF+FIUKr+IxFJqzE+z6UjA6bUSxCdfXjhLInlgDttMUu1oShPFjVud5Agwo0by5DE5zda3NbFLBXjiwJjD0hHghsBhQGy0JRL77n2LFSHcblUfYpzVfhamzb21CFubulfIC22vds7P04RXNgS364pjmJpqu+mKIXG0ul5S9uxPFyS5Y/lstfLC8lT41eIbMq0sOlEt8ZxrgT3Tuxb2UjFv3R85OEnfdthsEo+GtZdG7ZiUrHDqddCKsbZDLlRjse2wS38lMqpUG60sOlEt8bcQxf6PUN7daJY3Y1YLoZoV9u3pbkaYnOLMtdm91OBiE3TxD7GCnnTey7V5tucDuN9frB7zpZxueN60JzdvrQaLip4onzrG49eZFDuGkzm614arA2SVaqQdOMJEE+dgp805NRF9ZpvcgIaIKulz1U3XzkwxcSlxinAgnpFpCeDp8WuYcmH/2DBjRHzPtJG8anO8Qmmtpsn6UOwYVCPT4iqtqMxxo/0PJtHOkZOalnV36Q9Q1/VXZSrqZjUOE+2N4tjLX0eg+0t8yMGp9CyLTmHjqp/V7E0FNcS5YTVLuJXqO0je1rFud1xT5gVj4+qwz9k5+sNXtZhGp03oX6tUW7EHevLLQAunHbH8QVzih7aC6dCQRjZ06p95V6LMau3Q3W2cmEIjHzbdfsXZk14qb3aLb2J28xzzTjlyf+9aZuBbmLuB1oPHkBW+SZJOnJdP1X+1poCKsM97bP7awXHx3USA7nVB2hUL52cbge6wUg6VseFen0y33gk6uWNfc7FTuO0l8Huotew4z9P0LtN0pGB+MFnjbkTghHB8NI80iHnbyrj/NmkTVXdRbm6CZjCJ7fkZUufx8B6y/yIwSm0bEvOOR13FUtDcc1dTn7aU/wKNb9pnJOSu+NejfZPdVhfLpTB6+WM4+sYNS9h7WsRF/CghUExpv2kX28pAwJXytVVZPZsN/zXgo3T/8Ng37v8KM0dOxzp65e+aNSY8P2+9DdtDAFDIBCB4RNYGl8tURd9oHwrlgkB+Yc3zWmcP9TCU1Pmx+o9fSwFn/eupOlnCBgChsBGBCyWLgC40/i1YNHJodZxzrh34o6+N3Z6LRjHO6bVTjkVimhr7vl7Le6/lqb8DLXnpstdOVcv+vaK7OYa4FqYJsYoSHxR3XWSL9DxgRYeZAWne8ElraAhYAh0g4AueJ5SM0WH/9yzG91MkVME5CvekP5SC3NRf6PlV63zOXiXSbp9imLKrXPYpYdMKUPAEMiFgNo5i6ULYAqfXcWvBVNODsmu5nHOuHfiku439nYtGMe7p9Td3jgVimgn3OPlIuL7J6F6W7nLCFwrVy9bfDhyRXYPL7DJnuOAt9YZ6OYaIf11yLr+9bpHjU/YYHfXPjXlDIHLCKiR8n984Z/KXS5sR7pAQD6jgSawfKf1B10oda7ES+2yt7rPcbE9hoAhcIUIWCwNc+pO4leYMXd3XcQ5416ou/oqt5NrwTjeF20WtdkJpxZtGB3sgnvSh69n/ADZSD1b3YLAlXE1GIo92y3dGXv4Py3PtP6pFl6W5E3pn7WQmvzn4KHq4F+m5eL/6qJ0tcHuYHytoCHQJQIMeL/uUrMbV0qN8UOWGRh+cfu664BJX952I5BEPTWdsdF2GQKGgCGwJwQslo68tcf4NVJ/cbXDOGfcW/RY24N7vBaM4205s1b7Hjm1ZpM/3hn3+F8E5uzu7n7L49V7fs1cXcL+Gu2WTczdzZ+K8rU5f/zMvT5/Ak2KGkA+nFL9l6/iv4qt1Qa7YxGz8oZARwiooRr+xE/58EluR6rdtCryB29t/8ri1rvHQ3ryhJe3MZhuxZIhYAgYAjeDgNo/i6XO23uMX6FE7THOGfdCvVe/3B6vBeN4fZ7E1LhHToXa1xv3pA//acUgnt3XhDpxVO6auToy82z1xux+JAB+lM28+d1tkn4vnHJDXz1GURvsjkHLyhoCfSLwWGrxpI7BSksdIOCCBoHj+5kAQmAhdfPHotKRwXn+8OG51vfwdBf8LBkChoAhkBMBi6VC08Ws3cSvUAJ0HueMe6GOrFhub9eCcbwiORKr2hunQs3smHtM1/DC6RdqjpUTAtfK1TXnXqPdsumplr+1cL8/JK0zbsQX6Hxd1ntCx7kxlVW9bbB7FSIrYAj0jYAaK25Kmceo53mg+waxjHY0zCf/GCxf8QkOgYaHE/itl8Snfl9IJ5u+pBePmB6GgCFQFQGLpSdw7yl+nSi+sNFtnDPuLXit/aE9XQvG8fZ8CdFgT5wKsYcyXXJPbStvgnK/xZerluIRuEauhqBwbXYzsD39E0qmBGE8gj9z7TZJP6YhQn98Ep3eevXq1e86C0O7ecsw2go7wRAwBO50DdMQMDdZ143WLblKvqCBHn8+h48YVO6mvZUuDL4/6kmnW+KI2WoIGAJ9IaC20GKpXCIcuo9focyRLbuIc8a9UI/WLbeHa8E4XpcTW2vbA6dCbeyde9KPKRAY2Htb6z29aBQKcdNy18TVGCCvzW7Zw1cOPtEn4iXJbsYjvGLTXDry4iDzjY/HU6bFzrZVnuv+mQ12n0FjOwwBQ8AQMAQMAUPAEDAEDAFDwBAwBAwBQ8AQMAT2jIAGvvgfJeYmTno7dM+2m+6GwF4R0PXKSxcMdkc/qPKD3ff2arzpbQgYAoaAIWAIGAKGgCFgCBgChoAhYAgYAoaAIWAIXEDgufa/0AAYX25ZMgQMgX0gwPRIz3TdJn+RYYPd+3C0aWkIGAKGgCFgCBgChoAhYAgYAoaAIWAIGAKGgCEQiIAGy5ji83MtDJ5ZMgQMgc4R0DXLtcrXGJv+T8wGuzt3tKlnCBgChoAhYAgYAoaAIWAIGAKGgCFgCBgChoAhEI+ABs2+1Fl/KB/PXRwvyM4wBAyBogjoGn2qCt5THjVP95xSfs5uLvpvxwUkPPl18bEcWzcEDAFDwBAwBAwBQ8AQMAQMAUPAEDAEDAFDwBAwBFohoDEu5gDmz/m+bqWD1WsIGALzCOi6fKgjr7U8jhmPVln+dHOc+IPKj/yb3fxD7d+jxZ54jaGydUPAEDAEDAFDwBAwBAwBQ8AQMAQMAUPAEDAEDIFdIqBBsY+kOPMA8+d3lgwBQ6ATBNyAdfRAt1OfaU/OxrP/H14qQ++rikB2AAAAAElFTkSuQmCC\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": [
Martin Bauer's avatar
Martin Bauer committed
251
       "<lbmpy.methods.momentbased.MomentBasedLbMethod at 0x7fd8406872b0>"
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
      ],
      "text/plain": [
Martin Bauer's avatar
Martin Bauer committed
285
       "AssignmentCollection: d_8, d_5, d_7, d_4, d_0, d_2, d_3, d_1, d_6 <- f(f_6, F_1, f_2, f_5, f_7, f_0, f_3, F_0, omega, f_1, f_8, f_4)"
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
      ]
     },
     "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": [
Martin Bauer's avatar
Martin Bauer committed
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>38.68 ms</td> <td>114</td> <td>67</td> <td>0</td>  <td>181</td> </tr></table>"
314
315
      ],
      "text/plain": [
Martin Bauer's avatar
Martin Bauer committed
316
       "<pystencils.simp.simplificationstrategy.SimplificationStrategy.create_simplification_report.<locals>.Report at 0x7fd840687b38>"
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
      ]
     },
     "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
341
   "outputs": [],
342
   "source": [
Martin Bauer's avatar
Martin Bauer committed
343
344
345
    "simplification_strategy = create_simplification_strategy(method)\n",
    "simplification_strategy.create_simplification_report(collision_rule)\n",
    "simplification_strategy.add(ps.simp.sympy_cse)"
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
   ]
  },
  {
   "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
363
       "<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>"
364
365
      ],
      "text/plain": [
Martin Bauer's avatar
Martin Bauer committed
366
       "<pystencils.simp.simplificationstrategy.SimplificationStrategy.show_intermediate_results.<locals>.IntermediateResults at 0x7fd84078d828>"
367
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
      ]
     },
     "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
396
   "version": "3.8.2"
397
398
399
400
401
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}