path: root/notebook/seminar_talk.ipynb

{
 "cells": [
  {
   "cell_type": "raw",
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "<script>\n",
    "setTimeout(function() {\n",
    "\tReveal.configure({\n",
    "\t\twidth: \"90%\",\n",
    "\t\theight: \"100%\",\n",
    "\t\tmargin: 0,\n",
    "\t\tminScale: 1,\n",
    "\t\tmaxScale: 1,\n",
    "        fragmentInURL: true,\n",
    "\t});\n",
    "}, 1000);\n",
    "</script>"
   ]
  },
  {
   "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} \\\\ (\\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": [
    "## Symbolic computation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Experimental symbolic differentiation\n",
    "\n",
    "```Prolog\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": [
    "```Prolog\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": [
    "## Symbolic collision kernel\n",
    "\n",
    "- Discrete equilibrium distribution\n",
    "  - Characteristic constants\n",
    "  - Moments\n",
    "- Discrete BGK collision operator"
   ]
  },
  {
   "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