diff options
| -rw-r--r-- | notebook/seminar_talk.ipynb | 1360 |
1 files changed, 1360 insertions, 0 deletions
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": [ + "<IPython.core.display.Math object>" + ] + }, + "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": [ + "<IPython.core.display.Math object>" + ] + }, + "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": "iVBORw0KGgoAAAANSUhEUgAADqwAAAAqCAYAAAB8rpsbAAAABHNCSVQICAgIfAhkiAAAIABJREFUeJztnXm4JEWVt99e2ARpWWQbGltABrQbENlk84qgbM |