diff options
| -rw-r--r-- | notebook/seminar_talk.ipynb | 454 |
1 files changed, 213 insertions, 241 deletions
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,6 +1,28 @@ { "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": { @@ -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", @@ -204,6 +237,22 @@ } }, "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" ] }, @@ -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}$$" ] @@ -339,34 +388,6 @@ ] }, { - "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": { @@ -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" } @@ -545,8 +566,19 @@ ] }, { + "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 @@ "<IPython.core.display.Math object>" ] }, - "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 @@ "<IPython.core.display.Math object>" ] }, - "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])" ] }, { @@ -885,12 +810,57 @@ }, { "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, "metadata": { "slideshow": { "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/WohWLOFAh0fyt |