From 0dc1f594bd085eeb4d43dc816d7348d3bfdf2616 Mon Sep 17 00:00:00 2001 From: Adrian Kummerlaender Date: Wed, 10 Jul 2019 22:20:46 +0200 Subject: Update slides for talk --- notebook/seminar_talk.ipynb | 454 +++++++++++++++++++++----------------------- 1 file changed, 213 insertions(+), 241 deletions(-) (limited to 'notebook') diff --git a/notebook/seminar_talk.ipynb b/notebook/seminar_talk.ipynb index 3f9276d..ecca3ba 100644 --- a/notebook/seminar_talk.ipynb +++ b/notebook/seminar_talk.ipynb @@ -1,5 +1,27 @@ { "cells": [ + { + "cell_type": "raw", + "metadata": { + "slideshow": { + "slide_type": "-" + } + }, + "source": [ + "" + ] + }, { "cell_type": "code", "execution_count": 1, @@ -105,7 +127,7 @@ "## 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", + "$$\\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." @@ -119,9 +141,20 @@ } }, "source": [ - "## Experimental symbolic differentiation\n", + "## Symbolic computation" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "subslide" + } + }, + "source": [ + "### Experimental symbolic differentiation\n", "\n", - "```\n", + "```Prolog\n", "% derivative of x is 1\n", "diff(X, X, 1) :- !. % cut to prevent further rule applications\n", "\n", @@ -153,7 +186,7 @@ } }, "source": [ - "```\n", + "```Prolog\n", "?- diff(x*x,x,D).\n", "D = x*1+x*1.\n", "\n", @@ -196,6 +229,22 @@ "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": { @@ -282,11 +331,11 @@ } }, "source": [ - "Weights of third order Hermite polynomial:\n", + "Weights of third order Gauss-Hermite quadrature:\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", + "Weight formulation [Eq. (39) resp. (43), HeLuo97]:\n", "\n", "$$\\omega_i = \\frac{1}{\\sqrt{\\pi^{d}}} \\prod_{j=0}^{d-1} \\eta_{(\\xi_i)_j}$$" ] @@ -338,34 +387,6 @@ "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": { @@ -385,12 +406,12 @@ } }, "source": [ - "$$\\sum_{i=1}^{q-1} \\omega_i (\\xi_i)_a (\\xi_i)_b = c_s^2 \\delta_{a,b}$$\n" + "$$\\begin{align*}\\sum_{i=1}^{q-1} \\omega_i (\\xi_i)_a (\\xi_i)_b &= c_s^2 \\delta_{a,b} &&\\text{[Eq. (3.60), Krüger2017]}\\end{align*}$$" ] }, { "cell_type": "code", - "execution_count": 12, + "execution_count": 11, "metadata": { "slideshow": { "slide_type": "fragment" @@ -408,7 +429,7 @@ }, { "cell_type": "code", - "execution_count": 13, + "execution_count": 12, "metadata": { "slideshow": { "slide_type": "fragment" @@ -427,7 +448,7 @@ "3 " ] }, - "execution_count": 13, + "execution_count": 12, "metadata": {}, "output_type": "execute_result" } @@ -450,7 +471,7 @@ }, { "cell_type": "code", - "execution_count": 14, + "execution_count": 13, "metadata": { "slideshow": { "slide_type": "fragment" @@ -464,7 +485,7 @@ }, { "cell_type": "code", - "execution_count": 15, + "execution_count": 14, "metadata": { "slideshow": { "slide_type": "fragment" @@ -478,7 +499,7 @@ " f_next_6, f_next_7, f_next_8], dtype=object)" ] }, - "execution_count": 15, + "execution_count": 14, "metadata": {}, "output_type": "execute_result" } @@ -490,7 +511,7 @@ }, { "cell_type": "code", - "execution_count": 16, + "execution_count": 15, "metadata": { "slideshow": { "slide_type": "fragment" @@ -504,7 +525,7 @@ " f_curr_6, f_curr_7, f_curr_8], dtype=object)" ] }, - "execution_count": 16, + "execution_count": 15, "metadata": {}, "output_type": "execute_result" } @@ -516,7 +537,7 @@ }, { "cell_type": "code", - "execution_count": 17, + "execution_count": 16, "metadata": { "slideshow": { "slide_type": "subslide" @@ -534,7 +555,7 @@ "6 + f_curr_7 + f_curr_8]" ] }, - "execution_count": 17, + "execution_count": 16, "metadata": {}, "output_type": "execute_result" } @@ -544,9 +565,20 @@ "moments" ] }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "fragment" + } + }, + "source": [ + "$$\\begin{align*}u(x,t) &= \\frac{1}{\\rho} \\sum_{i=0}^{q-1} \\xi_i f_i(x,t) &&\\text{Velocity}\\end{align*}$$" + ] + }, { "cell_type": "code", - "execution_count": 18, + "execution_count": 17, "metadata": { "slideshow": { "slide_type": "fragment" @@ -562,7 +594,7 @@ "" ] }, - "execution_count": 18, + "execution_count": 17, "metadata": {}, "output_type": "execute_result" } @@ -583,12 +615,19 @@ } }, "source": [ - "## Discrete equilibrium" + "## Discrete equilibrium distribution" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "$$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)$$" ] }, { "cell_type": "code", - "execution_count": 19, + "execution_count": 18, "metadata": { "scrolled": true, "slideshow": { @@ -639,7 +678,7 @@ " 36 ⎦" ] }, - "execution_count": 19, + "execution_count": 18, "metadata": {}, "output_type": "execute_result" } @@ -665,12 +704,23 @@ } }, "source": [ - "## BGK Collision" + "## BGK collision operator" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "fragment" + } + }, + "source": [ + "$$\\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))$$" ] }, { "cell_type": "code", - "execution_count": 20, + "execution_count": 19, "metadata": { "slideshow": { "slide_type": "fragment" @@ -683,7 +733,7 @@ }, { "cell_type": "code", - "execution_count": 21, + "execution_count": 20, "metadata": { "slideshow": { "slide_type": "fragment" @@ -699,7 +749,7 @@ "" ] }, - "execution_count": 21, + "execution_count": 20, "metadata": {}, "output_type": "execute_result" } @@ -713,7 +763,7 @@ }, { "cell_type": "code", - "execution_count": 22, + "execution_count": 21, "metadata": { "slideshow": { "slide_type": "subslide" @@ -727,7 +777,7 @@ }, { "cell_type": "code", - "execution_count": 23, + "execution_count": 22, "metadata": { "slideshow": { "slide_type": "fragment" @@ -735,141 +785,16 @@ }, "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⋅⎜- ───────────────────────────────────────\n", - " ⎜ ⎝ f_curr_0 + f_curr_1 + f_curr_2 + f_curr\n", - " ⎜────────────────────────────────────────────\n", - " ⎜ \n", - " ⎝ \n", - " -f_curr_0 + ─────────────────────────────────────────────\n", - " \n", - "fₙₑₓₜ ₀ = f_curr_0 + ─────────────────────────────────────────────────────────\n", - " \n", - "\n", - " \n", - "f_curr_2 + f_curr_6 + f_curr_7 + f_curr_8 -f\n", - "───────────────────────────────────────────────────────── - ──────────────────\n", - "_3 + f_curr_4 + f_curr_5 + f_curr_6 + f_curr_7 + f_curr_8 f_curr_0 + f_curr_\n", - "──────────────────────────────────────────────────────────────────────────────\n", - " 2 \n", - " \n", - "──────────────────────────────────────────────────────────────────────────────\n", - " \n", - "──────────────────────────────────────────────────────────────────────────────\n", - " \n", - "\n", - " \n", - "_curr_0 + f_curr_2 - f_curr_3 + f_curr_5 - f_curr_6 + f_curr_8 \n", - "──────────────────────────────────────────────────────────────────────────────\n", - "1 + f_curr_2 + f_curr_3 + f_curr_4 + f_curr_5 + f_curr_6 + f_curr_7 + f_curr_8\n", - "──────────────────────────────────────────────────────────────────────────────\n", - " \n", - " \n", - "──────────────────────────────────────────────────────────────────────────────\n", - " \n", - "──────────────────────────────────────────────────────────────────────────────\n", - " \n", - "\n", - " 2 \n", - "⎞ \n", - "⎟ \n", - "⎠ 3⋅(-f_curr_0 - f_curr_1 - f_curr_2 + f_curr_6 + f_curr_7 \n", - "── - ─────────────────────────────────────────────────────────────────────────\n", - " \n", - " 2⋅(f_curr_0 + f_curr_1 + f_curr_2 + f_curr_3 + f_curr_4 + f_curr_5 + f_cu\n", - "──────────────────────────────────────────────────────────────────────────────\n", - " \n", - "──────────────────────────────────────────────────────────────────────────────\n", - " \n", - "\n", - " \n", - " \n", - " 2 \n", - "+ f_curr_8) 3⋅(-f_curr_0 - f_curr_1 - f_curr_\n", - "──────────────────────────── - ───────────────────────────────────────────────\n", - " 2 f_curr_0 + f_curr_1 + f_curr_2 + f_curr_3 + f_c\n", - "rr_6 + f_curr_7 + f_curr_8) \n", - "──────────────────────────────────────────────────────────────────────────────\n", - " 3\n", - "──────────────────────────────────────────────────────────────────────────────\n", - " τ \n", - "\n", - " \n", - " \n", - " \n", - "2 + f_curr_6 + f_curr_7 + f_curr_8) 3⋅(-f_curr\n", - "───────────────────────────────────────────────── - ──────────────────────────\n", - "urr_4 + f_curr_5 + f_curr_6 + f_curr_7 + f_curr_8 \n", - " 2⋅(f_curr_0 + f_curr_1 + f\n", - "──────────────────────────────────────────────────────────────────────────────\n", - "6 \n", - "──────────────────────────────────────────────────────────────────────────────\n", - " \n", - "\n", - " \n", - " \n", - " 2 \n", - "_0 + f_curr_2 - f_curr_3 + f_curr_5 - f_curr_6 + f_curr_8) \n", - "─────────────────────────────────────────────────────────────────────────── - \n", - " 2 \n", - "_curr_2 + f_curr_3 + f_curr_4 + f_curr_5 + f_curr_6 + f_curr_7 + f_curr_8) \n", - "──────────────────────────────────────────────────────────────────────────────\n", - " \n", - "──────────────────────────────────────────────────────────────────────────────\n", - " \n", - "\n", - " \n", - " \n", - " \n", - " 3⋅(-f_curr_0 + f_curr_2 - f_curr_3 + f_curr_5 - f_curr_6 + f_cur\n", - "──────────────────────────────────────────────────────────────────────────────\n", - "f_curr_0 + f_curr_1 + f_curr_2 + f_curr_3 + f_curr_4 + f_curr_5 + f_curr_6 + f\n", - " \n", - "──────────────────────────────────────────────────────────────────────────────\n", - " \n", - "──────────────────────────────────────────────────────────────────────────────\n", - " \n", - "\n", - " ⎞ \n", - " ⎟ \n", - " ⎟ \n", - "r_8) ⎟ \n", - "────────────────── + 1⎟⋅(f_curr_0 + f_curr_1 + f_curr_2 + f_curr_3 + f_curr_4 \n", - "_curr_7 + f_curr_8 ⎟ \n", - " ⎠ \n", - "──────────────────────────────────────────────────────────────────────────────\n", - " \n", - "──────────────────────────────────────────────────────────────────────────────\n", - " \n", - "\n", - " \n", - " \n", - " \n", - " \n", - "+ f_curr_5 + f_curr_6 + f_curr_7 + f_curr_8)\n", - " \n", - " \n", - "────────────────────────────────────────────\n", - " \n", - "────────────────────────────────────────────\n", - " " - ] - }, - "execution_count": 23, - "metadata": {}, - "output_type": "execute_result" + "name": "stdout", + "output_type": "stream", + "text": [ + "Eq(f_next_0, f_curr_0 + (-f_curr_0 + (9*(-(-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) - (-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))**2/2 - 3*(-f_curr_0 - f_curr_1 - f_curr_2 + f_curr_6 + f_curr_7 + f_curr_8)**2/(2*(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)**2) - 3*(-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) - 3*(-f_curr_0 + f_curr_2 - f_curr_3 + f_curr_5 - f_curr_6 + f_curr_8)**2/(2*(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)**2) - 3*(-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) + 1)*(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)/36)/tau)\n" + ] } ], "source": [ "lbm = bgk_collide(moments)\n", - "lbm[0]" + "print(lbm[0])" ] }, { @@ -883,6 +808,23 @@ "## Common subexpression elimination (CSE)" ] }, + { + "cell_type": "code", + "execution_count": 23, + "metadata": { + "slideshow": { + "slide_type": "fragment" + } + }, + "outputs": [], + "source": [ + "a, b, n, m = symbols('a b n m')\n", + "exprs = [\n", + " a*b + n,\n", + " a*b + m\n", + "]" + ] + }, { "cell_type": "code", "execution_count": 24, @@ -891,6 +833,34 @@ "slide_type": "fragment" } }, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAUIAAAAXCAYAAAB+mUZsAAAABHNCSVQICAgIfAhkiAAAB55JREFUeJztnHmoVFUcxz+a+jKfaFm0p622+uwZlZFwI1soIigijIqpjMoWKKEoCiwCw4osiMgSXoWm0WvRVoO6baapqWmbUj6lMpc2W7TS7I/vmebOnXOXuc9ZPR8Y5s35nW1+9/c72+/MA4fD4XBYeRpYD/Qznz1ge+D1ZRllHWIc0t2VKfKOMHmvipDvSfHz2L4jOthNPNLbiMNRKeJ8I8mvijgB+Be4JZDmmQp8YAJwQxllHeJJpMO2lPlfBNYCrRbZbug5TAC6qK+B0CfeRhyOSpLkG3F+VcQc4BegbyDNMxVOyFDWIRYDm4FeKfOfiHR+R0I+n/oaCCfUthsOx//4lPqG1a96hjIdAYwGnkNOWw7dKdvstADHAJ8CW1OW+RhtL68BdqlQv5qFHDJur7bd2CnJ0Vi6t/pVeCC8EugBzMzQQFLZOUhhF4TSewAdRnZfhnZrxSXANGAF8BvwM7AAuMKS9zigN7AIGA50Aj8CfyC9HBvRxgzgIDTBNANnoOc8EX3n6cA64HdgLnBSDfrUTHZ5GurvA0A78BLwE/Ar2hLuY/IdjXS/3sheQXZWbWql+0S/WohWLOFAh0fytieqbJ42YBvwOcUrnAdN3VNi6q43+gP/AB8BU5FjTwU2ou9yWyj/NSb9VTT4PQ9MAt436euBgZZ2RlMw7Ch8GmdrfKvJ04l2DbOQHl416euQbsslR/ZVSTPZ5XjU51nAn0jP9wPLTPprwHlo4pllZAsonO1mIUf96t7H7huxftUPDWTLLDKPeCOPKxukw9STM5/vMJ9nUro6rWdaKcyuQfZFq8NwxHQK+p4bgeNDsqeJPgscYGQfx/TFp3EGwhkUBv72kKzTyEZlaDtH97ZnHTSHXU5D/f4BGBZI7492LFtRoGBkQNYHWIOCnLtmaDNH/erex+4bsX51hBHOscg84o08rmyQA9BKoAtFFbcDb6CH0Sx8hRw9yEKir86cZGTPRtS3GRl2FD6NMxB+ZfKca5Hda2RnZGg7R/ecsVns8gvU97MtssVGdqlF9paRDcrQZo761b1PtG9E+tVICiNxGI94I48rG2Yihfs9H6Jwt41xwCpgCzpby7JSqBS7A3cB89BMu43ie0ufBPL2Bv5CD9oWMR5syrwY0dZ3xAdYfBpjIGxFelqNzoDCTDflD0lop4vSe2Jxr46U/W90u+yH9PtNhPwndC7d2yJbAWxK0UYXjaV7n2jfKPKroGPmI71ZlsfllN0Q+PsqdJYR5mLgYfTFPwCuA15Hh7xrMvRvRzIMrXz3RkvrGcjI/gEOBi4HlgbyH4dmt9nYB7TB5n11RHt9aY4o/HC01ckfkIdpRwf3qxLqmUzpeepw4HzgKeSsQZak7F+j22Ub0u9bFtkQNHm/gOw0SCtwKBqAkmgm3Uf61X7IQD+wyDziZ/u4skHGoLOItSb/YxH55gNPhNJWopmj1ixFM5Jnkd2Dvtf1gbSxJm18RH0PGfk5FllPpK+vY/rj0xgrwpuM/FqLrD9azbyTse0c3dueNYNd5reVV1tkFxrZ7RbZKCObnLHdHPWrex+7b5T4VfAwci0amYfGVBxFmrLnoFnjM7Sq+hINEkeG8vVBP4UJnzfOAU7J0LcdyYGo7z6lUbaBFK7OLAqkjzDve1jq2xdFlFdiP18diraRaWfWeiYfJFoYIetJsd6qRTPYJRSCTzYdjoiR5Z/LJxZZpamV7kv8KjgQbgfeQ7/VO6zMipPKnoqujHwLnIkGzbvQ1jx8T2hPFEpfF0pfR2mktoPiiFOl2WLeD6H4rGUQOh89AG1/gwNX3gjHUHy1qBVFjFvQasm2bT7ZvJe7UuqgunpJQzvwN7pUHiavo2o7Y6XsEqr/DPL6XW6Rxem3PUZWSSqp+yRK/Cocnu4072dlqDyqbBu6sPkrigiuNenPo9XB+dgPPcNL2h6WtHz/0/5ao7tsAN4GDkdL9UnAMygaugkttz+nMGD2QmeEi5GRLkF3tx5Bs9/pwM0oSmbjTLRlfLnMflZbL0m0AEchJ/3bIo9bsVSKStolVPcZtKCzsmXY9duOzqA3Rsg2o4hztai07pNI9Ks+KKQ8P5TukXw1wlb2MJP2M8X3mvLkLzbOC9WzFbgolPdR4N1Q2mI0AO0e068dzV5oJbcB3RmciwIkx6PvMjWQt43CxdD90WH1JvN6k/hzlQHIQF9K6I9PqTHUQi8e0TaS/68fj0eU/QLpMuu9sRzlnVNV2i6hus8gTr/5WwmdFlkLCp7Ms8jSkqP+dJ/Hp9Q30voVt5vCwYu/HskDYVTZLMyn9Fb5CooPRgeiUX1SN9uqV24k3QVjn+KHXSu9eKSzkUYmjV1C89tmLUir+yA+pQNhWr9iV7SMnh1I8yi+HxT1v+ZsZbNwMVrij0VbqsnoZ0GDA3nOQ1vQLGcE9U5f4Hu0VbAR9z/XaqUXj3Q20siksUtobtusFWl1H+cbkX5lu+C7BbgM/YC7H/ptbBdwdyCP7awhqmwWZqIAxJ0osrocRZiCd+1mk+3OYyMwBM1+HRHyPyl+HkFqpZcu0tlII5PGLqG5bbNWpNV9nG8MId6vHA6Hw+FwOBwOh8PhcDgcDsdOyn9L6N59ZKcf2AAAAABJRU5ErkJggg==\n", + "text/latex": [ + "$$\\left ( \\left [ \\left ( x_{0}, \\quad a b\\right )\\right ], \\quad \\left [ n + x_{0}, \\quad m + x_{0}\\right ]\\right )$$" + ], + "text/plain": [ + "([(x₀, a⋅b)], [n + x₀, m + x₀])" + ] + }, + "execution_count": 24, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "cse(exprs, optimizations='basic')" + ] + }, + { + "cell_type": "code", + "execution_count": 25, + "metadata": { + "slideshow": { + "slide_type": "subslide" + } + }, "outputs": [], "source": [ "subexprs, exprs = cse(lbm, optimizations='basic')" @@ -898,7 +868,7 @@ }, { "cell_type": "code", - "execution_count": 25, + "execution_count": 26, "metadata": { "slideshow": { "slide_type": "fragment" @@ -908,24 +878,24 @@ { "data": { "text/latex": [ - "$\\displaystyle \\begin{align} x_{0} &= 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}, \\\\ x_{1} &= - f_{curr 3} + f_{curr 5}, \\\\ x_{2} &= - f_{curr 0} + f_{curr 8}, \\\\ x_{3} &= f_{curr 2} - f_{curr 6}, \\\\ x_{4} &= x_{1} + x_{2} + x_{3}, \\\\ x_{5} &= \\frac{6}{x_{0}}, \\\\ x_{6} &= x_{4} x_{5}, \\\\ x_{7} &= 2 f_{curr 0} \\end{align}$" + "$\\displaystyle \\begin{align} x_{0} &= 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}, \\\\ x_{1} &= - f_{curr 3} + f_{curr 5}, \\\\ x_{2} &= - f_{curr 0} + f_{curr 8}, \\\\ x_{3} &= f_{curr 2} - f_{curr 6}, \\\\ x_{4} &= x_{1} + x_{2} + x_{3}, \\\\ x_{5} &= \\frac{6}{x_{0}}, \\\\ x_{6} &= x_{4} x_{5}, \\\\ x_{7} &= 2 f_{curr 0}, \\\\ x_{8} &= 2 f_{curr 8}, \\\\ x_{9} &= - f_{curr 1} + f_{curr 7}, \\\\ x_{10} &= \\frac{1}{x_{0}^{2}}, \\\\ x_{11} &= 9 x_{10}, \\\\ x_{12} &= - f_{curr 2} + f_{curr 6}, \\\\ x_{13} &= x_{12} + x_{2} + x_{9}, \\\\ x_{14} &= x_{4}^{2}, \\\\ x_{15} &= 3 x_{10}, \\\\ x_{16} &= x_{13}^{2} x_{15} - 2, \\\\ x_{17} &= x_{14} x_{15} + x_{16}, \\\\ x_{18} &= \\frac{1}{\\tau}, \\\\ x_{19} &= \\frac{x_{18}}{72}, \\\\ x_{20} &= f_{curr 0} - f_{curr 8}, \\\\ x_{21} &= f_{curr 1} - f_{curr 7}, \\\\ x_{22} &= x_{20} + x_{21} + x_{3}, \\\\ x_{23} &= x_{22} x_{5}, \\\\ x_{24} &= x_{22}^{2}, \\\\ x_{25} &= 6 x_{10}, \\\\ x_{26} &= f_{curr 3} - f_{curr 5}, \\\\ x_{27} &= x_{12} + x_{20} + x_{26}, \\\\ x_{28} &= x_{27}^{2}, \\\\ x_{29} &= - x_{15} x_{28} + 2, \\\\ x_{30} &= x_{24} x_{25} + x_{29}, \\\\ x_{31} &= \\frac{x_{18}}{18}, \\\\ x_{32} &= x_{27} x_{5}, \\\\ x_{33} &= - x_{15} x_{24}, \\\\ x_{34} &= - x_{32} + x_{33}, \\\\ x_{35} &= x_{11} \\left(- 2 f_{curr 2} + 2 f_{curr 6} + x_{26} + x_{9}\\right)^{2} + x_{29}, \\\\ x_{36} &= - x_{23} \\end{align}$" ], "text/plain": [ "" ] }, - "execution_count": 25, + "execution_count": 26, "metadata": {}, "output_type": "execute_result" } ], "source": [ - "aligned(subexprs[0:8])" + "aligned(subexprs)" ] }, { "cell_type": "code", - "execution_count": 26, + "execution_count": 27, "metadata": { "slideshow": { "slide_type": "subslide" @@ -941,7 +911,7 @@ "" ] }, - "execution_count": 26, + "execution_count": 27, "metadata": {}, "output_type": "execute_result" } @@ -963,7 +933,7 @@ }, { "cell_type": "code", - "execution_count": 27, + "execution_count": 28, "metadata": { "slideshow": { "slide_type": "fragment" @@ -971,27 +941,20 @@ }, "outputs": [ { - "data": { - "image/png": "iVBORw0KGgoAAAANSUhEUgAAACcAAAASCAYAAADYFMcrAAAABHNCSVQICAgIfAhkiAAAAjJJREFUSInt1cuLiGEUBvCf65RJuVNIjWFDiYVLyhgLJZKVhWRlYyFKCRs2CJMoSbMgpRRiRiSXWJH8BYYao+QyLsndMDMW75lmvPN9GFbkqa+37zzn+l7O4S/DUlzBI3xEM05jXoFuCzpLvqc/iDEBR/EYn8PPAQzvqdQvM9qDzXiJBrxANZZjINbgRJbcsHCc4x3qCuSTcQtj0Ii7mI1aNGF+xP8O49AuVTwm42ql3WjO5C3x9QWXw9f6TL4/5EeKjOYE2Vji9A3e/mFyVRHjAfpn3FBpt9+jUqZwH23SFo/KDBeE8bWCgBVYjW3YIO3ygJLkFsV6BR0Z9xY3MQRzi4w3hlEr6rEbp/ApHObH3aL4MTSjpsD/vuA3lSR/KPh1JbwVeJUFu49VBbrbpd0YK1U8XbozHfiAGZl+ffhbWxJ7Z/Bbi8jN+CpdzqoIOEv3Jd5bVlGGutA/18fkdgW/JScWBnG2wGiI1PfaI+mfoTp85S2hT8fa80Esi/VGgdEH3An9mb+QXGuslZm8KdapJXZTYr2XJ1cR6+gSwy552y8k1zVN8r7YVfhixa1kvjSVbucOV+oeO+Mzbol0yT9iZMimYURBYpOkB9QptZccv9WE++NqKLzBcWmcnY/EOqU+1oUdUou5hMOheyYK6MRFDC6IMxnPQqdBalfX479Jd/G9MEjqdbcjwa/S/bkgHUVP1OCkNBtf4wueR4Fr9J7bPTERx/BEuiYPcVDxSfzHv4VvYwymLpcPWpcAAAAASUVORK5CYII=\n", - "text/latex": [ - "$$850$$" - ], - "text/plain": [ - "850" - ] - }, - "execution_count": 27, - "metadata": {}, - "output_type": "execute_result" + "name": "stdout", + "output_type": "stream", + "text": [ + "476*ADD + 60*DIV + 9*EQ + 118*MUL + NEG + 40*POW + 146*SUB = 850\n" + ] } ], "source": [ - "count_ops(lbm)" + "print('%s = %s' % (count_ops(lbm, visual=True), count_ops(lbm)))" ] }, { "cell_type": "code", - "execution_count": 28, + "execution_count": 29, "metadata": { "slideshow": { "slide_type": "fragment" @@ -999,22 +962,15 @@ }, "outputs": [ { - "data": { - "image/png": "iVBORw0KGgoAAAANSUhEUgAAACcAAAASCAYAAADYFMcrAAAABHNCSVQICAgIfAhkiAAAAdZJREFUSInt1k+ITXEUB/DPTE+KlGKkJEUWs0RNKemZpERSZMWOrMhCyEKz9L+RLXqFlWlYWKEUCjsrkj9Nspg0JLOYTMazOPfldufd++6dsZnyrVu38+/3vef8zjmXOYS9uIqn+IEmblXwP5D4NHGwwG4HHuAzJvARd7CxKPirJPA43lQktxLfE98icucS/Riu4SyGMInf2J93wBasRRfqFch14RE+4EIBueWYwiiWtTm7KbIIahmDxyWItMNR9IsP6i+wW4VuvMSXNmePo6cl6J4hmTR6RWmu4EkH23eifH1YmtFtxiJRAUzPXFXUcBOfcLqE/TecxGW8xj18xRrswkMc/lfkzmAdNomuK4NBjOAGDqXk79GQKvdsytonsnUJzyv4nRDd2RAZW4gNohFu43yZIHX53VrDW1Ga+RndgPxubcUcbqNbIObeFFbPhtxif4dtp2cw5XcxkR3JOXM40e9h5nfuJ67n6NaLe/hMZDdd8laWe7JOGflkJwJ11dcXxWXdl+hGsSKj2y42xASWMD1zu5OHmObEvmsk72M4XpFsGkNijm0V6/FuQrQXO8WmOSXGyzQMKL4/IyUItGLk7dZ5OIYX4ufilxgf97GtRPz/mNv4A9lIhKfVHIfDAAAAAElFTkSuQmCC\n", - "text/latex": [ - "$$148$$" - ], - "text/plain": [ - "148" - ] - }, - "execution_count": 28, - "metadata": {}, - "output_type": "execute_result" + "name": "stdout", + "output_type": "stream", + "text": [ + "43*ADD + 5*DIV + 9*EQ + 49*MUL + 3*NEG + 8*POW + 31*SUB = 148\n" + ] } ], "source": [ - "count_ops([subexprs, exprs])" + "print('%s = %s' % (count_ops([subexprs, exprs], visual=True), count_ops([subexprs, exprs])))" ] }, { @@ -1047,7 +1003,7 @@ }, { "cell_type": "code", - "execution_count": 29, + "execution_count": 30, "metadata": { "slideshow": { "slide_type": "subslide" @@ -1105,7 +1061,7 @@ }, { "cell_type": "code", - "execution_count": 30, + "execution_count": 31, "metadata": { "slideshow": { "slide_type": "subslide" @@ -1152,8 +1108,6 @@ } }, "source": [ - "## Measuring performance\n", - "\n", "$$\\text{MLUPS} := \\frac{n_\\text{cells} \\ n_\\text{updates}}{10^6 \\, \\Delta t}$$\n" ] }, @@ -1161,7 +1115,7 @@ "cell_type": "markdown", "metadata": { "slideshow": { - "slide_type": "fragment" + "slide_type": "skip" } }, "source": [ @@ -1178,22 +1132,17 @@ } }, "source": [ - "| GPU | Bandwidth | D2Q9 |   | D3Q19 |   | D3Q27 |   | \n", - "| ------ | ----------- | ------ | ------ | ------ | ------ | ------ | ------ |\n", - "|   |   | single | double | single | double | single | double | \n", - "| P100 | 512.6 GiB/s | 7242 | 3719 | 3528 | 1787 | 2502 | 1262 |" + "### Measured performance..." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { - "slide_type": "subslide" + "slide_type": "fragment" } }, "source": [ - "### Measured performance...\n", - "\n", "| CSE | D2Q9 |   | D3Q19 |   | D3Q27 |   |\n", "| ------ | ------ | ------ | ------ | ------ | ------ | ------ |\n", "|   | single | double | single | double | single | double |\n", @@ -1209,8 +1158,31 @@ } }, "source": [ - "### …compared to theoretical maximum\n", - "\n", + "### …compared to theoretical maximum" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "fragment" + } + }, + "source": [ + "| GPU | Bandwidth | D2Q9 |   | D3Q19 |   | D3Q27 |   | \n", + "| ------ | ----------- | ------ | ------ | ------ | ------ | ------ | ------ |\n", + "|   |   | single | double | single | double | single | double | \n", + "| P100 | 512.6 GiB/s | 7242 | 3719 | 3528 | 1787 | 2502 | 1262 |" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "fragment" + } + }, + "source": [ "| CSE | D2Q9 |   | D3Q19 |   | D3Q27 |   |\n", "| ------ | ------ | ------ | ------ | ------ | ------ | ------ |\n", "|   | single | double | single | double | single | double |\n", @@ -1220,7 +1192,7 @@ }, { "cell_type": "code", - "execution_count": 31, + "execution_count": 32, "metadata": { "slideshow": { "slide_type": "skip" @@ -1275,7 +1247,7 @@ " plt.scatter(*zip(*list(data)), **kwargs)\n", " \n", "def mlups_per_size_overview_plot(data, title, **kwargs):\n", - " plt.figure(figsize=(11,7.5))\n", + " plt.figure(figsize=(12,8))\n", " plt.grid()\n", " plt.title(title)\n", " plt.xscale('log')\n", @@ -1291,7 +1263,7 @@ }, { "cell_type": "code", - "execution_count": 32, + "execution_count": 33, "metadata": { "slideshow": { "slide_type": "subslide" @@ -1300,9 +1272,9 @@ "outputs": [ { "data": { - "image/png": "iVBORw0KGgoAAAANSUhEUgAAAqUAAAHYCAYAAAB0jJGwAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzs3Xl8VeW59//PlRBASIhKIooCCREcGIwIiFoxSkEtaqPHGRUtWKMVbas8rcdTh1Zb+9BSHzl9pCooKI946hCHWg/IMag/FdJgQJRqCEkJgzIpJKASkvv3x1qJScgIO3vt4ft+vfLKvtd4rb0X5Nr3tMw5h4iIiIhIkBKCDkBEREREREmpiIiIiAROSamIiIiIBE5JqYiIiIgETkmpiIiIiAROSamIiIiIBE5JqYiIiIgETkmpSBwzswfMbJuZfR50LJHAzM4wsxIzqzKz3KDjkdAwszwzezPoOESkdUpKRaKImZWb2dd+0vSFmT1pZskHeKx+wB3Aic65I0MbadT6NfCfzrlk51x+05UN3v9KM/vKzN7zE56EBtv81MzWmdkuM9tkZn8ysy4N1puZTfeT36/NbL2Z/dbMujbYZrqZrfbPU2Zm0xus6+9//g1/nJnd0dwFmdl9ZvZMC+vavB5/u9Fm9rq/zQ4zW25mN7T2RprZYf6XntX+PuvM7DEzG9jKPrMbXNNeM6tuUP57a+frLGa20My+9WPYYWZvmNkgf90IM1tkZtvN7Jtm9k03s1fNbLf/OV7aZP1k//OvMrPnzSw1XNclEomUlIpEnwudc8nACGAU8B8dPYCfJA0Atjvnthzg/rFoAPBxG9tc6JxL8bd9CPgFMKfB+leBEc65XsBQ4CTgtgbrHwF+DFwHpADnA98HFjbYxvz1hwHnAbea2ZUAzrn1ftKc7N8Hw4Ba4IWOX27b12NmpwH/AywFjgV6Azf7cTfLzI4HlgNdgH8D0oFTgPeBRWY2obn9nHN5Da7rt8BzDa61xfOFwW/8mPoDlcDj/vJvgWeBvBb2ewz4EjgCmALMbZDQnox3L1wBHIX3mf+fzroAkWigpFQkSjnnNgJ/x0t8MLNUM5tjZpvNbKNfS5Xor7vezP4/v9ZuB1AALAb6+rU0T/nbXWRmH/s1YgVmdkLd+fxatV+Y2Spgt5l18ZdNN7NVfm3QHDPrY2Z/92vf3jSzwxoc469m9rmZ7TSzt81sSIN1T5nZn83sb/6+y8wsq8H6IWa22K+t+sLM/t1fnmBmvzSzUr/G6r/M7PCW3jczu9HM1vrHecXM+vrLS4GBwKv+e9Ktjfd/p3PuFbykYrKZDfWXlzrnvqo7HV7CeKx/jkHALcAk59z7zrl9zrmP8RK3iWZ2ln+M/+2cW+Gv/xR4GTijhVCuA952zpW3Fm9bWroeYAYwzzn3e+fcNucpcs5d3txxzKvxfQG4xTn3S+fcp865Gufcl865J4FxwH+a2aEHEqeZnenfG1+Z2QozO6PBuhv9e7LSr5m9rIVjPGpmG8yrzV5uZmPac27nXBXel4e6z/pj/5rWNHOOw4ALgV8553Y75/4H+G9gkr/JtcAL/n1QCdwDXGFm3dv9ZojEGCWlIlHKvOb3HwAf+ovmAfvwEqCTgQnA1Aa7nAqsw6u1GY9X07XJr4W63swG49X6/BSvZut1vASta4NjXAVMBA51zu3zl/2bf7zBeH+E/w78O5CG939Mw1rCvwOD/BhWAAuaXNZVwP14NYRrgQf9a00B3gTeAPr617jE3+c2IBc4y1/3JfDnFt6zc4DfAZfj1U79C7+G0jmXBazHr4l2zn3b3DGacs4tBzYAZzY4z9VmtgvYhldT+hd/1Thgg79Pw2NUAB/gfWZNYzb/2C3V4F6H99mHRMPrMbMewGnA8x04xFXAu865xWY2zMwKzazCvG4E7znn/uXHe01HYzOzDCAfuBs4HK+VIN+8rgKH4SXQ4/ya3zOB1S0c6n28GubeeAn/X80sqR3n7+Vf34dtbQscD1T611tnJVD3RWyIXwa8BBdIBLIQiVNKSkWiT76ZfQW8i9ek+lsz64OXZP7Ur5XZAvwJuLLBfpucc7P82revmznuFcDfnHOLnXPVwB+AQ4DTG2zziHOuosn+s5xzX/g1t+8Ay5xzH/pJ3Ut4CTIAzrm5zrlKf919wElN+tG96Jxb7ie8C4Bsf/kFwOfOuT86577xj7HMX3cTcLdzbkOD417aQheDScBcvxbyW+Au4DQ/2TkYm/CSpLrr/H9+8/1gYDbwhb8qDdjcwjE2430ZaOo+vP+rn2y6wszOBPrQsaSxPequ5zD/3C3F3JzxfNcV4Qm8pu5M/5h9/eXFeElbR03Gu0fedM7VOudeBz6hcTI/1My6O+c2Ouf2q8EEcM7N92tuq/G6CfTGqyVvyd3+v7lP8RLHqa1sWycZ2Nlk2U68Lhstra9ssF4k7igpFYk+uc65Q51zA5xzt/gJ4gAgCdjsN2t+hVc7d0SD/SraOG5fvJpDAJxztf4+R7dxjC8avP66mXIygJklmtlDfjP7LqDc3yatwfYNZwHYU7cv0A8obSHuAcBLDa57DVCDl6w11fQaq4DtNL7GA3E0sKPpQudcCV4N5//1F23Dq6FtzlHA1oYLzOxWvJrQiS3U3E7GawKuOsC4W1J3PV/idT9oKebmHAFs9F8PA57xv2Q0HGzVr8E2HTEAuKbus/Y/75FAX+fcl3hfOm4DPve7Zhzb3EHM7C4z+9TMduJdY3ca34dNPej/mzvKOXdxk9rPllQBvZos64WXeLa0PqXBepG4o6RUJDZU4A26SPP/eB7qnOvlnBvSYBvXxjE24f3RB+qbjZsmD20dozVXAz/EG9STCmTUnaod+1bQcrNmBXB+g+s+1DnX3a+5barpNfbEqyU7kASp7hij8JK4d1vYpEuD2P8H6Gdmo5scox8wBq/mu27Zj4Bf4jVHb2jmvIcAlxHCpnv/uPXX45zbg9fU/W8dOETDxPsjvCQyEb+53sxOAaYB/+8AwqsAnmjyWfd0zv0JwDn3N+fcOLwvH+uBR5u5vvH++S8GDsWrEf6a9t2HHfFPoJeZ9W+w7CS+64bxsV+ui+tEvC9TLX35Eol5SkpFYoBzbjOwCPijmfXyB/9k1Q2caaf/whtsM87vX3cHXqL7XojCTPGPtx3ogdds2l6vAUeaN91SNzNLMbNT/XWzgQfNbADUT8PzwxaO8/+AG8ws2x/I9Fu87gblHb0Y/32+AK+p+hnn3Ef+8qlmdoT/+kS8LgJLAJxzn/nxLjCzMX7t8RC8gUHv4fWbxcwm+bGNd86tayGEi4GvgLfaEW6CmXVv8LPfIK6Wrgf4X8D15g1o6+1ve5KZLWx6DN//AHVTH00FbsSrnT4W2A38Bri2nbWNTc0DLvPv0UQzO8R/faSZHW1mE/1+sN/i1UTWNHOMFKAar1a6K940YAc0uMg83f3j4L+3XQH8mtvXgF+bWQ8zy8GbSaGuH/UzwL/590EyXl/q55xz+00tJRIvlJSKxI7r8P44foLXJPk8HWh29Ud5XwPMwqvtuhBv0M/eEMU3Hy852ejH+EEHYqvE66t4IV4Tfwlwtr/6/wCv4E01VOkf99QWjrME+BVeErgZrwbzyua2bcWr/nkq8AbczAQaztl5BvCRme3GGyz2Ot7Arzq34vW1fAavi8JqvPcl1+8yAfAAXg1uoX03T+fsJnFMBuY759pTe30VXm1g3U/D2rhWr8c59x5wjv+zzrzZGx7zr6s5zwDjzews59xHzrlRzrljnHP/y6+5v8g5t6IdMe/HT9D/DS+B24b3vt2O97csEe8LwOd4X3xG4dWINvUq8Lb/Hqzzj7O1me3a4zi897MI6Oa/XtVg/Y143QK2AU8BU/wuHfjvwU/x/p1+4V/D7QcYh0hMsPb9fyYiIp3BzH6NN3vAWPfdVFJRzcyG4Y1qfwyvZnAj3mCnu4Ba59yPAwxPRCKUklIRkYD5A5rWOufeCDqWUPFnhPgl3hRiffCmmXoO+KNzbneQsYlIZFJSKiIiIiKBU59SEREREQlcTD6/Oi0tzWVkZIT9vLt376Znz55hP69IW3RvSqTSvSmRSvdm6BQVFW1zzjX3cJBGYjIpzcjI4B//+EfYz1tQUEBOTk7YzyvSFt2bEql0b0qk0r0ZOmbWring1HwvIiIiIoFTUioiIiIigVNSKiIiIiKBi8k+pc2prq5mw4YNfPNN5z3BLTU1lTVr1nTa8cXTvXt3jjnmGJKSkoIORUREREIkbpLSDRs2kJKSQkZGBmbWKeeorKwkJSWlU44tHucc27dvZ8OGDWRmZgYdjoiIiIRI3DTff/PNN/Tu3bvTElIJDzOjd+/enVrjLSIiIuEXN0kpoIQ0RuhzFBERiT1xlZSKiIiISGRSUhqwqVOn8sknnxzQvuXl5QwdOvSAz11VVcVNN91EVlYWQ4YMYezYsSxbtgyABx98kCFDhjB8+HCys7Prl+fk5HDccceRnZ1NdnY2l1566QGfX0RERKRO3Ax0ilRPPPFEYOeeOnUqmZmZlJSUkJCQwLp161izZg3vv/8+r732GitWrKBbt25s27aNvXv31u+3YMECRo4cGVjcIiIiEntUUxomu3fvZuLEiZx00kkMHTqU5557DvBqHuseiZqcnMzdd9/NSSedxJgxY/jiiy8AKC0tZcyYMYwaNYp77rmH5OTk/Y5fU1PD9OnTGTVqFMOHD+cvf/lLq/GUlpaybNkyHnjgARISvNtg4MCBTJw4kc2bN5OWlka3bt0ASEtLo2/fviF7L0RERESaUlLagvXb9zB+5lKy7nqd8TOXsn77noM63htvvEHfvn1ZuXIlq1ev5rzzzttvm927dzNmzBhWrlzJ2LFjefzxxwG4/fbbuf322yksLGwxOZw