From b0285f2b570c722fef9e181f6ecc8085fd9bf31c Mon Sep 17 00:00:00 2001 From: Adrian Kummerlaender Date: Wed, 10 Jul 2019 00:16:10 +0200 Subject: Add basic talk slides --- notebook/seminar_talk.ipynb | 1360 +++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 1360 insertions(+) create mode 100644 notebook/seminar_talk.ipynb diff --git a/notebook/seminar_talk.ipynb b/notebook/seminar_talk.ipynb new file mode 100644 index 0000000..3f9276d --- /dev/null +++ b/notebook/seminar_talk.ipynb @@ -0,0 +1,1360 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 1, + "metadata": { + "slideshow": { + "slide_type": "skip" + } + }, + "outputs": [], + "source": [ + "from sympy import *\n", + "init_printing(wrap_line=True)\n", + "from itertools import product\n", + "from IPython.display import display, Math" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": { + "slideshow": { + "slide_type": "skip" + } + }, + "outputs": [], + "source": [ + "from sympy.printing.latex import LatexPrinter\n", + "\n", + "class Printer(LatexPrinter):\n", + " def _print_list(self, l):\n", + " items = []\n", + "\n", + " for expr in l:\n", + " if type(expr) is tuple:\n", + " items.append(\"%s &= %s\" % (self._print(expr[0]), self._print(expr[1])))\n", + " else:\n", + " items.append(\"%s &= %s\" % (self._print(expr.lhs), self._print(expr.rhs)))\n", + "\n", + " return r\"\\begin{align} %s \\end{align}\" % r\", \\\\ \".join(items)\n", + " \n", + "def aligned(expr):\n", + " return Math(Printer().doprint(expr))" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": { + "slideshow": { + "slide_type": "skip" + } + }, + "outputs": [], + "source": [ + "# copy of `sympy.integrals.quadrature.gauss_hermite` sans evaluation\n", + "def gauss_hermite(n):\n", + " x = Dummy(\"x\")\n", + " p = hermite_poly(n, x, polys=True)\n", + " p1 = hermite_poly(n-1, x, polys=True)\n", + " xi = []\n", + " w = []\n", + " for r in p.real_roots():\n", + " xi.append(r)\n", + " w.append(((2**(n-1) * factorial(n) * sqrt(pi))/(n**2 * p1.subs(x, r)**2)))\n", + " return xi, w" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": { + "slideshow": { + "slide_type": "skip" + } + }, + "outputs": [], + "source": [ + "def eq2assgn(m):\n", + " return list(map(lambda expr: (expr.lhs, expr.rhs), m))" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "slide" + } + }, + "source": [ + "# Symbolic generation of GPU-based Lattice Boltzmann implementations\n", + "\n", + "#### Adrian Kummerländer" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "slide" + } + }, + "source": [ + "## Lattice Boltzmann Method\n", + " \n", + " \n", + "$$\\begin{align*}\\Omega(f) &= \\left( \\partial_t + \\xi \\, \\partial_x + \\frac{F}{\\rho} \\, \\partial_\\xi \\right) f && \\text{Boltzmann equation} \\\\ \\Omega(f) :&= -\\frac{f-f^\\text{eq}}{\\tau} \\Delta t &&\\text{BGK collision operator}\\\\ (\\partial_t + \\xi_i \\cdot \\nabla_x) f_i(x,t) &= -\\frac{1}{\\tau} (f_i(x,t) - f_i^\\text{eq}(x,t)) &&\\text{BGK approximation} \\\\ f_i^\\text{eq} :&= \\omega_i \\rho \\left( 1 + \\frac{u \\cdot \\xi_i}{c_s^2} + \\frac{(u \\cdot \\xi_i)^2}{2c_s^4} - \\frac{u \\cdot u}{2c_s^2} \\right) &&\\text{Discrete equilibrium} \\\\ \\overline{f_i}(x+\\xi_i,t+1) &= \\overline{f_i}(x,t) - \\frac{1}{\\overline\\tau} (\\overline{f_i}(x,t) - f_i^\\text{eq}(x,t)) &&\\text{Discrete LBM BGK equation}\n", + "\\end{align*}$$\n", + "\n", + "…nothing special." + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "slide" + } + }, + "source": [ + "## Experimental symbolic differentiation\n", + "\n", + "```\n", + "% derivative of x is 1\n", + "diff(X, X, 1) :- !. % cut to prevent further rule applications\n", + "\n", + "% derivative of a constant is 0\n", + "diff(C, X, 0) :- number(C).\n", + "\n", + "% sum rule\n", + "diff(U+V, X, D) :- diff(U, X, DU), diff(V, X, DV), D = DU + DV.\n", + "diff(U-V, X, D) :- diff(U, X, DU), diff(V, X, DV), D = DU - DV.\n", + "\n", + "% product rule\n", + "diff(U*V, X, U*DV+V*DU) :- diff(U, X, DU), diff(V, X, DV).\n", + "\n", + "% quotient rule\n", + "diff(U/V, X, (DU*V-DV*U)/(V*V)) :- diff(U, X, DU), diff(V, X, DV).\n", + "\n", + "% some special derivatives\n", + "diff(sin(U), X, DU*cos(U)) :- diff(U, X, DU).\n", + "diff(cos(U), X, -(DU*sin(U))) :- diff(U, X, DU).\n", + "diff(log(U), X, DU/U) :- diff(U, X, DU).\n", + "```" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "subslide" + } + }, + "source": [ + "```\n", + "?- diff(x*x,x,D).\n", + "D = x*1+x*1.\n", + "\n", + "?- diff(3*sin(x*x),x,D).\n", + "D = 3*((x*1+x*1)*cos(x*x))+sin(x*x)*0.\n", + "\n", + "?- diff(x*x*x+2*x+log(x),x,D).\n", + "D = x*x*1+x*(x*1+x*1)+(2*1+x*0)+1/x.\n", + "```" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": { + "slideshow": { + "slide_type": "fragment" + } + }, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAIMAAAAYCAYAAADZPE7mAAAABHNCSVQICAgIfAhkiAAAA91JREFUaIHt2VuIVlUUwPGfF0LC6GJlQT0UChKWEfQQFI1aFEEXy5d8qBOWBEVEFFEUDUFplFRKFj3EFPVgBgUWVtCFTBFG0hho6MpQlDM42ZWKbvaw9jBfZ777d+b75uH84TDf2WfttdbZa529195DSUmPuAeD+BkHsR1Le+pRSc94CzeIBDgTr2IUx/XSqZKZwXz8g8t77UgJs3P3j+AdfIPfcQj78AAWTIP9o5IPh6ZBdz0W4EYxM30h3vUnfIi1po5Lt1iNzdgpltLDeLHJvgvFh7WpKGf+xB48hw3JscHk1Lc4tShDia0i2eYUrLcRN4t3+g4vYb145x9T+yuY1WWfYH+y/wuGtZYM65L88qKcmVej/aFkaEtRhvAoDmBRm/0z4VNfG31XiKUpPwOchK+T3mt64NdyLBaJ2Ke1ZNiBcR18WPnB+KOG3Mvp7+Jc+9vC4atz7bMwkJ5tqKJvI67DSjFNd5t3xU7m31z7KJ5Jv/u66VDiPXwuxq0VjhYJvl0sFRO0G5+63Jc6bsy1L0vGP/H/jNyY5J+tomsTxnBGq07kyLT/BdbjrqT38Tb7Z4rxq0/zM8OaJJsvxNuJzxTuRL8YkJ2p48c4oYrsQHqepft70/1WU2eeLaIwWiGm5IlrfjNO5cgUnwxzMZT0XtKmjkz3k2EbflV9mR/QfHyqMpo6TFw7RLVajVNENT6CW5P8mziiiuzhGld/M07lyBSfDI8lnW90oCPT3WSYJwrObTWetxKfuizEKnwqKu9zasitNxnYXTiyVUMNGFE7kapdA23YuC31Hdb8Idh0+tWnuWS4IsmtqSPTVHzmNjA0JvbiH+EzvKD68fHBit9r8VsDva3yBI7JtZ2NK/G8CEol+1vUfwueFGvrSs2fe0y3X82wShwJ1JvNCo/PPpFZx+farxVV+YH0/OlODTVJppjp+PakZwgndqiL7i4Tc8R2ckcdmWmJz1hSdmxF22UiK4dEcTmMv7CkCIMNyHQ+6HcnHftMTfJ2yXQvGZYnmXU1nrcdnyWiss8z2+Sh066K9vPFdPMVTk5tq5Pca42MFUCms0G/P/Xfq9h/lGW6lwybxdaxWnHfcnwqa4ZLxangB/gS3ycjF+J0scO4Kckuw+viPP9iMQURx7h7xZp5gdiWzkSux4NiIHeK4jHPiPYK0U64Kl1MfpjnVfgxLrb9lfK7xaxdScfxWYqnRJEzjr+TskGx9Zv4ehaJxPgBZ1XRc5HIvj21DBVEpv0vsF/jyv/9GejXSIXsuantjpyOmRKfki7ysAjqab12pKT3DJuerWpJSUlJSUlJSV3+A/N7UzIz/6I6AAAAAElFTkSuQmCC\n", + "text/latex": [ + "$$3 x^{2} + 2 + \\frac{1}{x}$$" + ], + "text/plain": [ + " 2 1\n", + "3⋅x + 2 + ─\n", + " x" + ] + }, + "execution_count": 5, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "x = symbols('x')\n", + "diff(x**3+2*x+log(x),x)" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "slide" + } + }, + "source": [ + "## Characteristic constants" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": { + "slideshow": { + "slide_type": "fragment" + } + }, + "outputs": [], + "source": [ + "d = 2\n", + "q = 9" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "metadata": { + "slideshow": { + "slide_type": "fragment" + } + }, + "outputs": [], + "source": [ + "c = [ Matrix(x) for x in product([-1,0,1], repeat=d) ]" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "metadata": { + "slideshow": { + "slide_type": "fragment" + } + }, + "outputs": [ + { + "data": { + "text/latex": [ + "$$\\left [ \\left[\\begin{matrix}-1\\\\-1\\end{matrix}\\right], \\quad \\left[\\begin{matrix}-1\\\\0\\end{matrix}\\right], \\quad \\left[\\begin{matrix}-1\\\\1\\end{matrix}\\right], \\quad \\left[\\begin{matrix}0\\\\-1\\end{matrix}\\right], \\quad \\left[\\begin{matrix}0\\\\0\\end{matrix}\\right], \\quad \\left[\\begin{matrix}0\\\\1\\end{matrix}\\right], \\quad \\left[\\begin{matrix}1\\\\-1\\end{matrix}\\right], \\quad \\left[\\begin{matrix}1\\\\0\\end{matrix}\\right], \\quad \\left[\\begin{matrix}1\\\\1\\end{matrix}\\right]\\right ]$$" + ], + "text/plain": [ + "⎡⎡-1⎤ ⎡-1⎤ ⎡-1⎤ ⎡0 ⎤ ⎡0⎤ ⎡0⎤ ⎡1 ⎤ ⎡1⎤ ⎡1⎤⎤\n", + "⎢⎢ ⎥, ⎢ ⎥, ⎢ ⎥, ⎢ ⎥, ⎢ ⎥, ⎢ ⎥, ⎢ ⎥, ⎢ ⎥, ⎢ ⎥⎥\n", + "⎣⎣-1⎦ ⎣0 ⎦ ⎣1 ⎦ ⎣-1⎦ ⎣0⎦ ⎣1⎦ ⎣-1⎦ ⎣0⎦ ⎣1⎦⎦" + ] + }, + "execution_count": 8, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "c" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "subslide" + } + }, + "source": [ + "## Derivation of weights" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "fragment" + } + }, + "source": [ + "Weights of third order Hermite polynomial:\n", + "\n", + "$$\\eta_{-1} = \\frac{\\sqrt{\\pi}}{6},\\ \\eta_0 = \\frac{2 \\sqrt{\\pi}}{3},\\ \\eta_1 = \\frac{\\sqrt{\\pi}}{6}$$\n", + "\n", + "Weight formulation:\n", + "\n", + "$$\\omega_i = \\frac{1}{\\sqrt{\\pi^{d}}} \\prod_{j=0}^{d-1} \\eta_{(\\xi_i)_j}$$" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "metadata": { + "slideshow": { + "slide_type": "subslide" + } + }, + "outputs": [], + "source": [ + "def weights(d, c):\n", + " _, omegas = gauss_hermite(3)\n", + " return list(map(\n", + " lambda c_i: Mul(*[ omegas[1+c_i[iDim]] \n", + " for iDim in range(0,d) ]) / pi**(d/2),\n", + " c))" + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "metadata": { + "slideshow": { + "slide_type": "fragment" + } + }, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlMAAAAVCAYAAABi32uwAAAABHNCSVQICAgIfAhkiAAABzlJREFUeJztnWuoFVUUx3++wqKHWZEEkhZC4Y1I0jJMt0VShqFSXwIhpPoSFIhF9KBjlFlRaFpEVJr5IeghfQgrE063J0pZUWkvtRdpXtOyNM1HH9Ye7pxhzjkzc/acM/vM+sFh7t2zZ85ea/9n1r77dUFRFEVRFEVxggGOhj6bO1oaRVEURVGU4nAqte2ko8GJgTGZ3wUWAMsi6dcCS4H3gL/sTVYlLMDpwGHgiVDaw8A64GdgP/AHsBG4Dzilwb0uBV4FfgMO2OPbwPSEZWkHLn01AJgLfAzsBfYhfroVGOSuyE5waXdWfUDxNeL6WfJJI2Hm0P9CurFJXp/tLmt9a8xIjsaMdHa71Ack18g+pH20APix3s2MNaRS5/xn9vxeYBPpKvtmm39qKO0gUtnPA4sQh26w+X4FRsbc5x57fiewHFgIPGOveyRhWdqBS1+ttL/vAJ4FlgBf2bRXkAenKLi0O4s+wA+NuH6WfNJIwEhgD+KDJI0pn+0ua31rzEiOxox0drvSB2TXSJVQz1QYQ+PG1FRgDFIRQd6klb0G6KO2VTy0Tt4H7b2fiqRfZ9PXAifEXDckYVnagStfzbTXbkG6FwOGAKvtuRtcFNgRLjWSVh/gj0Zc+sk3jYDY/Q7wA/AoyRpTPttd1vrWmJEcjRnp7HahD2hNI1UyNqbi8iYx+iSk22x5grwA59NvXMBARBz/AKclvE9RMGT3VfAXxi0xeXvsuU9aL2IuGPLRSJw+wF+NGFrzk48auQ04AkxG3jfNGlPdYjeUs75BY0YaDBozGuFCH9C6RqqEGlODM9wgLVcDxwCvJcw/wx6/CKVdAoxGuih323v2AP8C64GPnJS080R9NcIet8TkDdLGAcOQIRNfSaOROH1AOTQS5yffNHIu0gW/BOgFLktwTTfYnQW1OxkaM8oVM1zoAxxrpB2NqVlIyy/aKgyYDxyPtDYvBCYhRi8K5RlvjzuAT4HzIvfoRSay7XRT5I4R9VWfPY6OyXtW6OdzkLFiX2mkkST6gHJoJM5PPmlkMPAi8BNwV4rrfLc7K2p3PBoz+iljzHChD8hRIwb3w3xDkYllLzfIs53aZYZrkJn6YR6y5w4B3wGXI84aC7xpz1UTlLsTGLL76np77ffA8FD6YGTlQeCzqxyV1SUGNxpJog/wVyOG1vzkk0buR1bgTAylVWg8zNcNdocxlKe+wxg0ZiTFoDGjHq70Aa1rpEpomC9uawSXTEMKt7pBnhHI5LMRwGyk9bwR6YoMCCaZDUBaiuuAv5EVCrOAX4Ap1L6kfSPOVy8hQjgb+BpZYbAYWQExHREASIDylWYaSaIP6H6N1POTLxqZgPRGPUa67nPf7c6K2l0fjRlCGWOGK32AY43k3ZiahSxXfCNB3h2Ig6Yhe0KsDJ3bbY9bgM8j1+0H3rI/T8hc0s4T56sjwDVIt+V2ZF+euUglTwJ22Xy/t6+YzkmqkUb6gO7XSD0/+aCRYHjvW+DelNf6bHcrqN3N0ZhRvpjhSh+Qo0YMbof5BiHjt2vSFgRpRR6lf2nnbPv7hjr5g+XVd2b4rrwx5OOrY5EK30exlvgGGPLTSFQf4K9GDPn5qSgaGUZtt3ujz+LQdb7bHYeh++s7DoPGjKQYNGbE4VIf0LpGqrRpNd9kpDXYqDuuHmfYY9AV2YuMa45BZvEfjOTvscdtGb6rCGTx1Rxk/PgF4L88CtUGsmokqg/obo1k9VNRNHIAeK7OuXHABcD7wDfUDgH6bndW1O70aMxojuqjTTHD4LZnailS8LiJX+fQv4QzzED6N9j6IHJulU1/IJJ+BdK1uQf5CzjMCjq/SZmhNV+dGJM2Htkqfy+1KzQCVuC33Vn0AX5qxNCaPsBfjUDjCejdaLeh/fW9gs7XtUFjRlIMGjPicK0PyKaRgCoZe6Zm2g/0F3oiUgkg3W/zI/k/RMYto1yJdKH1Irsg70IcNAWp6O3ATZFr5gEXAXcjLdT1wJnIGOphmz+6b0YwJ+xQc/Oc4tJXa5Gu2S+RB2EsMpHwANJNGbefiO92Z9EH+KMRl/oAvzSShm6xu9P17fv7ADRmaMxwrw/IppGmGBr3TFVoPNdhWyjveJs2r869eoAnkRUGfUgF/omMXVaoXdIZZjjwOLAV6ZLbBbwOXFwn/0bknyeeXOd8XlRw56vbkR1r9yAPw1bgaWBUg+/33e6s+gA/NFLBnT7AL41EqRDfM9VNdlfobH37/j4AjRnbQnnLGDPy0gek10hAlVDPVBhD8mG+Ziy094rbOKxdDENalkX6Z5ZxuPZVWe3Ogg++ysNPandx0fdB5yirr3ywuwj6iFKlSWMq+Gxu4Us2IS3ETjID2RY+bhy1SLj2VVntzoIPvsrDT2p3cdH3Qecoq698sLsI+gBZDRjtPQNks6qAUdROPusDluVfNkVRFEVRlMJzHHBHJK3SgXIoiqIoiqIoiqIoiqIoiuV/+NxET1RKZ7IAAAAASUVORK5CYII=\n", + "text/latex": [ + "$$\\left [ \\frac{1}{36}, \\quad \\frac{1}{9}, \\quad \\frac{1}{36}, \\quad \\frac{1}{9}, \\quad \\frac{4}{9}, \\quad \\frac{1}{9}, \\quad \\frac{1}{36}, \\quad \\frac{1}{9}, \\quad \\frac{1}{36}\\right ]$$" + ], + "text/plain": [ + "[1/36, 1/9, 1/36, 1/9, 4/9, 1/9, 1/36, 1/9, 1/36]" + ] + }, + "execution_count": 10, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "w = weights(d, c)\n", + "w" + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "metadata": { + "slideshow": { + "slide_type": "fragment" + } + }, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAA0AAAASCAYAAACAa1QyAAAABHNCSVQICAgIfAhkiAAAAHZJREFUKJFjYKACCGFgYJjMwMBwmIGB4RMDA8N/BgaGJYQ0XYAq/MzAwHCdWE2ODAwMqgwMDIwMDAwOuDSxoPH3EzKVgYGBgYkYRaOaBlwTeuQGQDEDAwODBJS2ZGBgWABlv2FgYChBN6SBAZJ0cOEH5LiMzgAA6XoX52TB9a4AAAAASUVORK5CYII=\n", + "text/latex": [ + "$$1$$" + ], + "text/plain": [ + "1" + ] + }, + "execution_count": 11, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "sum(w)" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "subslide" + } + }, + "source": [ + "## Recovering the speed of sound" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "fragment" + } + }, + "source": [ + "$$\\sum_{i=1}^{q-1} \\omega_i (\\xi_i)_a (\\xi_i)_b = c_s^2 \\delta_{a,b}$$\n" + ] + }, + { + "cell_type": "code", + "execution_count": 12, + "metadata": { + "slideshow": { + "slide_type": "fragment" + } + }, + "outputs": [], + "source": [ + "def c_s(d, c, w):\n", + " speeds = set([ sqrt(sum([ w[i] * c_i[j]**2\n", + " for i, c_i in enumerate(c) ])) \n", + " for j in range(0,d) ])\n", + " assert len(speeds) == 1 # verify isotropy\n", + " return speeds.pop()" + ] + }, + { + "cell_type": "code", + "execution_count": 13, + "metadata": { + "slideshow": { + "slide_type": "fragment" + } + }, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAADEAAAAdCAYAAAAD8oRRAAAABHNCSVQICAgIfAhkiAAAAi9JREFUWIXt1t2LTVEYx/GPGTTFNJkhk5SoKeRCopAkaoSkhD/BtXArLrwkF8pbjQvuKcqUuZIkbpTxFlFKTSF5yUthzMTF2jNzzmnvffY5Z585Ls63dmu113p+63n2evazFk3+D6ZUOL8d3xq4fi7swuZGLJwn5zC90U6U0lLB3NboGa6TL1VTSRBrcb9ejkwWJ9DVaCfiqGQnuvCpXo7UQtYgevAyYewkbmEIP/EZgzgs+87NxSjOVKObtU4fQD9exYwN4yGe4wNmYDVW4m3UHyqjvxd92IjbOeoW0Zcy1pbw/hj+4kIG/QF8FKpfxbqF6dSCo9hSYtQp/V/4lfD+StT2pNhCh7AD/UJKVaxbGEQv1mB3idFW3CzjSBzbo/ZJmXnbhAP0Wl66y4X8KwzuvGwF4CCO4DTuClv+GHPK2F3FD8npU5XuG6yL+tMUV4w03kcLjD0DQtVJow3fhUCq1o37wjewI+pvwJ0yjozRLVS7buzEIqEkrkix6cVMXM9Z1yYTpfRUtEg1LMBvPEuZczma05Gzrqn4giWylcc0BoUUmB0z1iqU1YFadePSaUSoRofwtIoFCpkXtaMxY+uFkzctlarRHWePEOn8MmKLhVwtpcXEoXQvwfZs5ETcz1+L7jjtsl279+GPcMe5KNx0L+F1tNA7LE2wHRJKZt66RczKMGeZcI48EvJ7BF/xQKjtnQl2qyJn9uesO6kcF4JY2GhHauGF8JWbNGlSB/4BwBaSc+fGZlwAAAAASUVORK5CYII=\n", + "text/latex": [ + "$$\\frac{\\sqrt{3}}{3}$$" + ], + "text/plain": [ + "√3\n", + "──\n", + "3 " + ] + }, + "execution_count": 13, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "cs = c_s(d, c, w)\n", + "cs" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "slide" + } + }, + "source": [ + "## Moments" + ] + }, + { + "cell_type": "code", + "execution_count": 14, + "metadata": { + "slideshow": { + "slide_type": "fragment" + } + }, + "outputs": [], + "source": [ + "rho = symbols('rho')\n", + "u = Matrix(symarray('u', d))" + ] + }, + { + "cell_type": "code", + "execution_count": 15, + "metadata": { + "slideshow": { + "slide_type": "fragment" + } + }, + "outputs": [ + { + "data": { + "text/plain": [ + "array([f_next_0, f_next_1, f_next_2, f_next_3, f_next_4, f_next_5,\n", + " f_next_6, f_next_7, f_next_8], dtype=object)" + ] + }, + "execution_count": 15, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "f_next = symarray('f_next', q)\n", + "f_next" + ] + }, + { + "cell_type": "code", + "execution_count": 16, + "metadata": { + "slideshow": { + "slide_type": "fragment" + } + }, + "outputs": [ + { + "data": { + "text/plain": [ + "array([f_curr_0, f_curr_1, f_curr_2, f_curr_3, f_curr_4, f_curr_5,\n", + " f_curr_6, f_curr_7, f_curr_8], dtype=object)" + ] + }, + "execution_count": 16, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "f_curr = symarray('f_curr', q)\n", + "f_curr" + ] + }, + { + "cell_type": "code", + "execution_count": 17, + "metadata": { + "slideshow": { + "slide_type": "subslide" + } + }, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAosAAAAXCAYAAABj0QUGAAAABHNCSVQICAgIfAhkiAAACedJREFUeJztnXusHFUdxz+XtlCkhRrrA2xtRcWKF3u1qYDYuqlPrMaiVRtfbKJEbH2gYLTVaBUViBR5aGMqaovYkCgJsY01SMj1gcVQoDyspdVya1SeQqSg9kpb//iesXPnzMzOzs7rbn6fZLJ3z8ycPedzZ+ecOa8FwzAMwzAMw8hACzgU2nbWmpruORa4EhgBRlEePldnghqAOfExJz7mxMec+JgTH3PiY07iGQ9epjO2Hngo2DEx5uBfAcPAo1WkrEB+DLwN+DlwLfA08LNaU1Q/5sTHnPiYEx9z4mNOfMyJjzmJZzx4+RfwFfd3G5gVd1AL1SJXV5GigpmD0v6LuhPSIMyJjznxMSc+5sTHnPiYEx9zEs949DJMqGXxiPrSUSiL3Ov1taaiONron9TqIQ5z4mNOfMyJT785gd69mBMfc+JjTuLpKy8tsrcsvhbVOp9ATal7gS9SfeXzXfj968E2p+K0FEmb/BenOfExJz7mxKdfnUB+L+bEx5z4mJN4xrOXYTqMWezE14BVwN3Ad4GpwDLgQmDAvVbFA6h/fTnwTODrLvwQsLvCdDQJc+JjTnzMiY858TEnPubEx5zE05deWnRuWfy0O+YSVDEMGHLhj5Peunieiz/rtiRDuiegQZl3Zzh2vNCmt2Zvc+JjTnzMiU8/OoHevJgTH3PiY07iGa9ehsnZsngCajXcCqwMRwJsR0vtzAFmom7pOM4jYXZNAhuAGzocczJwNHBHF/H2O+bEx5z4mBMfc+JjTnzMiY85iafvvLRIb1k83+1fmrB/m9s/o+iEdeBs97mfrPhz41gO3A/8B7gdWJDhnBGSxzTEbeszxNkUJwuBTcDfSb92oozQv05WAreh8b6PID+DGc4boX+drEBP3U+4bSuwOMN5I/SvkzCrUJq+nfH4EYr10hQnq/HT/mDGc0dizu0HJwDHo4aVR4B/A/cAr8tw3gj962SE+Dx8J+d5/XBPmYAa/II6yv1oWGFSo+EwOVsW3wEcRGsExfF84ClUMaiSV7nXOyv4rEnAfxPC3wlcgSqMvwU+BmxBTxV/SYnzcmBaJGwI+d6ALt4w2zOksylOjgHuAn5Id7PA+tlJC1iLKowDwFeBm9B18lhKnOPdCcR7mQT8Ffg8GsMzgG6uNwDzSO+66WcnQdhpwDl014VVtJcmOAnC72Nsd+CBjHH2q5NnA7egMmcxqjCeCDycIc5+dTIJmI8qRwGDwC+Bn3SIc7zfU9KcfBY9mJ+NHihegSq7++lyrkmL5JbFiagm+lDCufPduZs6fEYZYxZ/gyqxU2P2zUAyHkTpvxd4o9tGkcCAmS4PLwqdewh4H/BrJPRDKeG/B74X+fzdwEUZ8hClTW9jJJriJEw3LYtxtKneCXT20osTgCmowHt7jjy1aaaT4Py4/GdxAqo4f7T7LPWVk+OAP6NlN4bJ3rIYR5v8XupwEhe+x8VfFG36w8ktOdKfRJtynEB9ZQ+oEvgnxs61yEqb/iiPN+O3hG5w4XEMk6Nl8eXAUcCRaALLwcj+893rug7xFD1mcQCYiypl+yL7ZqLurDtQJeVhtOTPPtQ9vIOxNfAht29P6D3ABagbaDeawPOamPB9wA+ASyNpuDF0fFU0xcnjRWWoAPI6AeU5zUvQXZrXyVT0nUprVSyDMp0E72Fs/s+ICYs6mQC8B1Wif5crZ/lpmpN1wE+Bm4Ev9ZSz/NThJOn78yngI8DfUEF6KxrWMdJTDrunSU42IQ8bgTegnr2rUXfr/wv6CkhzAvWWPUcCHwAuo1on0KzyeAZwLppbshP1Zi0iY4NW1sriPPc6DXgLY7uiVwDvRSuTd2pZnJ3x87JyEips47rG16GJNks4XLnd5V4/gd98/ErU1RNcTHPR+I+lHP7nJIWfgAq5aMvrQ+gLXCVNcdIk8joBfSHTvPTq5HIX/60Z8lEkZTqB+PwviwkLOAXdOCejG+IS1FVSJU1ycg7wYuCDOfJRJHU4SQrfgn6OdhfwXFRR3IoaM6p82GqSkxNRGfwt4M0u/qvINj6vSNKcQL1lzxJUd1nfORuF06Ty+GKXlh2oN2siWspnbZaMdFtZ3ISedDeisREL0JPxNtTkWTXBWIDoLKNZqFJ7Kn4rKEh6tMt4CI2vC7/fjH8BJoWD/9QyEBNWNk1z0gTyOoHOXnpx8k00CWgB2cdeFUWZToL30fynObnP7Z+GFrK9BnX7FNnt2ImmOHkp8A10XYxmTHtZ1OEkKXxL6O97UEVxD+omvCwpAyXQJCdHoMmUK937O1EFZQXVVhaTnED9Zc+H0bVT9XwKaFZ5vBR4P6qr/cEddwWa6PL9DvnI/Isr89BNa5mL/EzgM8B09MstC6in2zHpHzGECt9tMeccjb5M0Vr7/EjYXNRnHyUu/FH3ec+LhD+H5HGeZdEUJ00ijxPI5iWvkzVoHMnr0ViaqinTCcTnP83JKPKwDRV829GwlSppipPT0b31XvQLWU+j2a3L3d9HJWehcOpwkhYe5inUSvKSDscVTZOcPIAK/jA76G64VxGkVRbrLHtmod69q1OOKZMmlcdr3HYdetj6EXrIWpWU+CRaxE9wCRaUvL3bCGvkrai2fmzMvkGUz+NDYQtd2Kvd+2PQP/L0yLlJ4aAJLtExm7vIN8GlDOpwEqbXCS5lkOYEOnvJ6+RK9BBxcq5Ul0uvTiA+/1mvk4CbgWszHls2VTuZ5uIMb7ehHp1B8g3SL5qynKSFR5mMKkt1jeeMUoeTjWgCRZgLUYWxKdRZ9qxG10ieX6srkzqc/AP4eCRsJZpEF8cwCT2jLcauH7TThZ/i3kebRJvMs9AYlo1oPMscNAZoELX+HQTOcseeivJ6ANXoQaIPIPFhksJB4zZH0QDsl6GxaE9S/RNeEnU4mYKeoIbQNXSB+/sFRWSoANKcQGcveZysRWsJLnLxB9uUIjJUAL06gfj8pzm5GPVOzEb3m4vcZ5xZQH6KoA4nUYbpbTZ00ZTlJC38UtTC+kIX32b0XRoP91gox8l8NAniC2iM67uBf6Ju6KZQR9kD6jndi+4vTaMOJ+vRMmWL0b32LDSccE3omOn4a0l6zGbssjVBDbTtTlged1KDOQMtKfAk6iK/EXULg9Z0ewzNqrsG+DLwx9C553K4skyG8IDlaGbeftQSuzB36suhaict8i9kWhVpTiDdSx4nSYu7ru4hD0XTixOIz3+ak/Xopr4fzQi8CQ3WbxJVO4kyTLMqi1COk7Tw69C4s1EX5/U0r3W+aiegwv8utPTKLrT4cxNan8PUUR6/Cd1bT+ol4SVStZOpqBFrL5oAsweNjZ4cOuYZ+EsYZiaYWXVaNycZhmEYhmEYhmEYhmEYhmEYhmEYhmEYhmEYhmEYhmGE+B8VQQZW2OW0owAAAABJRU5ErkJggg==\n", + "text/latex": [ + "$$\\left [ \\rho = f_{curr 0} + f_{curr 1} + f_{curr 2} + f_{curr 3} + f_{curr 4} + f_{curr 5} + f_{curr 6} + f_{curr 7} + f_{curr 8}\\right ]$$" + ], + "text/plain": [ + "[ρ = f_curr_0 + f_curr_1 + f_curr_2 + f_curr_3 + f_curr_4 + f_curr_5 + f_curr_\n", + "6 + f_curr_7 + f_curr_8]" + ] + }, + "execution_count": 17, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "moments = [ Eq(rho, sum(f_curr)) ]\n", + "moments" + ] + }, + { + "cell_type": "code", + "execution_count": 18, + "metadata": { + "slideshow": { + "slide_type": "fragment" + } + }, + "outputs": [ + { + "data": { + "text/latex": [ + "$\\displaystyle \\begin{align} \\rho &= f_{curr 0} + f_{curr 1} + f_{curr 2} + f_{curr 3} + f_{curr 4} + f_{curr 5} + f_{curr 6} + f_{curr 7} + f_{curr 8}, \\\\ u_{0} &= \\frac{- f_{curr 0} - f_{curr 1} - f_{curr 2} + f_{curr 6} + f_{curr 7} + f_{curr 8}}{f_{curr 0} + f_{curr 1} + f_{curr 2} + f_{curr 3} + f_{curr 4} + f_{curr 5} + f_{curr 6} + f_{curr 7} + f_{curr 8}}, \\\\ u_{1} &= \\frac{- f_{curr 0} + f_{curr 2} - f_{curr 3} + f_{curr 5} - f_{curr 6} + f_{curr 8}}{f_{curr 0} + f_{curr 1} + f_{curr 2} + f_{curr 3} + f_{curr 4} + f_{curr 5} + f_{curr 6} + f_{curr 7} + f_{curr 8}} \\end{align}$" + ], + "text/plain": [ + "" + ] + }, + "execution_count": 18, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "for i, u_i in enumerate(u):\n", + " moments.append(\n", + " Eq(u_i,\n", + " sum([ (c_j*f_curr[j])[i] for j, c_j in enumerate(c) ]) / sum(f_curr)))\n", + "aligned(moments)" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "slide" + } + }, + "source": [ + "## Discrete equilibrium" + ] + }, + { + "cell_type": "code", + "execution_count": 19, + "metadata": { + "scrolled": true, + "slideshow": { + "slide_type": "fragment" + } + }, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAADm0AAAAcCAYAAAA58Bo0AAAABHNCSVQICAgIfAhkiAAAH7dJREFUeJztnXuwLUV1hz/wiiAQiRIhJAoaNaIYiYCRJOAmikZTGm7UwiIaR0VN1IgxWr7jKaNi8AWSWPgoudFoxYrEB1GCxqqriBGQh0YlouKJliiCiCAgz5s/enadObNn7+mZ6Z7px++rurXvnpk907tnrV+vtc70bhBCCCGEEEIIIYQQQgghhBBCCCGEEEIIIYQQQgghhBBCCCGEEF54N3AmcBrwtInbEgqvBC4ArgOuwvTPQZO2SISKbEWIMJAvilSRbaeN7q/ICdm7sEW2IkQYyBeFWI78YxjqPyHiRf47DPWfyAnZ+zDUf0KEgXxR2CJbEcIN8iWRKrLttNH9FTkhexe2yFaEGI78KA6OxczD/CRw+qoDtwNv9tCAx3o451icDTwTY9gPAT4G/Bi4+5SNEkEiWxEiDOSLhkcCd5m6EQkzRWwj204b3V+DtMse5VgiB2QrQoSBfNGgOMUvscY28o9hqP/iQPrnF+lfnqj/DNKXbkgv8kT9FwfSM3ukZSJ1ZCtCuEG+ZFCM4Rc99yNco/trkHbZE2t+BLJ3YY9sRYjhyI8MscQYa8CXVx2wHfeTNp9K3IFFnT2A24EnTN0QETyylbz4APATYPepG5IBhwA7gGdbHp+rL+4FvHXqRiRKKLFNrradC7neX2mXHaHokCtytXfRHdlKfijPGgflWHYoTvFHSrFNrv7hCvVfmEj//CH9E3Ny7T/piz3SCzFH/Rcm0jM7pGUiR2Qr+aGath2qSduhGMMfocQludp2LuR6f6VddoSiQ67I1d5Fd2QreaH8yA7lR3bEEmOsUZm0ufMIF7wvJqg428O5/wH4HPAD4CbgGuBi4HXAPTxcb86emL67Zsn+fTBO8M7y/T2A4zEzmr9TtvXnwBcxjjXGfWhiqv6bAtnKMHKyFYCdgGdhxPJ64EbM930RcKclnzkUeBpm0vsNI7Qxdy4EPg68ARN4tJGKLz4ZOBU4B7PE+Q7gX1Ycfy3wFeB5Fueu9wEM8/0jgDOAHwE3l6+fAR5v0ZbQWRbbTKGVbbbtk9zGhlDvb9V3Q9Uu6KZfXbQLFvWrzzgeG8qxwrV3aeP02ihbCZO+2qw8azyUYynHmpLUYpuucTyE6yOh9p8vchuflafZ8XRM3+zA+KlLpH/SvyqqA9mRkr50ITe9kFasRvHSOIydL+agZ7lpGcQzdufk2zHYCoRrLznZimra/lFNWjXpKdFzP4acdB3CvL+xaBeo3uOa1PIjSMfepY3SRltysxXNrfCL6/wI4qlLjRljBJMfbcftSpsfAg53eL4qt2Ac//2YNp8KXIC5UT8E7uXpuh/B3JxlAvPcsg1Hle//snx/BaY/TizbfG25/aMYIRubqfpvCmQrw8jJVsD8qsMO4ErgfcApwDdYfQ8+g7lPu43URgEPx9yTV1kcm4ovXlJe/3rgUtqLt2C+80W0/0pJvQ+gv++/pjzmKuB04E3Ae8rPntTSjhhYFttMoZVttu2T3MaGUO9v1XdD1S7orl+22gWL+tVnHI8N5Vjh2ru0cXptlK2ESV9tVp41LsqxlGNNRWqxTdc4HsL1kVD7zxe5jc/K09q5F8YPr8fPpE3pn/SviupAeelLV3LTC2nFahQvjcPY+WIOepablkE8Y3dOvh2DrUC49pKTraimPQ6qSasmPRV67seQk65DmPc3Fu0C1Xtck1p+BOnYu7RR2mhLbraiuRX+cZkfQTx1qTFjjLHyozUqK202sR13kzZ/B/imo3M1seuS7W/EdMq7PFzzLZhZsfdbccxZwNVsOMEfYZaerc9A3hf4PqatT3LbTCum6L+pkK0MIzZbKTDtmvX47DHlZy8H9q5svzPm1wV2lOev8gDgDoz4inG5FPg/VgcdofliQX/7PAq4PyYwmmFXvAUTELQFcfU+gH6+/5Ry32cxv+RR584t7QidVbHN2FppY9s+iW1sGEqo97fqu6FqF/TTLxvtgs190Gccjw3lWAbFzWEQojbGEtfFZisF4+ZYkHaeVTBsXPWJcizlWGOTWmzTJ44Hfz5SMExvQu0/X+Q0PoPytDZ2Av4L+C7GNl1P2pT+GZTbGVQHyktfupKjXsSST4RaR1tFQVzx4RAKxtcy6Jcv5qBnOWoZxDN2x+TbQ4nBVkBxsisK+vmGatrjEltNuiCeGANUk24ited+CuLKaaYkxPsbS/4Lqve4JLX8CNKyd2mjQdrYTmy2UqC5FXMK0n9mB+KpS40VY8B4+dEaDiZt/mF53HXAbRjDeA2LN/FdwDtazuWDh7LRMXVeVe7b2rBv/3Lfvy8579sws8MftOLad8Msd3q6ZVvn7TnV8vgx8NV/ISJbGUaotlLQX/znvwTxgoZ9B5X7Lqxtf3O5/VE9rmdLzL7ns+2vKz//2CX7Q/TFAjfB3gz74u0ME1Qt+xWMrn2wzPd3xgTlNwC/ZnmuELCNa6BfbONDK21seypCHRt8MeX97eK7oWgX2OvXjNXaBYt90GccDwXlWMtR3JyONqZmKwXui3ih2krBuDkW+M+zYu1PUI5Vp0A5Vmj4zrFcEEscX21PHx8p8PMHJ+W5G8Q8nlSZoTytzgmYPyYfifnD2g7aJ21K/1aj3E51oLbjctEXiFsvlNu3M2UdrY2CcePD3GKl+bFd88VY9SxmLQON3TaE6Nu+iMVWqu1JPU52QUE/3/BV047Zd1ST3qAgnhhjFapJdyPkemhBPDlNqCj/7c4M1Xua0HM/y4nZ3udIGw2ylXZCtZWCcZ/7CTk/Kkj7mR2Ity41w1+MsQrX+dEalUmbTUlJG28AvgDsBZwGvBf4FeDvgVfXjt2KWXZ1bJ5Qvn6tYd/DytemoO3Q8rWpze8EnoaZybvqFy7+BNgFe0e7tXy9zfL4MfDRf6EiWxlGirayb/l6ecO++baHYTRwzqOB22mZET+QWPsT/Lb93PL16IZ9OfliG5dglux++JL9Xftgme//PnAf4NPAz8rzvhzzkNvhHdo7Jl3iGugX27jWSlvbnooUx4ZVTHl/u/huitoFi33QZxwPAeVY6Y/V0kaDbKWdFG2lrzb7zrNi7U9QjhUCyrGWM0aO5YJY4ngI00eU524Q83jSh1zytAMxf0w+BaNpNkj/pH8h9V+IfddGLvoC8euFcvt2pqyjTYVre4mZPvlijHoWu5aBxm4bcvLtWGwFwrSX1PJmXzXtGPtijmrS06Oa9HL03I97UtP1NpT/+kP1Hj33MycFe5c2GmQr7aRoKz7mVsTaFxB+fgRp1KVW0SfGWIXX/GiL7YElf4MJHk4CXoGZTQrwboxhvQSzNOgdmCVK92UcZ3kpsAdmNuyhmF+r+BrNK4Yeglni9PsN+5Y5ybswBn4McA0bwvOL8l+VrZiZtE2z4+tsAf6i/P9/WhzvC9/9FxKylWHkYCtXl6/3adh338r/H4gJJHYHDsYsR32Dx3bF2p/gt+0XlK9H1ran7otduRb4HqafzmvY39YHtr5/WPl6JXAR8JDa/i8ATwau6tZ8b3SJa8A+tvGplV1seyxyGBuqhHR/bfUrVe2CxT7oOo6HgHKsNMdqaWPz95WtLJKDrfTR5jHyrFj7E5RjhYByrGZ85VguiDGOh3B8JKT+G4Mcxue+5JCnbQE+iLmnr7L8jPRP+gfh9F8ofdeVHPQF0tAL5faLhFRHGwuf9hI7ffLF2PQsBS0Djd1N5OTbMdoKhGMvqefNvmraMfbFHNWkp0c16Wb03I8bUtf1Osp/x0P1Hj33MydGe5c2ShttycFWfMytiLUvIPz8COKsS3WhT4xRZfL8aHvDxQD2w9zQL9G8QuelmEBj//L9MeX7vbs2oAc/Lq81/3cWsE/DcXcv9y8zqM+W+/erbd+x5N9a7bhdgeuBf7Ns91vL83zK8nhf+O6/kJCtDCMWWynov8zyceVnv4P5HnO2AGew8d0fV25/QPn+M/2aasXU/TmEMdp+E8Y2q4TsiwVulgGfYbcE+JxzaV7i26YPbH3/xHL/bcC3Mcva7wE8GGMDOzBxRgh0jWvAPrbxqZW2tj0msYwNrgjl/nbRr5C0C7rp1zLtguY+6DqOT41yrDDHahdIGxe/b4q2UjBcG2OxlYLxcizwn2fF3J/KsRYpUI4VAj5zLBfEGMfDcB8pcOMfofTfWOQwPteZoTxtzusxv/xb/XXUNUwbj284XvoXpv65Isb+Cym+AulLlRT0Qrl9M6HU0WwoGC8+zDVWgu75Ykx6loKWgcbuZcTg266I0VYgjzjZJQX9fMNHTXvqvhiCatKbKYgnxqiimrQhxed+CsLPaUJD+e9wZtjrl+o9m7UolfwI0rd3aaO00ZZYbKWg/3d1Pbci5r4IPT+CeOtS4DfGqOI7P1qj5cclttM8afNvy5M+ecnnvlLu/83y/YvK93deca11lhtQ07+2zt8HMyP2W8AVbCw9O+fR5XnesOTzP8XMgu3LE8vzH2dx7Lx/LmWzeC1jHbd91cTU/VdlHb/fd+rv6tNWIP3+q7JOt++6reV8O2OWMN6BEeT3ACcDX8cMgpeV+x5THn94+f4jDttYvx8h92eb7YzR9h/Sf0nu0Hxxm3XLuxdvPw18rmF7lz5o8/2TynPdDjy0tm834AflfuslwT3SNa4Bu9imSkhaCRobUoubbH03NO2Cbvq1TLuguQ+6juNToxyrndDGamnjalZ939xsZZt1yw0x28q2lvP10ea2PKtrG2PKsba1nE85lj9fnNFN/5RjbcZHjrWO+7Frar0NqX68rVPLDVP3X5V1po1tIOzvu63j+WfY91PKedrDMWPgSbXta5g2Nk3alP7ZEZL+xagHqgPFry+Qhl6EYusQpr1PrRV11pk2Psw1VoLu+WJMepaClkE4ehailkG8vp1ynAeKk9tYx51v+Khph9wXeu5nOevEHWM0oZr06rikTkg6B3nlNOtMO47FHANAWPqles9mLdJzP6sJ0d6n7r8q60gb58hWVrOO2+/qem5FzH0Ren4EcdelZviLMZrwlR+tUZm0uaWlEVX+FLM096eX7P8NzPKhV5Tv9wR+Cdy64pzfLY+x5YqW/VcCH8MsPXoZ8AHgoMr+Q8rXCxs+e1+MsQ1Z1nUrcAvtM41fAJwCfBMz2/Yai3O77qsmfPTf84GXAb8OfAN4MXCORVt8f9+UbQXi7L8jMbZyCMZengJ81KItJwN71bYdjNGsf8YMDlUuaTnfHRihPgF4evnvVswv4TwD+Efg/sBPyuNvKl93XXHOofejrz328T/XttO17X3sYDc27kNXfPuia/scwnU0L01v2wfQ7vs/K18vB75a++xNwNnAszEPwP23bcM90TWuAbvYpsrUY02dGMcGCEPLmpj6/tr4bsraBc190HUcnxrlWO0obk5HG2O3Fd/aqBxrsza35Vkx5Viu+1M5VvhxinKsZvrkWD7GrpRjmzH8Y+r+qzJ1bANxj89DSDVP2wJ8EHOvX9vhc9I/O/T3M9WBbEhVX+akoBch2DqEa+9Ta0WdqePDXGMl6J4vxqRnKWgZhKFnoWoZ+BkLYqxpQxi2AvnEyX3HAXDrGz5q2jH7Ts416dhjjCZUk477uZ9Yc5oQtKwJ5b8bqN7THz33007s9i5tlDbakvLfFV3PrcjpmR3ofh+H5EcQf13Klq4xRhOT5UfbWVxpcwtG5JfN+j0MM0P0zMq2VwM3drmwYy7GtKm6hPi/ltv2bzj+heW+N/W83p2AqzFLo67ixeV1/ge4Z89rjYGL/jsWI8jPAQ4E3olZBv7eHto7BNnKMFz13+MwvzrwZ6z+5RkbivIcswHnaGI+CN7Ixi/d7Fde64uOr1WlT3+G4n9d297VDnbGBIPf7dG2qXyxwI19zuj2azgfBs6vbbPtgyaafH9+3y5Y8pm3lPtf0eN6LukT18Cw2GbssWZqXH3fULSsjRBjiVC1C7rpV5N2QT/9ahrHp0Q5VjuKm9PSxhRtpcBP3pF7jgX+86ypfa+gf38qx1qkQDlWjjmWC3KIbQr86Dcoz4W0xucZytP2wv5XYU8uPyP9s0O5nepAuesLpKMXU9s6hG3vVUKNlQrGiw9zjZXAbb4Ykp6lomUwvZ7FomXgxrdjrWnD9LYCecXJLscB8OMbfWvaMfuOatKbKUgjxgDVpLuSW4yvv22nFwMUqN7jGj33004K9l5F2ihbsSX1vytW6TO3Ymq/KRjvmR3odh+H5EcQf11qxjQxBrjNj9aorLTZxHYWJ20+tDzZzzCGUGdufE+obDuh3DZVYHRlef1frWy7BDNjts5dMEa3zAmeD3wPE1xdCBzRcMxR5eefu6JNLy+PuZjNNzNEXPTfecB7a8d+GzjRXTOd4MpWjsQE11cs2T9HttI+4IQaWDy3PO+2yradML8McZXja1Xp05+h+N8QW7CxgwPL487o0bapfLHAjX3O6Fa8/RSLy4Db9MEymnx/b0xQfC2wS8Nnzio/89Qe13NJn7gGhsU2rrTSJiYJAVffNxQta8PF97WNI6Ddd0PWLuimX03aBf30q2kcnxLlWMqxIC9tTDHHKvCTd+SeY4H/PGtq3yvo35/KsRYpUI6VY47lgrHztCl8pMCPfsP4edDUpD4+z1CethvwviX/LsK085zy/bHlZ3LXP1BuN0d1oOXMkL5AOnqRolYUhJ3bu46XCsaLD3ONlcBtvhiSnqWiZTBt7heTloEb3461pg1p1rV94TpvDnXc61vTjtl3VJPeTEEaMQaoJt2VUJ/7KQg7pwlFy9rQ37a7MUP1HtBzP6r3GKSN0sYmUv+7YpU+cyum9puC/n3h+z4OyY8g/rrUjGliDHCbH63RY9Lms8qT7QAeX9v3gnJ7fSbqfFbpPVZdbAAPBPZt2L4z8Mby2ufW9p1fbn9AZdvuwIfY+H6/VfuM7azsU4HbgX2WtPe15fm/gln2dmp8998uwG2Y5Xur/BPw+SEN78FYtmI7C1620tx/VaYOLH6lYdthmKWfr8csX13lo+X17tfzem107c+Q/G+ILdjYwTPL417Yo21T+WKBm+BkRrfi7RdZLFis6oM+vk/Znh0YPaxyNObXP65l89Lo28rjixVtd02fuAZWxzZjaGVIv7A0xvcNScvG+L5dfk1nle+Grl3QTb+atAtW90HXcRzi0SLlWJvJMW6OWRtTzLEK+mmjcixDmzb7zLOm9r2C/v2pHGuRAuVYe9X2bWPc2MZHjuWC0PK0KXykoL9/hJYH+SbH8bnKDOVpq1jDtOX42vac9Q+U21VRHWg5M6QvkI5epKgVBWHn9iGtODaGvVSJNVaC7vkidNezbUjL5sSU+4WmZeDft2OuaUOade2+jJ03Tz3uua5px+w7qklvpiCeGEPP/SyS4nM/BeHmNCFpmf62rXqPD/Tcj+o9IG3MVRtz/Luiy7kVU/tNQf++8H0fh+RHEF9dqs4MfzHGmM/srFGZtLml4aRNHFK+nolxoA9jZj4fAfwB5sYdV/vMpeXrvsBPLa/ThT/GLCf6Bczyrz/FdOwjMU7/Y0wwUOVsjDh8HvgYsAfwKOBrwI+AuwKX1z7zEszNnM/MfhFmAPkr4JWV444BvkTzMufPAF6PufnnlOeos864v/7hu//2xiwtW++PK4FHu/0qrYxlK2dht4yubKW5/0Lis5ilur+OCSQejEmqbsYEjvW2nwE8CXgs8B0P7enanyH5n29beAzGXz7R47Ox+SKYNh9T/n8eOBxeacfVwEuXfHYf4JMN51vWB318H8y4+XvAqzG/knM+Zhn4rZj+fA4mQJkz/5Wp25a02wd94hpYHduMoZW2MckYjPF9Q9KyMb6vbRwBy303VO2C/vrVpF3z8y3Tr67jOMSjRcqxNgjR3qWNq7+vcqwNlGPZabPPPCtm31OO5ZYUcywYP7bxkWO5ILQ8LTYfCS0P8k2O47PytOHkrH+g3K6K6kCbkb4skopeSCs2CO3vz2Mwhr2ExJj5InTXM2nZBrHkfiFqGfj37ZDqaqFpt8a+sMcC1zXtmH1HNWm3pFiTTiUuCS2m8I3+tq2/bXdF9Z5F9NyP6j0gbcxVG3PLj8Dt3IqY/Sbk/Ajiq0vBeDHGmM/stLKdxZU2v4xxqLtilpD9Yfn+f8uL79pwnp0wHfTnXRtgyUGY2dGXlNe5Dfg5cAFmZmrT7N9dgVMwSyzfiAmInoeZ2XoH5rtXsZ2VfRhm9uxLlrR1jY2Z08v+1a/tG9/9tx/me9WXO38dxm7GZAxbqbNsFrxsZVj/2VIwbMb+y4ALMQJ6M/A94DTggCXH74IR6vN6Xq+Nrv0Zkv8NsYU2O7gbJgD8eI92TemLBf3ts61d60s+tyemvw+vbGvrgz6+P+fuwNsxvnMLJrD5BPCIhmMvBq5j83LivukT18Dq2Ma3Vob0C0swztgQkpaNPRau0r9VvrtGmNpl07b1hs80aRe061fXcRzi0SLlWBuskWfcHLM2pphjFfTTRuVYdtrsM8+a2vcK+vencqxFCpRj1Rk7tvGRY7kgpDxtKh8p6O8fIeVBY5Dj+LyG8jRb1jDfr77SZq76B8rt6qgO1K1t6w2fSV1fUtGLFLWiIJ/cHsaND6f+vgXjaxn0yxehu55JyzaIJfdbIzwtA/++HXNNG9Ksa/clpDjZloL+vuG6ph2z76gmvZmCeGIMPfezSIrP/RSEm9OEpGWh5DMx5r+2bVtv+IzqPYvEnh+B6j11pI15amNI+aUtBcPGAZdzK6b2m4L+feHzPg7JjyDOupRN29YbPtMnxhjzmZ01KittNrGdzZM274QxqAtXfWgJ7wfe1uNzoTB38CNr2/8O+Fbl/ZvK4+4zUrtiILTJJWOzTFBlK3a4KL6OzSsx7f7dqRtCOv7XZgd/TXMQZkNuvngkpiiwc2VbCH2wF+ZXJk4a8ZpD4hqYLraxjUlSIhUt68Mq/QvBd8eiSbvAfR/EpkXKsfJG2riIbMWOGHMsCCfPSsX3lGO5I9QcC8aPbWLNsVzQJU8LxT5CJtaxaggxfueU87Su5Kx/oNzOJaoDGVLWl5z1QlrhlhhjhyHE+n3HyBelZeOi3G8YqdTV+qK6dn9iHQdc1bRT8R3VpN0Rak1acYk9eu5ng9i0rA+KAQyq9zQTc34EqvcMRdq4iGzFjhhzJOVHi6y6j0PyI8jLl8aKMfqyRmXSZr2RTTwI2A24qMfFTsMs0Ro7O2rvd6pt2wp8FTNzVhhuwQSjR9e2H41ZbjZXZCvp8g7g+5ilo6cmB//bDRPMnYFZprsrufniY4D3Yn5VYk4IfXAEcCvm1yfGYkhcA9PHNm0xSUrkoGV9CMF3x6JJu8B9H8SmRVPrkCuUY/VD2riIbCVtQsmzcvA95VjdCDXHgvFjm9hzLBfY5Gmh2IcQQ0k5T+uK9M+g3M4vOfVfyvoivZBWiLwYI1+Ulk2Dcr9+5FBX64NsJV1c1bRz8B3VpLsRak1acUl39NxPWlrWlRD8dkxU72lmah1yheo9/ZA2LiJbSRflR/YMzY8gL18aK8ZwwhaLYw4pXy/ucf7zgcuAQzFLv8bG1Zhf39i3tv2ewJWV9weO1qK4eDvwQYwdnItZ+nc/TMCZK7KV5ewB3K/y/gDgYOAazIAdOr8Eng4cBewO3DBtc6L1P1s7OAB4D7Ct53Vy8sWdgccDj6xtD6EPzsQsFT8mQ+IamC62sY1JUiNWLfNJCL47Bsu0C9z3QWxapBxLSBs3I1tZTuw5FoSVZ8Xqe8qx3BNyjgXjxzax5lgu6JKnhWIfYnpiHp9Tz9O6krP+gXK7scil/1LXl5z1QlohuhJzrATj5YvSsnFR7jecWOtqPpGtNBP7OABua9qx+o5q0u4JuSatuMQePfcTl5b5IgS/HQvVe5YTc34Eqve4QNq4GdnKcmLPkZQfGWzu4wEMy48gH18aM8bwxnbgzZX3p2J++eARPc/325hZrLFyHsYBqlwGnDhBW2Lk+cA6cDNmhnt9OfRUiXEJ6qmZYfqt/m/bdE2Knhj9b4bswDVbgROmbkRADI1rYLrYJteYJEYtG4riiPS1SzlWnnrmEmmjsGGGYmvXxOh7M2QHrkk9TulKzDmWCxTXuCOXsX5GvLos/dtM7voH0kBX5KJ/q0hdX3LXC2mFO3LQixnxxkqQtp5Jy6RlQ4mxruaCHLTbJTPiHgd8EKPvzNB9dE3KMUYfYo5Lco0pYtSyoSgGSF+79NxPnnrmEmmjsGGGYusqsfrNDN1Hl8QQY6wBX151wHY2T9p0wTMxxhYjx2KW1D0eM/P2ZOAXwP5TNkoEyR6YWe8HY4T0peX/7z1lo4QQWbMncArmVyWEW6aIbRSTpI3iiA2kXXYoxxI5IG0UQoSG4hR/xBrbKK4Zhsb6eJD++SNW/QNp4BCkfxtIX+yJVS+kFcOQXsSD9MwOaZnIAWm3EMIlijH8oed+hGsUA2wg7bIj1vwIpGfCHmmjEMIlscQYa7RM2nw38CngfcAzHF74iQ7PNTaxzsoW4zJDs+CFEGHxKOCuUzciYaaIbRSTpMsMxRFzpF32KMcSqTND2iiECAvFKX6JNbZRXNOfGRrrY0H655dY9Q+kgX2ZIf2bI33pRqx6Ia3ozwzpRSxIz+yRlonUmSHtFkK4QzGGX/Tcj3DJDMUAc6Rd9sSaH4H0TNgxQ9oohHBH6DHGcZh5mP8BnD5xW4QQQgghhBBCCCGEEEIIIYQQQgghhBBCCCGEEEIIIYQQQgghhBBCCCGEEEIIIYQQQgghhBBCCCGEEEIIIYQQQgghhBBCCCGEEEIIIYQQQgghhBBCCCGEEEIIIYQQQgghhBBCCCGEEEIIIYQQQgghhBBCCCGEEEIIIYQQQgghhBBCCCGEEEIIIYQQQgghhBAiSP4f570Ehom6B9IAAAAASUVORK5CYII=\n", + "text/latex": [ + "$$\\left [ \\frac{\\rho \\left(- \\frac{3 u_{0}^{2}}{2} - 3 u_{0} - \\frac{3 u_{1}^{2}}{2} - 3 u_{1} + \\frac{9 \\left(- u_{0} - u_{1}\\right)^{2}}{2} + 1\\right)}{36}, \\quad \\frac{\\rho \\left(3 u_{0}^{2} - 3 u_{0} - \\frac{3 u_{1}^{2}}{2} + 1\\right)}{9}, \\quad \\frac{\\rho \\left(- \\frac{3 u_{0}^{2}}{2} - 3 u_{0} - \\frac{3 u_{1}^{2}}{2} + 3 u_{1} + \\frac{9 \\left(- u_{0} + u_{1}\\right)^{2}}{2} + 1\\right)}{36}, \\quad \\frac{\\rho \\left(- \\frac{3 u_{0}^{2}}{2} + 3 u_{1}^{2} - 3 u_{1} + 1\\right)}{9}, \\quad \\frac{4 \\rho \\left(- \\frac{3 u_{0}^{2}}{2} - \\frac{3 u_{1}^{2}}{2} + 1\\right)}{9}, \\quad \\frac{\\rho \\left(- \\frac{3 u_{0}^{2}}{2} + 3 u_{1}^{2} + 3 u_{1} + 1\\right)}{9}, \\quad \\frac{\\rho \\left(- \\frac{3 u_{0}^{2}}{2} + 3 u_{0} - \\frac{3 u_{1}^{2}}{2} - 3 u_{1} + \\frac{9 \\left(u_{0} - u_{1}\\right)^{2}}{2} + 1\\right)}{36}, \\quad \\frac{\\rho \\left(3 u_{0}^{2} + 3 u_{0} - \\frac{3 u_{1}^{2}}{2} + 1\\right)}{9}, \\quad \\frac{\\rho \\left(- \\frac{3 u_{0}^{2}}{2} + 3 u_{0} - \\frac{3 u_{1}^{2}}{2} + 3 u_{1} + \\frac{9 \\left(u_{0} + u_{1}\\right)^{2}}{2} + 1\\right)}{36}\\right ]$$" + ], + "text/plain": [ + "⎡ ⎛ 2 2 2 ⎞ ⎛ \n", + "⎢ ⎜ 3⋅u₀ 3⋅u₁ 9⋅(-u₀ - u₁) ⎟ ⎜ 2 3⋅u₁\n", + "⎢ρ⋅⎜- ───── - 3⋅u₀ - ───── - 3⋅u₁ + ───────────── + 1⎟ ρ⋅⎜3⋅u₀ - 3⋅u₀ - ────\n", + "⎢ ⎝ 2 2 2 ⎠ ⎝ 2 \n", + "⎢─────────────────────────────────────────────────────, ──────────────────────\n", + "⎣ 36 9 \n", + "\n", + "2 ⎞ ⎛ 2 2 2 ⎞ ⎛ 2 \n", + " ⎟ ⎜ 3⋅u₀ 3⋅u₁ 9⋅(-u₀ + u₁) ⎟ ⎜ 3⋅u₀ \n", + "─ + 1⎟ ρ⋅⎜- ───── - 3⋅u₀ - ───── + 3⋅u₁ + ───────────── + 1⎟ ρ⋅⎜- ───── + 3⋅\n", + " ⎠ ⎝ 2 2 2 ⎠ ⎝ 2 \n", + "──────, ─────────────────────────────────────────────────────, ───────────────\n", + " 36 9\n", + "\n", + " ⎞ ⎛ 2 2 ⎞ ⎛ 2 ⎞ \n", + " 2 ⎟ ⎜ 3⋅u₀ 3⋅u₁ ⎟ ⎜ 3⋅u₀ 2 ⎟ \n", + "u₁ - 3⋅u₁ + 1⎟ 4⋅ρ⋅⎜- ───── - ───── + 1⎟ ρ⋅⎜- ───── + 3⋅u₁ + 3⋅u₁ + 1⎟ ρ⋅\n", + " ⎠ ⎝ 2 2 ⎠ ⎝ 2 ⎠ \n", + "───────────────, ─────────────────────────, ──────────────────────────────, ──\n", + " 9 9 \n", + "\n", + "⎛ 2 2 2 ⎞ ⎛ 2 \n", + "⎜ 3⋅u₀ 3⋅u₁ 9⋅(u₀ - u₁) ⎟ ⎜ 2 3⋅u₁ \n", + "⎜- ───── + 3⋅u₀ - ───── - 3⋅u₁ + ──────────── + 1⎟ ρ⋅⎜3⋅u₀ + 3⋅u₀ - ───── + \n", + "⎝ 2 2 2 ⎠ ⎝ 2 \n", + "──────────────────────────────────────────────────, ──────────────────────────\n", + " 36 9 \n", + "\n", + " ⎞ ⎛ 2 2 2 ⎞⎤\n", + " ⎟ ⎜ 3⋅u₀ 3⋅u₁ 9⋅(u₀ + u₁) ⎟⎥\n", + "1⎟ ρ⋅⎜- ───── + 3⋅u₀ - ───── + 3⋅u₁ + ──────────── + 1⎟⎥\n", + " ⎠ ⎝ 2 2 2 ⎠⎥\n", + "──, ────────────────────────────────────────────────────⎥\n", + " 36 ⎦" + ] + }, + "execution_count": 19, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "f_eq = []\n", + "\n", + "for i, c_i in enumerate(c):\n", + " f_eq_i = w[i] * rho * ( 1\n", + " + c_i.dot(u) / cs**2\n", + " + c_i.dot(u)**2 / (2*cs**4)\n", + " - u.dot(u) / (2*cs**2) )\n", + " f_eq.append(f_eq_i)\n", + "\n", + "f_eq" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "slide" + } + }, + "source": [ + "## BGK Collision" + ] + }, + { + "cell_type": "code", + "execution_count": 20, + "metadata": { + "slideshow": { + "slide_type": "fragment" + } + }, + "outputs": [], + "source": [ + "tau = symbols('tau')" + ] + }, + { + "cell_type": "code", + "execution_count": 21, + "metadata": { + "slideshow": { + "slide_type": "fragment" + } + }, + "outputs": [ + { + "data": { + "text/latex": [ + "$\\displaystyle \\begin{align} f_{next 0} &= f_{curr 0} + \\frac{- f_{curr 0} + \\frac{\\rho \\left(- \\frac{3 u_{0}^{2}}{2} - 3 u_{0} - \\frac{3 u_{1}^{2}}{2} - 3 u_{1} + \\frac{9 \\left(- u_{0} - u_{1}\\right)^{2}}{2} + 1\\right)}{36}}{\\tau}, \\\\ f_{next 1} &= f_{curr 1} + \\frac{- f_{curr 1} + \\frac{\\rho \\left(3 u_{0}^{2} - 3 u_{0} - \\frac{3 u_{1}^{2}}{2} + 1\\right)}{9}}{\\tau}, \\\\ f_{next 2} &= f_{curr 2} + \\frac{- f_{curr 2} + \\frac{\\rho \\left(- \\frac{3 u_{0}^{2}}{2} - 3 u_{0} - \\frac{3 u_{1}^{2}}{2} + 3 u_{1} + \\frac{9 \\left(- u_{0} + u_{1}\\right)^{2}}{2} + 1\\right)}{36}}{\\tau}, \\\\ f_{next 3} &= f_{curr 3} + \\frac{- f_{curr 3} + \\frac{\\rho \\left(- \\frac{3 u_{0}^{2}}{2} + 3 u_{1}^{2} - 3 u_{1} + 1\\right)}{9}}{\\tau}, \\\\ f_{next 4} &= f_{curr 4} + \\frac{- f_{curr 4} + \\frac{4 \\rho \\left(- \\frac{3 u_{0}^{2}}{2} - \\frac{3 u_{1}^{2}}{2} + 1\\right)}{9}}{\\tau}, \\\\ f_{next 5} &= f_{curr 5} + \\frac{- f_{curr 5} + \\frac{\\rho \\left(- \\frac{3 u_{0}^{2}}{2} + 3 u_{1}^{2} + 3 u_{1} + 1\\right)}{9}}{\\tau}, \\\\ f_{next 6} &= f_{curr 6} + \\frac{- f_{curr 6} + \\frac{\\rho \\left(- \\frac{3 u_{0}^{2}}{2} + 3 u_{0} - \\frac{3 u_{1}^{2}}{2} - 3 u_{1} + \\frac{9 \\left(u_{0} - u_{1}\\right)^{2}}{2} + 1\\right)}{36}}{\\tau}, \\\\ f_{next 7} &= f_{curr 7} + \\frac{- f_{curr 7} + \\frac{\\rho \\left(3 u_{0}^{2} + 3 u_{0} - \\frac{3 u_{1}^{2}}{2} + 1\\right)}{9}}{\\tau}, \\\\ f_{next 8} &= f_{curr 8} + \\frac{- f_{curr 8} + \\frac{\\rho \\left(- \\frac{3 u_{0}^{2}}{2} + 3 u_{0} - \\frac{3 u_{1}^{2}}{2} + 3 u_{1} + \\frac{9 \\left(u_{0} + u_{1}\\right)^{2}}{2} + 1\\right)}{36}}{\\tau} \\end{align}$" + ], + "text/plain": [ + "" + ] + }, + "execution_count": 21, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "collide = [ Eq(f_next[i],\n", + " f_curr[i] + 1/tau * ( f_eq_i - f_curr[i] ))\n", + " for i, f_eq_i in enumerate(f_eq) ]\n", + "aligned(collide)" + ] + }, + { + "cell_type": "code", + "execution_count": 22, + "metadata": { + "slideshow": { + "slide_type": "subslide" + } + }, + "outputs": [], + "source": [ + "def bgk_collide(m):\n", + " return list(map(lambda expr: expr.subs(eq2assgn(m)), collide))" + ] + }, + { + "cell_type": "code", + "execution_count": 23, + "metadata": { + "slideshow": { + "slide_type": "fragment" + } + }, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAADqwAAAAqCAYAAAB8rpsbAAAABHNCSVQICAgIfAhkiAAAIABJREFUeJztnXm4JEWVt99e2ARpWWQbGltABrQbENlk84qgbMrigjoipaKoLMp8KIICjTOIDjijiOCoM1/iMCIoKosiKsgmqwg2CAoqVxsFZFEBWRv6++NkfjdvVlZVVmZGZkTV732efO6tXH8nIjJPRmTECRBCCCGEEEIIIYQQQgghhBBCCCGEEEIIIYQQQgghhBBCCCGEEEIID5kNbNW2CMHMtgUIIYQQQgghhBBCCCGEA+YClwO3A78A9mtVjRBCCCGEEEKI0FC9UgghhBChMaNtAUIIIYQQIjhe0bYAIepiGeC/gfXbFjLm7APs1LYIIYQQQgghhBBCCCGEcMDawGbx/2sAi4HntSdHCCGEEEIIIURgqF4phBBCiNA4HptQSAghhBBCiKIcCezVtggh6uCLwO5tixhzdgdOcnyN5R2fXwzHQcDvgCXAGS1rKUPo+iF8G0LXD+HbELp+IZpiFO6V0G0IXT+Eb4P0i7pYBkXhFcVZBbgf2KDieUbhGRC6DaHrh/BtkP5yfAv45z7bFwHrNaRFiEHIb04Rug2h64fwbZD+cgzym0K0gfyjEbp+CN8G6W+fNmwo4htVrxQ+Ib9phK4fwrdB+tvHV78p2mFL4KttixBCCCGEEEExE7gI2KRtIUJUYX/g7LZFjDnrAz8FlnV4jf2AXR2eXwzHxliD1BuxCKArtitnaELXD+HbELp+CN+G0PUL0RSjcK+EbkPo+iF8G6Rf1MnzsSi8Yrw5BOuM90i8XAvsmbPfycCZFa81Cs+A0G0IXT+Eb4P0l2dT4GFgTs62LYE7UCAG4R75zeEI3YbQ9UP4Nkh/efr5TSHqRv6xOKHrh/BtkP72acuGQb5R9UrRFPKbxQldP4Rvg/S3j69+U7TLQuADbYsQQgghhBBBsQlwC7BC20KEKMMc4B5gw7aFjDmXArs5PP+rgBMcnl8Mz0eAm9oWUYHQ9YOfNkRY41QRQtcP4dvgo34hfMTHeyUi/OfVMPioP0J50DYR8nkhswdwVNsiRKvsjZWDlwAbAScCz2CdERKeB/wV2LHitUbhGRC6DT7qjwjfl0eE7QsjwtF/E9axM81qwO3Ads3LEWOI/OZwhG6Dj/oj5DfbJiIc/Xl+UwgXjLt/jAjnuVAHPtoQMT554KP+iHDeT3r5RtUrRZPIb47PMxv8tCFifPLAR/0R4ftN0T4rAL/GfIkQQgghhBBFOQ2NBROB8hk0u2rbvBW4weH5XwhcBSzn8BpiOO4ElqaW89qVMzSh6wd/bYgo1sAZun4I3wZf9QvhG77eKxHhP6+K4qv+COVB20TI54XOWcAObYsQXvEwcHDq95vjdVVmmBiFZ0DoNviqPyJ8Xx4Rti+MCEf/8cBPU7+XA64EDmhBixAJ8pv5hG6Dr/oj5DfbJiIc/Vm/KUSTjJN/jAjnuVAVX22IGI888FV/RDjvJ3m+UfVK4QPym934qn8YfLUhYjzywFf9EWH7TeEP7wcubluEEEIIIYQIinnAY8CLW9YhxFDMAR4Bdm5byBgzA7gDeIfDa5yJDYoV/rAGFi3r48BawMrtyhma0PWDvzZEFGvgDF0/hG+Dr/qF8A1f75WI8J9XRfFVf4TyoG0i5PNCZz3gF8AyPbYfDdyI1fsfAC4E5jcjTTTMLOBtwNPAgtT6zwM/qnjuUXgGhG6Dr/ojwvflEWH7wohw9O+GPaNWwNpDz2a4GQmEqBP5zf6EboOv+iPkN9smIhz9ab8pRFOMo3+MCOe5UBVfbYgYjzzwVX9EOO8nWd+oeqVoG/nN3viqfxh8tSFiPPLAV/0R4fpN4RfLA38GXtO2ECGEEEIIERTfBr7UxIVmNnERMRa8D/gL8JO2hYwxewPrYA8QF2yOzbhzrqPzi3I8AqyPRTO7L/5dN3tgjV93AR+s+dxN6L8Aez59y8G5wa0Nc4HLgduxAQT71XjuhCbyYB5wGWbHL4HVaz6/SxsWALeklieAfWo8PzSTB4cAt8XLf2EfvoQIjdB9HshnFOEaTP9twHEOzu/ahklgEeYzXNSNQvfb8nn+8gfs3ju0x/YJ4HRgOyxQ1RLgx8CqTYgTjbAAi6D3FJbX+wC3prbPA/5U8Rry5YMJuf4Ko/E+Mo9w66/gvg7rgy//ExZgYR1ge2B/zMbE5gUI4R75zWLIb/ZHfnMw8puDGcZvCuGaUfCPofvGUfAtobcPg9s24tDbh8G9DcP6RtUrRVvIbw5GfnMw8pv9kd8cjOqUYfMkFnjj6LaFCCGEEEKIoDgL6GBBcVphZeBUrEL4NLAUOKotMR6i9OlmEXBGg9dTHnTzfeAbDs9/NvAJh+cXxkKsPPdbJlL7bxmve0FN1882PC2HNS7PBZ6HNbCtXYPuBNf6ZwGvBl5P/45LCymnH+q1Iav/H4DN4v/XABZj+ZDHMdgHjWR5Bns+ptftmHNcE3lwBbBT/HsOVq6ylNUPbvMg/Xsl4EFgxR7H+poHawK/xSLrzcBmRNurpmsJUZWFFH/++uLzfH1eDeMzFuKv31459f+1WNCULL7mwSysfrJSgWMX4m8euPTb8nl+sy324Tkvz7OsBDyLvWeL0WBZYEPsPj0Je+9Nz6J7CfCfOcctJDxfPozmLK79SJH6K/jrR4q+j/jsy4v4QfDXFxatw/qqv4gvf0msQTN9izaR35TfLIP8pn9+R35TiHop4x8X0t530bK+Efx9LhT1LQvx178XaR8Gf/OgaBvxQvx8Pwn9u658owgJ+c32n9nym+3ngfym/KaozjZYHr1swH5fw2Zj7dXnTf2xu1GadKM06UZp0o3SJB+lSzdKk26UJt2EkiavwLS9p20hBVkO+DtwZFsCLsQS7HvAv2KVvk3aEuMhSp/pbIKlR90zzvVDeTCd1bCHcGfAfjOAdwPXAY8CjwM3A4fTf+ahNbAZdDasKlQMZHVg4wFLuoH0IMwJp5kLnI/l7S+BF8f/vzDe/grgvNT+FwKfAa4HDs78Pj4+V8KJwHtr0D2sfvrY0E//wfE+E/TvuFRW/zA2VNGfsAhYr4eOVbF7NFnOw17S0utWqKCfkjZ8CJt9axBl9Q9jQ9U8+Cf6BwZwnQdl7+OjgN9jDdjLYPmxVR87hGiSYZ6/rp9XRX1eCM8r6O8zfPfbYGl4E/kfVn3Og0mKDVh1nQdlfYZrv+26DMnnVed3DK5XgnU8WYrNQJDHmtiA1lPj36th+f8d4DfYzEV/A67GGutmllYsXPFj4Kup3/8LnJuzX4i+3Pf63wSDB96E4Muh9/uIr768qB+sYoNrX160DhuyL0864KzZwzYh2kB+Mx/5TfnNqjbIb8pvirAp4h+b/i5ah28Ef58LRX2L7/4d+rcPg9++ZZLBbcS+vp+E/l1XvlGEjPxmN/Kb8psJ8puDbZDfDJc3AV8ArsJm0F2KzXhVlBlY0INP9dlnS+A54J/77KP+2N0oTbpRmnSjNOlGaZKP0qUbpUk3SpNuQkqT7wD30r9+k+1PB/Z+fikWROkJ4GHs/f54rL9dP3bE3v3vBZ6K//4Q2KOA3h9i9ch+5Okditk56zbGIuVcAuxZ9sQjjNKnm93jv9c3dD3lQTd7YQ0HVwzY70zgACxa0jnYyPhdgM9j0b7ejD3Us7wZuBvrWCzc8mC8FGVz4JbU72Wx2XYPB36CRW97ElgFeCDeZwFwW+qY+cDlWCMTwMdSv9/EVIMWwD1YlMGquofVP7uPDf30F6WsfqieB0X1b4nd54t76Hg4XhIejX8Pum9d58E+sZbzscb9C4HjatQ/jA1V8+AtQNRHh+s8KHsfgwUo+AMWqfFs4MYBmoRoimGev66fV0V9XgjPq0E+w3e/fT3wUuD0zPUSfM6Dj2Hv5M8Bn8M6FOThOg/K+gzXfruJMiSfV42LgXfQ/50HrHzfggVDymNvbBDqd+LfbwbOwBrGfoLl0ZrAflgnnN3pXR8V7TATi6qdcDP5g5lD9OUh1P8G4bsvh/7vI7768qJ+sIoNTdT/EvrVYUP25QuwGcHvH6BViCaR36xmg/ym/KYr/YNsSCO/KUT9FPGPTX8XrcM3gt/PhYR+vsV3/z6ofRj89i1F2oh9fT8Zhe+68o0iVOQ3y+uX35TfrKpfflN+s00+gc32/Bj2LN14yOOXYvfQW7DZfPP4FDYY9owe29UfuxulSTdKk26UJt0oTfJRunSjNOlGadJNaGlyElY/O5zewUSy/ekAjgB+DvwIG+O1IrAtNjj3ffH/eXXVTwD/gtWVLsL65K0OvBwL+Pv9AXqvwAYB/wPwxyH0DkXeDBY7x3/Py9kmRi99OlilZaLCOXbAKrz31qCnCMqDbnbCZqe5u88++2CDVe8GXoZF4PoQ1rDxXeCNwIF9jr2qgj7hjmzD1L7ADVijFFi52IDpjVybMtUw9XysMevfe/yekXPNOjuRF9H/FLAR+TYM0t8EVfKgqP7VgK9hs0/V3YnfdR7Mxp5vhwNbY5H89m3Bhqp5MAd4JTaIo25c38erYC/r84B1sQbpiRr1C9EUrp9Xrn0eyGcUtWEbrCK8OfbMqhPXebA95uvegH0QWlCvfOc+w7Xfdl2G5POqcxnwavpHaDuZqYBHz/bYZ1/gIeDK+Ped2H2xLjZz0dHAu7HGxcVYfXS/itpFeT6NRd2bhz23TsLunXTHkEuwKIWrV7iOfHkY9dc6fHnb7yNl88D3+mue5qbrsD748p2AH1Q1RIgKyG8WR35TfrNN/Xma5TeFcIdP/rHsc80X3zjOviXB1/ZhaL+NWN915RvFaCC/WQz5TfnNqshvym82QYdqfX+PwMrgysAHSp7jZ1g5WDdn20bYhDbnYjN25aH+2N2MWppA9XRRmnSjNOlGaZLPqKVLB6VJlg5Kkywdxi9NbgB+BRwMzOqxT7Y/Hdh74LZYP7qPAYcBW2GDXtfB+thleTM2WPXHwPrAu7D61PviYz9eQO8v4r879NknT+9QpAesvhErFF+Mf385/r2U4aO2jCJKn95sg91crlEe9GZ7pkcnyyPp/PtZpkcdewY4Nv7/sJzjZmMf6ntFYRPtMQNrZErnzQLM4aWZz/TysWXq93zgGqYaXbO//wjMTR27LhY5rQ6K6k905dkwSL9rquZBEf3LYZEpToq31UkTeXAPNmX877HnzUVYg2pdNJEHYAP3f4BFJqyTJu7jXYDfYhEanwC+hzWWCxESTTyvXPo8kM8Y1m8/gg2c262S6uk0kQdJmbkXi1L1ijqExzThM1z67SbKkHxedW7G2op6NUZ9Fngn8Bp6R36egzUaXsjUgNbL4t/PZfa9D/hS/P9EKcWiDtYCzgJ+DVyKNV7uzvQO+7dikQDfWvIa8uXh1F8h/PeRUa2/5mlusg7rgy9fAbPtKxXsEKIq8pvFkN+U36yK/Kb8pggLn/xj6L4x0TaOviWNj+3Dyba22oj1XVe+UYwO8puDkd+U36yK/Kb8Zij8BLiLam1eSf7vmLPt3VhZOidnm/pjd6M06UZp0o3SpBulST5Kl26UJt0oTboJOU2+AayHvWtnyetPB72/TZ0b/31JZv1M4DPA48DbgUdzjn2mgNbkHXK7Htvz9K6O9f1bOmB5ClgebCBawr3ACcAHseg5J8brl2IvxOOO0iefFbGR26VHTQ+B8iCfZYANsY7F/Vgr/vu7nG3Jui2AFwB/TW3bFMvnOypoFG5YikVVSHM/UxHzZmL3yqpMRcjaFmvcS+6Z+cCi1PHZ3zdgM/LOxQY670u+Ey1DUf0P0duGdw3Q75qqeTBI/wwgwhqW/2dIbZ0C+zSRBzdiESlXwxo5X0WxqeGL6B/GhqplaH/gtIKaEjoF9mniPl6MBR5YHnsJncBe4IUIiSaeV2V9XqdmG9rwGUVoIg/mYFFdH8CeWa8F/qOAtk7NNpTVv2J8jkeBlbAK+7nURxM+w6XfbqIMyedV527g71hj1PmZbadi70SvBm7vc449sXv52wWvmTSSLSkuU9RMp+B+JwBfAM6g9+y6vfDZlxelqbqHK5rIg7LvI52C+7nOg7J+EPyp/yUMW4ftFNjHB1/+HuC6eBGiLToF95PflN+U3+xNpwH9eZrlN4VwR6fgfk34xzZ8Y6dG/W34lqK4zoOy7cPgj29x2Uas77ryjWJ06BTcT35zsH75TfnNXshvym+OE3fGf+fnbNsF8yF5+aj+2N0oTbpRmnSjNOlGaZKP0qUbpUk3SpNuQk6Tn8Z/dwUuyWwbtj/d6+O/2e9a2wEvBr4F/CU+73xs4OsNwLUFz78YG3y6UY/teXpXAj6Z+j0POBAL5HNRav0DsZ5pA1avwSJzfRQbmLawoFBf+TA28K4otwDf7bN91NKnLjaI/9YZJa0XyoN85mHTRg/Kgwfjvy/O2bZ+6v+NmV5BTR5CfywjTjROhEXEug1rTDoMi8J4EVZW7oqXpDP4AixiIz1+LwGOiNfNBD6P2/s9olv/1fS2YZB+sOnON8MaOu/BpkEv6ozLEFE8Dwbp3x7rLLMIi1oHcAAWTdMVEfXmwbPAUVhEuBnAFdQ7eCePiPryAKxxdwvgRw41p4mo9z6+Lj7+Zuzl8lLgAlfihWiQiHqfV037PJDPyP5eBTgPC8gyE/MX6YqsCyLqy4M1sQr6DOz9/CvYh0qXRNTrM5r22xH1liH5vOosBf4A/GNm/enAO7Dny8NMBUR6LF7S7IsNei3y7jQbm7EVbEYj4TeXYB3518UihlclQr7cx/prlTzw5X1kVOuveZrbrsNGNOvLn46vIUQIyG92EyG/Kb9ZjQj5TflNETpN+MfQfeM4+xbf24fzNLfdRhyh77ryjWKUkd+cToT8pvxmNSLkN+U3R5N747/rZ9aviM36ewf27TaL+mN3ozTpRmnSjdKkG6VJPkqXbpQm3ShNugk5TZL6zE452wb1pzsSGxA6B9gS2AGrr346s99W8d/7gZ9j7/lprgTehA0a7cez8T7Zd8h+eieZnh8HYQNWvwGcMuB6EItdilVsQmeSwVPNppeowDlHKX0SOphNEyWPf218/HE16RmE8qCb18THf2zAfm+P9/sN9gE+YTbWgJXcC7tnjjsmXr9KSX1CCCGEEEIIIcLhUuAXmXW92lIWZvZbHouE/c2C1zolPs/3SmoVQgghhBBCCCGEEEIIIYQQQog26FCt72+aifhcZ5U49nGmZvRK2Cg+3w/7HKf+2N2MYppAtXRRmnSjNOlGaZLPKKZLB6VJlg5KkywdxjdNngDuy6wr0p/uPqb3ybsYC/KT5aR4+xIsWM1rsIGuL8Mmi1gKXF5Q663A33LWF+3/d1p8vV0KXo8D4wMOL3rAmBF6+kxS/yDefeN9P+rg+nkVL+VBN3vH+354wH4zge/H+94HfBn4HBZx6wngznjbazPHfTZev8yA808OsKVI/gohhBA+Mol8nBBCiPHh29gMVWV4A+YL315g38Pjfe9gelClfkwinyyEEELUzSTyr0IIIcaDSeTzhBBCiLqZRP5VCCHEeDBJ/X1/00xQ3lc+hM3AleaV8fnO6XOc+mN3E3qaQP3pojTpRmnSjdIkn9DTZRKlSZZJlCZZJlGapPkjNpg0zTD96dbExuj9GvgTsEVm+7/F53oW2CyzbQVgcbz9lQWudWOO1mH0XhXv98JeO8zO/E6MubmAOBecAmwC7JlZ/0HgI8DawC+xgXlXNSsNaD990pRJk88BL8is2xwb8Hgm9rBIc0sBHSvGf58ssO9vC+6X8KecdcqDbormwXPYw+NDwAHx8gw2bfaB2Aj3lwB/zjn/c/G+/agjf4UQQggfkY8TQggxTjzJVD1zWPYFnmbwjKmHAJ8HbscivT1c8PzyyUIIIUT9yL8KIYQYF+TzhBBCiPqRfxVCCDEuuOj7Wxd533efiP8u3+c49cfuxpc02QlLj1dgafJm4FsFj607XXxJk6OB/YB/BJ4CrovX3Vbg2FFNk0OAg4F58e9fAv/K4P4aMLppkuUY4ETgi8ChA/Yd5WfKQuD4zLr7gbUGHDfKaQL2fP00sAc2Q+VvsHJyRZ9jRjlNJoEX5aw/HXve9GKU02QWdv+8Aysv9wL/G6/LG+gJNmj0icy6ov3pwO7N7wA/xyYk/BowP7X9L/Hf3wG/yBz7BHAJ8B5ga+DaAdd6ErNxOcy3DqN3BrAp1s7zwIDr/H+uwgamPb/oATVzBXBCZt3+2EC592KDWU8FHgPWG3CuD2MFoeiyTwF9TaZPr9ksl6F8muTRwUY1T5Q4FuBt8fFHlDx+WJQH3fxTfPzBJY+HqQfj43Tb/VXsgSOEEEIIIYQQYvT5Gla/HZZZwIPAxQP2+zBWh70VWKPEdYQQQgghhBBCCCGEEEIIIYQQom06VOv7m2aC8jOsLsYGLaRZJz7f1X2Oa3rMQF6fbN/6Y/uSJrtjAw/3w+x5U8XrdCifLj6kCcClwLuwQTsLsME89wGrlrxOh/DTZG9soN1LgI2wgZnPYAOIytAh/DRJr98WuBsb0HVayet0COeZ0i9NFgK/wgaoJkvPGREH0GE00uQFmO/8Gjaw78VYwP1NSlynw2ikyQuZXkZ2obxdnQrHgh/PlGWwQe8PA6/HggO8If59bI/zzMR0/za1rmh/ujxuxtJx9dS65N3gxh7HnBxv/1iB818T75ueCLWo3vXjY79f4DqAjXB9BJs6Nsu8+GT7YCNu/44l4s6pfdbFRgs/DPwVOI/pHQ/fjI26TY+6/jxwF1a4n2b6VL+3xvtcD3wlo+cu4KQB9kxmzld1auF+6QNmf4S9/DyJRezYNd62K2ZfuiDPja+7QXxsMmXulVg6vbPP+rJpkkeHag+DZLrfIgW6KsqDfPah+pTX76P3fXBqvG1mhfP3Y5j7VIsWLVq0aHG5NEHbNmrRokWLFi1L6c83sWhtw/Lq+Nzv67PPUfE+NzO9Ma1p2k5/LVq0aNGipemlCdq2UYsWLVq0aEmWJmjbRi1atGjRoqXppQnatlGLFi1atGhZynB04mMmhjwuj4n4XGUGrD7IVJ/7hBnAn+k945XL/tjJ8Uvp7nudt86X/tg+pUmapbQ3YLWtNOm1Ps1KwLPYIKIydBi9NAEbx1N2AqwObtIEeqeLqzSZw9RYp8tpZ8Bq2bLiIk0WUmw24iJ0GI00+RTw0xI25NHBT99T9XnyOWzW2Rkl7OrgZ5okxxf1xxfRPb7qzHh9HpvE5zkvta5If7pe3B8fu0pq3epYcIK/AsvmHHNxfMxbC5z/51gapimqd494v38rcB3ApkhfCnwjZ9ve8bYrsJHSGwIXMDXF7vpYYpyEJfLm8b7p6ednAD9j6qXyyPiYDbCBcFvH19gaG5G9CpaAS7DBrmm+SP+pll3QL33mAvdgabIDFqXi3VhkBoCP0D118euxG2kGsFd87p8Du2FpsmqP9WtQb5p0qFZhSgrkv5Q8fhiUB/kk0QuKDBpeOWfdVtgL6qPYvZzlxPj8K5XUJ4QQQgghhBAiHC7EGhyH5QvYB6E1e2w/Fqtb/ozyUU6FEEIIkc9c7IP77Vik6P1aVSOEEEIIIYQQIiRUpxRCCCHK0cGPAauPAdfmrP9WfM4Nc7a57I8N+X2v35mzzqf+2L6kSfZbepsDVttIk15997Ppsna8z/ZD2pTQYbTSZBbwNmyg1oIhbUroUH+aQP90cZUm5wCfif+/nHYGrJYtKy7SZCHwOPBHbNbZs7EJDcvQYTTS5Hbgs8DXsQAPtwCH0vzgTF/GR2WfsctiATGOKWET+JsmMJw/Pgb4A7BxfOxLsZntP9hD+7vi8xyaWtevP93G2NjJLDOZGsOVN7D6rHjbv2bW74rN8PpXbBbhQdxBd3CTQf3/Et4aa/hEgesA5qSWAh/N2XYclklrp9YdiGU2wA+xUeZpXgf8JbPutdho3o/F59sytW0fphcEgHViTTvl6OkXhcEF/dLnYqwg9JqB8n/pHll9HHB1/P/HMSeQHSyYt77uNOlQrcKUTOWbjW7jAuVBPhvHx/97gX2vZ+rF6yTsYb0EmzX5dT2OeX98/rzBrHVzEDa9+hLgjAau54LQbZD+9gndhtD1Q/g2SH84hG6r9LdP6DaErh/Ct0H628dHG64DLitx3GLgqh7bDsTqlUuA/8Aaw7NLp8Q1hRBCCGGsDWwW/78G5pef154cL99xhiV0G0LXD+HbIP3tE7oN0i+EEEKIpvCtTgnhv0uErh/Ct0H62yd0G6S/fUKwoUO1vr/7YP2aI+AH8bl+m1p3SoFzLB8fd27OtqTP9SF9trnojw35fa9974/tS5pkaXPAahtp0m99mnOwiddm9dmnHx1GI00WYIPWl2Djdvboq74/HepPE+ifLi7S5L3ATUzNPHg57QxYLVtWXKTJ7sAbsfKyC3ApcC/lgs13GI00eTJeTgJejg00fIx8nzmIDn76nirP2Ldgz5V1+uzTjw5+pgkM549nYmXkOWwcZN4g0TRnY+k2N7WuX3+6D8fnvRT4cnyt/8beB5di9+lLc45bA5uNfik2I+wpwDfjaz9DdyCQXjyETTyRpp/eNFvE138QG8f2tkEHnBwfsEvOtvOw0eNpTsAe4C+Kj3scu0mT5QlsBtUs12AJsVvO+bKGJS+iO2bWHw/8qqclbuiVPon9W/c59nbgiMy6b2MRYMAKR15lIW993WnSoVqFaTZWqL9X8vhhUB7ksxz2EDy7wL4fwV7C/opNVX038CX6R8nYlXx762Zj7NnwRqwxfEXH13NB6DZIf/uEbkPo+iF8G6Q/HEK3VfrbJ3QbQtcP4dsg/e3jqw1/YPigVFth9cZ/7rF9Yby933L58FKFEEII0YNFwHotXdvXd5xhCN2G0PVD+DZIf/uEboP0CyGEEKJN2qxTQvjvEqHrh/BtkP72Cd0G6W+fUGzoUK3v70L6fz+dLHCOZNKhk3K2LQvch010k8Vlf2zI73vte39sX9IkS5sDVttIk37r07ruJX/24KIdtdFGAAAYiklEQVR0GI00WRZLhy2x58CDwPw+GvrRod40gcHpUnea/CM2W+DGqXWX086A1bJlxeW9k7AiNs6rVx+ffnQYjTR5mu7ZyU+KrzUsHfz0PVXKySXAhQP26UcHP9MEhvPHb8H60b0VG/B9APAw8J6cfedgYyi/m1o3qD/d/FjbLdjzewnwN+BG7D2x36DyVbGBondj5fkh4HymZpsdxHJ0Bz0ZpDfLxzF//Bw2Y3FpfgMclll3ATYrxhuwwW8b5ixzM8fsjM3k+Bw2ojbN+cCpmXXLYomeHeH7ReCKYY1wxN6Yxl6jtFeIt786s34xFsEB4E7ypwXOW+9jmizCRmi3hfLA7tEbHZ17deyefe+gHSuSDKb1iQh72BcldBukv34iVIbaJkJ50DYRYet3hY+2RoSdVxHjox/CtyF0/RC+DdJfPxHhl6HlsfrusBEDP4U1WL24dkXtswrWUL9BzrYQIjv3I3T9EL4N0t8+odsg/e3Thg3fov9Hoi2BO4AZzcjpwsd3nIjw39Mi9K7cJhFh6weVIR+IGJ888FF/m/SqV+pdqH2kv31Ct0H62yd0G1SnzMfHd4mIsN+FIsJ+n44IWz+oDPlAxPjkQej6wU8bfGUn7DvtO3tsPzre/vKC56ujPzbk970OpT92lqbTJEsdA1brxmWa9FsPNiDmfvJnnGuTNtMkzY+BrxbYryn6pYuLNOlg98yS1LIUGwexBBug1TZNp0kvfoI/9eg20uT3dN8rB2Bj3Xyg6fFRaV4EPBtr8Ik2/PFi4EOZdZ/AZkDNchjdQTh87k+3IabtuNS6VvSuhD2kd8isX4y94O6OZezzB5xnM2xg64HYqOGLM9t/hzXEZbkem942zZ3kR4Npgz2w9Fm5x/b5WKatnVqXVBC2xiIUPAu8MnNcr/XgX5qcgWldoaXrKw+sQv93ej+Aq3I78AVH5wZLu3RkqPMcXmsYIoo3lIRug/S7IUJlqG0ilAdtExG2fhf4amtE2HkVMR76IXwbQtcP4dsg/W6ICL8MvQLTs+mQx92BRXsLiaOxwE+PYBE2LyQ/0ujJwJk560OJ7NyL0PVD+DZIf/uEboP0t09bNmyKRW6dk7NtNawtd7uGtGTx9R0nIvz3tAi9K7dJRNj6QWXIByLGIw981e+CKvVKvQu1j/S3T+g2SH/7hG6D6pT5+PouERH2u1BE2O/TEWHrB5UhH4gYjzwIXT/4a4OvHIKlU17wXbCAxb+n+CxpVftjQ37f65D6Y2dpI03S+Dhg1VWa9FsPNjGaj4NVob00yXIZcFYxyY3QL11cpMkL4vOmlxuBr8f/txmYJ6HpNMljeWxWxOMG7dgQbaTJ14GrMuv+hXIzrLqgjfFRCQux8jF7WNGOacMfPwQcmll3NN0DVlcA/oQFCUvjc3+6N2Jps3NqXSt6t8MyYKXUutUwcQuwiKAPYDOkvhx76d0Vi26SPNRfBPwRm/IVrPHrOaaPHp7EXiTWwZxFwv7YFLUHAZsAnwMei8/pA6thDXlfB16GNSi+l6mPUWthtu4b/94G+BVTAzxfGf+fbXjstR78S5P9sfKwfUvXVx7Au7A8cPUSfgpuHz5rAL/GnhFr0duRNE1E8YaS0G2QfjdEqAy1TYTyoG0iwtbvAl9tjQg7ryLGQz+Eb0Po+iF8G6TfDRHhl6H3Y42OroIh+cQlWF16Pta+9R3gPmDV1D7Pw4Kv7dh1tJ+RnSMUXbtNIsLWDypDPhAxPnkQun5o14ab6J4RfTngSiz6b1v4+o4TEf57WoTeldskImz9oDLkAxHjkQe+6ndBlXql3oXqJyJs/aD3aR+IGJ88CF0/hG+D6pT5+PouERH2u1BE2O/TEWHrB5UhH4gYjzwIXT/4a4OvfBnrn9+PnYDjKRYgo2p/bMjvex1Sf+wsbaTJSsDm8bIUODL+f706DKoBV2nSb/3pWNCwnePzJ8tK+EEbafJprM1pHtY2dVJ8jd1rsKcu+qWLizTJ43LgtLIGOKCNNDkFeBU2a+I2wEXY/RTCc9ZVmmwFPIO9b2yIzfT9N7rrym3RxvgosP5iv8eeL77Rhj+OgHuAPbFn7b7Y2MnPZvbbBHvXnTe0Ve3xSeAp7BuGM4p0QNwMuAt7+Ut4OSbuDuAvmGN7PjY19C1YlNDF2EvSqsAPsAfbifHxi7DRw+noJx8H9sMy9MTU+nOAD2NT596COdY9sBvBBx4CXo89sK8HrsUeWH+Ot98HHAP8F1YpOAQ4G4v+8gRT6ZudPrrXevAvTb4HPM700dVNojyACzCn6SoPvo4NNF/T0fkfAdYHforl1yMOrrEH1pBxF8Wnux8G1zZcgD1vs5EX6sK1/rnYS//twC+w532dNFGGrsG034abqDJN2DATixbkohw1oX8Se4e4BXvnqBvXNszDIljdDvwSWL3m87vWvwBL+2R5AtinxvM3UYYOwe7h27D3glkOrlEE17aG7vPArd8L3eeBextGwedB2H5vkrB9Hrj1e6H7PHBvQ1mf9zrgXKwBb9R5HfB/sTS6FeuI9UKmB+PaE0uLqzPH3gn8G7AFYUZ2Dl0/hG+D9LdP6DZIf/u0bcMFwNtTv2dgH8wuA/6nYS1p1M47GLXzDkZ13mK4qvOqnXcw81A7bz/GqZ23bL2y7feIOgjdBulvn9BtkP72Cd2GtvX7WqcEfUctQsj1StUpBxN6nRL0HbUI89B31H74+h3VV3YEvjtgnyuBE8jvC52lan9syO97HVJ/7CxtpMmWwM3xAjYu42ZscIkPuEqTfus/gI1VuRQLwp0sR9ZkU1XaSJO1sNlUf42ly1bYmJ6La7KpDvqli4s0CYE20mTd+Dy/Br6NjQPbljCes67S5EbsHeYt2DvBicCx2OB4H2hjfBTALlhwhP+uyY46acMfH4bVmU7Hxk5+FvgKUxN5JtyBDVidLGdaK0wAP8TGAQrhPWdjFUTRHhfGiytuwF3UiC2xRvgXDNqxINlGhOUwRzIXiwJwO9On+044BhucnyzPYJGb0uvyZtaBem3I6p8FvBpzsoMaCsva4DoP/gFz6GBR0BaTH5HBV/2zmIrWNgt7ydk85zifyxDYPfwN+pcjn/NgksERsXzOgyuwiHUAc7BnUx4+50HCSsCD5EeU8VX/msBvgeWxj64XAnvVdK1hcVnWivo88Devivo9X/W79nng/nlVxOeBv3lQ1O/5nAeTFIsC6XMeFPF7PutPcOHzwG0ZKuvznod98N2iBk0hsjaWJ+mOxZ8HfpSzr6+RnSMUXbtNIsLWDypDPhAxHnkQun5o34bdsPeOJBrsDthgmHRnsQUNawK182ZRO6+feTAK7bzgrs6rdt7pqJ3XP/0+tfNmKVqvbPs9ohcR4bwL5RERtn7Q+7QPRIxHHoSuH8K3oW39vtYpQd9RszT9HRXc5sG41CnB7zIE+o7adh7oO2p431F95UVYWu3QthAhhBBCCBEMc7D383e0LUSIorwSa7z0ZWrwceS1WKV+FUfn3xcbtOqCg+iOKDAXOB+LCvRLbBp64t8vjP9/BVORJi8EPoNFTDg48/v4+FwJJ2JTgGdZFZvWPVnOA07NrFsh57hhbCij/+B4nwkGd2Qqa0NV/UVtSFiEReBoWj99bCiqfwXgJvIbnX0uQ2tgUZR2pn858jkPJhnc6OxrHnwI+PEA7VVtaPI+/ifs40Ub+uljQz/9R2ERoV4ALIPlx1Y9bHCNy7JW1OeB/2VtgnCfV2nq9nnD2ODS51WxoYk8KOL3fM6DSYp9aPX1Pi7q93wuQwkufN4wNjTp896J3TfjyjlYeqc/XJ8PnJmz7/JYw92EIy1lo/xHFOsQ51o/lIvyH1G8U6JLG8pG+Y/wQ39CmUj/Ef6UIRg+yn+EP3kwSbko/xH+5ME8ho/yH+GH/rJR/iP8KUODovxvinVE2sDR9cuidl6181bVn2eD2nm79bus86qdV+28Temnjw2htPNmKVqvbOJdrky9MsKfd6FxrFOCP+/T4LZOCc3YUGbmuAh/8mCS4euVEf7on4e7OiW4taHKzHER/uRBv3qlr3VK0HfUQTa4/o6aZ4O+o+o7qm95MIm+o7ZdhhLG6Tuqr3wYm9FrRttChBBCCCFEMByAzVibF0BKCG/5CcU/Vgg3/Ax3s6DOAK4GXuXg3KcB3039Xha4FYsKCFORxGYzvTGiA5wQ/3838H9S29K/3xRfI+EDFPtIEBXcD4rbMJnap0Mx/QkTDPdBDYrbUFU/FLMBLApa0YaSiHr1ly1DCdcDj2KNWEWI8KcMnQlsx/DlKMKfPLgba/C/EWv0LEKEH3mwD/AdphpPP1lQE/h5H5+PBTIoQoQ/ZegwbOa4h5juF5rGZVkr6/PAv7I2QbjPqwQXPg+aee8Y1ueBX3lQxu9F+JMHZXwe+HMfl/V7Ef6UoQQXPg/cl6EyPu86pqI5jxsnA/diH8fTXAL8Z87+ivI/HUX590t/8tvl7HGK8p+vN/17EkX5b1t/wqhG+X9JfP35NVy/TtTOq3beXkSUywO18+brd1nnVTuv2nl7EeFPGfKlnTfNMPVK1+9CReuVvr4LjUudEvx+nx6FGckH1SnB7zyYxN2M5KHXKaG5We371SnB3zwYVK/0tU4J+o4K7X5HhWbqZKDvqL1o4n1a31H74/o+1nfUwbguQz7WKcuyCHh32yKEEEIIIURQXMZw37dKM7OJi4ix4SgsQmuvyEfCPR8FjqC7MboOlmKNtS4GJW+ORaZMSGZzTSKF/g14CtgI+E1qv02xaJDPxxom/j1en/2d17i5tA7hKYrYsAHl9DdBFf1Q3IbVgK8B76HePHBdhhK2wT40b079H29clqGdsPS+pmbNaZrIg+2xiHhvwD7+LajVArd5MBtr7D8c2Bqzo2ijbRP68zT3ug/mYDOrX1yX8BjXZWgVYE8sovO62D08Ua8JhXFZ1nz3edC+32vKZ7jyedDMe0fbPq9KHvji96rkgc8+L09z036vqfvYlc8Dt2WojM/bLb7WlcMaMgJ8Fptd9jVMT2+wzmir5ByzORZ9+a+pdWUjOx+U+f0xrDPxYuBxrNNCdpBSwpdiLclyQc66nznQP8iGg7Cy/GgP3VX159lQJTp4Vv9e2CwyAH/GZvVZ3bH+XjaU0X9QvM8j8d9l4yUP12WICjasgc1E9OUe2qvoH8aGKnlQBF/v40OxTraJX0h8a9P66WND0TzYG4uS//ca9Q9jQxn9Hex9avn47wrA/Znrrxr/faCHvrZQO6/aeesgbYPaebv1u67zqp1X7bxVGad23oRh65Wu3+WK1it9rZO1UafsZYPLOmUVG5p4n3ZZpxzGBpd1yio2NNG2UgRfy5DrOuUwNrisU1axwXUedOhfr/S1Tglu3+V8r1OCf/VKfUf1s06p76j+1inzNOs76vC4LEM+1inLshMWVOOstoUIIYQQQohgeCkWcPHUJi6mAauiTm7AGl8/2LaQMeYyLHrYOx2d/1Ysjw+o8ZwzsAaDdCPDAqw8ZZnPVOMOWLS92+L11zDViJn9/UescT9hXeBPVYWnKGpDWf2uqao/2TbIhuWwj+InUW/DZxNlKM0j2L22WyXV03FdhrYDdsUiy30D2B34r9rUN5cHyX17L/B9rNG2LlznwT1YhMffYx9oL8IaWOuiqfsYLNrjD4AnK6ueookytAsWSfhh4Ange1jjf9O4Lmu++7xkW1t+r6nnlSufB82+d7Tp8xJtIfu9UfV5yba2/F6T734ufB64L0PD+rzlgA9hH8bHjVOBd2ARmW/P2X4z1lCXJfuhfFnsXv0c8HLsOfQn7KN30ulrAdPvoT9jnVr+M/P7l1in4oR7sI4veTyMfYxPlkdz1j3hQP8gG/Jmpa1Tf9aGXvpn97GhqP4tgWWYnid16+9nQ1X918frf5y5XlUbiuqvkgcnA8cCz+Vcv6r+YWwoq38pNpvMoCj/vt7Hv4+1DIry73MZSngLcE7N+oexoYz+zwCnAH8A7gN+hZWlNAvi82cHsraJ2nnVzlsHWRvUztut32WdV+28auetyji18yaUqVe6fpcrWq8MoU7WD1/f5YrWKavY0NT7tKs65TA2uKxTVrGhibaVIvVKX8uQ6zrlMDa4rFNWscF1HgyqV/pYpwR9R022+VSv1HdUP+uUiTZ9R/WvTpls03fU8vj2HdVnjsVmi326bSFCCCGEECIYjsHeIx9s4mKzm7iIGCuOxBpj/gdrCBXNcygWueo8piKf1snnsYhTVwN313C+pcDKmXX3Y40FYAPrVwEewqI8Jo3522INTncB7wIWpY6fn/l9A/AyrOH5QSzq1i4FtHVqtqGs/ip0CuxTVT90a87+ngFETD0fitIpsE8TZWgO9pHoASwK6WuB/6hJ/zA2lNX/6XgBiwp3KBapsgidAvs0kQcrxud5FFgJ2Bk4tyb9w9hQVv+NWJTQ1bBGz1dhH2GK0GlAf57mXs+i/YHTCmhK6BTYp4kytBiL7rg81ug/weDI2i5wXdbK+jzwr6wNS6fAPk2UNZc+D9yXobI+D/zJg7J+r4h+cJ8HZX0e+HMfl/V7deqvUoYSXPg8cF+GhvV5C4GPYnk1TpyOdSreB7N9rXj9Y/ECcAnWuWt1pjfUbY511kzIi+z8UrojO1+LP1H+y+oHf6L8Jzbk6YfyeZDgOsq/qzKUZhvsefNNujtvVKGIfihvQzrK/0RNmrO4zoPtsY5Ta8fXWYQFgvNBf57mXlH+N8fs+G58jaL12Dr008eGYaP871+T7jQuy1A6yv8TWFvvBHB56lw7YR3CfELtvGrn7Uen4H5ZG9TO212GXNZ51c6rdt5+dArsM07tvFC+Xun6Xa7peqWrdzmXNFEna7tOCdVtcFWnBLd54HudEtqvV7ouQ67rlEVtCL1OCe7qlT7WKUHfUfM0N/kdFdzXyfQdtR79+o6q76hV9es7qh91yirsjQVi+H7bQoQQQgghRDBsi7U1fqmpC2qGVVE3jwAHUbwxRtTPQ9jMOCc7vMaRWKOJKyJgA+yj103AJvH6i4G9gLOBPbBGhiVYFK10w0L29xLgCOBS7CPOGdQbJTGPiG4byuoH+1jxzfi4e7CGE5dEFNefpzn7e3usoWofLALaLfE+roiotwytgn2wWQT8DItme5Ez9UZEvWWoaSLqzYM1gauAXwDXYR/5s7OL1E1EfXnwLHAU9pFxEfYBrGijeVki6r2PwRp8twB+5Ep0ioh6y9B18bE3x+t/i80a7gMR9ZU1331enua2/V5EvWWtaZ8H9ZYhX3wehO/3RtXn5Wlu2+9F1F+GmvR5UG8ZGsbnrYp1CqlzEFUofADryHUpFoU7WY5M7XMrNpvIW1PrFOVfUf6r0vTsfY9Qb6R/RflXlP+qKMp/9Sj/K2C2faWiHU0QoXZetfNWI0LtvD7Wd2F067y+1XdB7bw+t/OWqVc28S7XdL3S1bucK5qoD/hQp0y0hTx73KjWKfM0N1mvbKIM+T4jeSh1ykRb3fXKkOqUoO+obdcrI/Qd1cc6JYRfrxzVOmWe5rbrlRH6juprnbIsa2IDc9/fthAhhBBCCBEMKwLHYd8qnmtZixCV2RV4S9sixpy9gFe3LUIIIYQQQgghhAiY1wF3ArP67HMY8IX4/5lYgKf3Yx1TwSLU/R2bXeK9wLGpY7O/Z2OzJszFOo/dDqxTyYLBDKMfBtuQMAF8q2ateeTph/J5MAPr4LDQjdxc6ixDYJH+Xxj/vzzWCW+v2lVPUXcepJmgvXJUVv+K2MAFsCj/NwFbOVE9Rd338Sysw85q2D3xDdzMKJPgqgx9H+us1AR1lqFtgZ9j9+8srGPx3qnthwI/rFe+EEIIIRwyqF5Z97tc0/VKV+9yEzRTF4B63+V8qVNCeRuarlNC/fXihAn8q1OCf/XKustQ03VKcFOG2q5Tgpt6peqUQgghRNiciA0aFkIIIYQQoiifBDZs+qKaYVW44kfAhW2LGHMuwiKKCSGEEEIIIYQQohyXAKdhM9L0IkJR/hXlvxoR9c681nSk/whF+VeU/2pEKMr/MFH+n8Y6MwshhBAiDAbVKyPqfZdrul4ZUf+7nM8zkodSp4TwZ48b1Tol+FevjKi3DPk+I3kodUpwU69UnVIIIYQIm+OxWdSFEEIIIYQoyqexQJdCCCGEEEIIIYQQQgghhBBCCCGEEEIIIYQQQgghhBBCCCGEEEIIIYQQQgghhBBCCCGEEEIIIYQQQgghhBBCCCGEEEIIIYQQQgghhBBCCCGEEEIIIYQQQgghhBBCCCGEEEIIIYQQQgghhBBCCCGEEEIIIYQQQgghhBBCCCGEEEIIIYQQQgghhBBCCCGEEEIIIYQQQgghhBBCCCGEEEIIIYQQQgghhBBCCCGEEEIIIYQQQgghhBBCCJ/5f30SHG4uoPELAAAAAElFTkSuQmCC\n", + "text/latex": [ + "$$f_{next 0} = f_{curr 0} + \\frac{- f_{curr 0} + \\frac{\\left(\\frac{9 \\left(- \\frac{- f_{curr 0} - f_{curr 1} - f_{curr 2} + f_{curr 6} + f_{curr 7} + f_{curr 8}}{f_{curr 0} + f_{curr 1} + f_{curr 2} + f_{curr 3} + f_{curr 4} + f_{curr 5} + f_{curr 6} + f_{curr 7} + f_{curr 8}} - \\frac{- f_{curr 0} + f_{curr 2} - f_{curr 3} + f_{curr 5} - f_{curr 6} + f_{curr 8}}{f_{curr 0} + f_{curr 1} + f_{curr 2} + f_{curr 3} + f_{curr 4} + f_{curr 5} + f_{curr 6} + f_{curr 7} + f_{curr 8}}\\right)^{2}}{2} - \\frac{3 \\left(- f_{curr 0} - f_{curr 1} - f_{curr 2} + f_{curr 6} + f_{curr 7} + f_{curr 8}\\right)^{2}}{2 \\left(f_{curr 0} + f_{curr 1} + f_{curr 2} + f_{curr 3} + f_{curr 4} + f_{curr 5} + f_{curr 6} + f_{curr 7} + f_{curr 8}\\right)^{2}} - \\frac{3 \\left(- f_{curr 0} - f_{curr 1} - f_{curr 2} + f_{curr 6} + f_{curr 7} + f_{curr 8}\\right)}{f_{curr 0} + f_{curr 1} + f_{curr 2} + f_{curr 3} + f_{curr 4} + f_{curr 5} + f_{curr 6} + f_{curr 7} + f_{curr 8}} - \\frac{3 \\left(- f_{curr 0} + f_{curr 2} - f_{curr 3} + f_{curr 5} - f_{curr 6} + f_{curr 8}\\right)^{2}}{2 \\left(f_{curr 0} + f_{curr 1} + f_{curr 2} + f_{curr 3} + f_{curr 4} + f_{curr 5} + f_{curr 6} + f_{curr 7} + f_{curr 8}\\right)^{2}} - \\frac{3 \\left(- f_{curr 0} + f_{curr 2} - f_{curr 3} + f_{curr 5} - f_{curr 6} + f_{curr 8}\\right)}{f_{curr 0} + f_{curr 1} + f_{curr 2} + f_{curr 3} + f_{curr 4} + f_{curr 5} + f_{curr 6} + f_{curr 7} + f_{curr 8}} + 1\\right) \\left(f_{curr 0} + f_{curr 1} + f_{curr 2} + f_{curr 3} + f_{curr 4} + f_{curr 5} + f_{curr 6} + f_{curr 7} + f_{curr 8}\\right)}{36}}{\\tau}$$" + ], + "text/plain": [ + " ⎛ \n", + " ⎜ ⎛ -f_curr_0 - f_curr_1 - \n", + " ⎜9⋅⎜- ───────────────