08_tutorial_shanchen_twocomponent.ipynb 42.8 KB
 Michael Kuron committed Oct 17, 2019 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 { "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Shan-Chen Two-Component Lattice Boltzmann" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "from lbmpy.session import *\n", "from lbmpy.updatekernels import create_stream_pull_with_output_kernel\n",  Martin Bauer committed Oct 21, 2019 18 19  "from lbmpy.macroscopic_value_kernels import macroscopic_values_getter, macroscopic_values_setter\n", "from lbmpy.maxwellian_equilibrium import get_weights"  Michael Kuron committed Oct 17, 2019 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42  ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This is based on section 9.3.3 of Krüger et al.'s \"The Lattice Boltzmann Method\", Springer 2017 (http://www.lbmbook.com).\n", "Sample code is available at [https://github.com/lbm-principles-practice/code/](https://github.com/lbm-principles-practice/code/blob/master/chapter9/shanchen.cpp)." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Parameters" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [  Martin Bauer committed Oct 21, 2019 43 44 45 46 47  "N = 64 # domain size\n", "omega_a = 1. # relaxation rate of first component\n", "omega_b = 1. # relaxation rate of second component\n", "\n", "# interaction strength\n",  Michael Kuron committed Oct 17, 2019 48 49 50  "g_aa = 0.\n", "g_ab = g_ba = 6.\n", "g_bb = 0.\n",  Martin Bauer committed Oct 21, 2019 51 52 53 54 55  "\n", "rho0 = 1.\n", "\n", "stencil = get_stencil(\"D2Q9\")\n", "weights = get_weights(stencil, c_s_sq=sp.Rational(1,3))"  Michael Kuron committed Oct 17, 2019 56 57 58 59 60 61  ] }, { "cell_type": "markdown", "metadata": {}, "source": [  Martin Bauer committed Oct 21, 2019 62 63 64 65 66  "## Data structures\n", "\n", "We allocate two sets of PDF's, one for each phase. Additionally, for each phase there is one field to store its density and velocity.\n", "\n", "To run the simulation on GPU, change the default_target to gpu"  Michael Kuron committed Oct 17, 2019 67 68 69 70 71 72 73 74  ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [  Martin Bauer committed Oct 21, 2019 75 76  "dim = len(stencil[0])\n", "dh = ps.create_data_handling((N, ) * dim, periodicity=True, default_target='cpu')\n",  Michael Kuron committed Oct 17, 2019 77  "\n",  Martin Bauer committed Oct 21, 2019 78  "src_a = dh.add_array('src_a', values_per_cell=len(stencil))\n",  Michael Kuron committed Oct 17, 2019 79 80  "dst_a = dh.add_array_like('dst_a', 'src_a')\n", "\n",  Martin Bauer committed Oct 21, 2019 81  "src_b = dh.add_array('src_b', values_per_cell=len(stencil))\n",  Michael Kuron committed Oct 17, 2019 82 83 84 85 86 87 88 89 90 91 92 93  "dst_b = dh.add_array_like('dst_b', 'src_b')\n", "\n", "ρ_a = dh.add_array('rho_a')\n", "ρ_b = dh.add_array('rho_b')\n", "u_a = dh.add_array('u_a', values_per_cell=dh.dim)\n", "u_b = dh.add_array('u_b', values_per_cell=dh.dim)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [  Martin Bauer committed Oct 21, 2019 94 95 96 97  "## Force & combined velocity\n", "\n", "The two LB methods are coupled using a force term. Its symbolic representation is created in the next cells.\n", "The force value is not written to a field, but directly evaluated inside the collision kernel."  Michael Kuron committed Oct 17, 2019 98 99 100 101 102 103 104  ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The force between the two components is\n",  Michael Kuron committed Dec 19, 2019 105  "$\\vec{F}_k(\\vec{x})=-\\psi(\\rho_k(\\vec{x}))\\sum\\limits_{k^\\prime\\in\\{A,B\\}}g_{kk^\\prime}\\sum\\limits_{i=1}^{q}w_i\\psi(\\rho_{k^\\prime}(\\vec{x}+\\vec{c}_i))\\vec{c}_i$\n",  Michael Kuron committed Oct 17, 2019 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173  "for $k\\in\\{A,B\\}$\n", "and with \n", "$\\psi(\\rho)=\\rho_0\\left[1-\\exp(-\\rho/\\rho_0)\\right]$." ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "def psi(dens):\n", " return rho0 * (1. - sp.exp(-dens / rho0));" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [], "source": [ "zero_vec = sp.Matrix([0] * dh.dim) \n", "\n", "force_a = zero_vec\n", "for factor, ρ in zip([g_aa, g_ab], [ρ_a, ρ_b]):\n", " force_a += sum((psi(ρ[d]) * w_d * sp.Matrix(d)\n", " for d, w_d in zip(stencil, weights)), \n", " zero_vec) * psi(ρ_a.center) * -1 * factor\n", "\n", "force_b = zero_vec\n", "for factor, ρ in zip([g_ba, g_bb], [ρ_a, ρ_b]):\n", " force_b += sum((psi(ρ[d]) * w_d * sp.Matrix(d)\n", " for d, w_d in zip(stencil, weights)), \n", " zero_vec) * psi(ρ_b.center) * -1 * factor" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The barycentric velocity, which is used in place of the individual components' velocities in the equilibrium distribution and Guo force term, is\n", "$\\vec{u}=\\frac{1}{\\rho_a+\\rho_b}\\left(\\rho_a\\vec{u}_a+\\frac{1}{2}\\vec{F}_a+\\rho_b\\vec{u}_b+\\frac{1}{2}\\vec{F}_b\\right)$." ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [], "source": [ "u_full = [(ρ_a.center * u_a(i) + force_a[i]/2 + \n", " ρ_b.center * u_b(i) + force_b[i]/2) / (ρ_a.center + ρ_b.center)\n", " for i in range(dh.dim)]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Kernels" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [], "source": [  Martin Bauer committed Oct 21, 2019 174 175  "collision_a = create_lb_update_rule(stencil=stencil,\n", " relaxation_rate=omega_a, \n",  Michael Kuron committed Oct 17, 2019 176 177 178 179 180 181 182  " compressible=True,\n", " velocity_input=u_full, density_input=ρ_a,\n", " force_model='guo', \n", " force=force_a,\n", " kernel_type='collide_only',\n", " optimization={'symbolic_field': src_a})\n", "\n",  Martin Bauer committed Oct 21, 2019 183 184  "collision_b = create_lb_update_rule(stencil=stencil,\n", " relaxation_rate=omega_b,\n",  Michael Kuron committed Oct 17, 2019 185 186 187 188 189  " compressible=True,\n", " velocity_input=u_full, density_input=ρ_b,\n", " force_model='guo', \n", " force=force_b,\n", " kernel_type='collide_only',\n",  Martin Bauer committed Oct 21, 2019 190 191 192 193 194 195 196 197 198 199 200 201 202  " optimization={'symbolic_field': src_b})" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [], "source": [ "stream_a = create_stream_pull_with_output_kernel(collision_a.method, src_a, dst_a, \n", " {'density': ρ_a, 'velocity': u_a})\n", "stream_b = create_stream_pull_with_output_kernel(collision_b.method, src_b, dst_b, \n", " {'density': ρ_b, 'velocity': u_b})\n",  Michael Kuron committed Oct 17, 2019 203  "\n",  Martin Bauer committed Oct 21, 2019 204 205  "opts = {'cpu_openmp': 1, # number of threads when running on CPU\n", " 'target': dh.default_target}\n",  Michael Kuron committed Oct 17, 2019 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220  "stream_a_kernel = ps.create_kernel(stream_a, **opts).compile()\n", "stream_b_kernel = ps.create_kernel(stream_b, **opts).compile()\n", "collision_a_kernel = ps.create_kernel(collision_a, **opts).compile()\n", "collision_b_kernel = ps.create_kernel(collision_b, **opts).compile()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Initialization" ] }, { "cell_type": "code",  Martin Bauer committed Oct 21, 2019 221  "execution_count": 9,  Michael Kuron committed Oct 17, 2019 222 223 224  "metadata": {}, "outputs": [], "source": [  Martin Bauer committed Oct 21, 2019 225  "init_a = macroscopic_values_setter(collision_a.method, velocity=(0, 0), \n",  Michael Kuron committed Oct 17, 2019 226  " pdfs=src_a.center_vector, density=ρ_a.center)\n",  Martin Bauer committed Oct 21, 2019 227  "init_b = macroscopic_values_setter(collision_b.method, velocity=(0, 0), \n",  Michael Kuron committed Oct 17, 2019 228 229 230 231 232 233 234  " pdfs=src_b.center_vector, density=ρ_b.center)\n", "init_a_kernel = ps.create_kernel(init_a, ghost_layers=0).compile()\n", "init_b_kernel = ps.create_kernel(init_b, ghost_layers=0).compile()" ] }, { "cell_type": "code",  Martin Bauer committed Oct 21, 2019 235  "execution_count": 10,  Michael Kuron committed Oct 17, 2019 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261  "metadata": {}, "outputs": [], "source": [ "def init():\n", " dh.fill(ρ_a.name, 0.1, slice_obj=ps.make_slice[:, :0.5])\n", " dh.fill(ρ_a.name, 0.9, slice_obj=ps.make_slice[:, 0.5:])\n", "\n", " dh.fill(ρ_b.name, 0.9, slice_obj=ps.make_slice[:, :0.5])\n", " dh.fill(ρ_b.name, 0.1, slice_obj=ps.make_slice[:, 0.5:])\n", "\n", " dh.run_kernel(init_a_kernel)\n", " dh.run_kernel(init_b_kernel)\n", " \n", " dh.fill(u_a.name, 0.0)\n", " dh.fill(u_b.name, 0.0)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Timeloop" ] }, { "cell_type": "code",  Martin Bauer committed Oct 21, 2019 262  "execution_count": 11,  Michael Kuron committed Oct 17, 2019 263 264 265 266 267 268 269 270 271  "metadata": {}, "outputs": [], "source": [ "sync_pdfs = dh.synchronization_function([src_a.name, src_b.name])\n", "sync_ρs = dh.synchronization_function([ρ_a.name, ρ_b.name])\n", "\n", "def time_loop(steps):\n", " dh.all_to_gpu()\n", " for i in range(steps):\n",  Martin Bauer committed Oct 21, 2019 272  " sync_ρs() # collision kernel evaluates force values, that depend on neighboring ρ's\n",  Michael Kuron committed Oct 17, 2019 273 274 275 276 277 278 279 280 281 282 283 284 285 286  " dh.run_kernel(collision_a_kernel)\n", " dh.run_kernel(collision_b_kernel)\n", " \n", " sync_pdfs()\n", " dh.run_kernel(stream_a_kernel)\n", " dh.run_kernel(stream_b_kernel)\n", " \n", " dh.swap(src_a.name, dst_a.name)\n", " dh.swap(src_b.name, dst_b.name)\n", " dh.all_to_cpu()" ] }, { "cell_type": "code",  Martin Bauer committed Oct 21, 2019 287  "execution_count": 12,  Michael Kuron committed Oct 17, 2019 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311  "metadata": {}, "outputs": [], "source": [ "def plot_ρs():\n", " plt.subplot(1,2,1)\n", " plt.title(\"$\\\\rho_A$\")\n", " plt.scalar_field(dh.gather_array(ρ_a.name), vmin=0, vmax=2)\n", " plt.colorbar()\n", " plt.subplot(1,2,2)\n", " plt.title(\"$\\\\rho_B$\")\n", " plt.scalar_field(dh.gather_array(ρ_b.name), vmin=0, vmax=2)\n", " plt.colorbar()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Run the simulation\n", "### Initial state" ] }, { "cell_type": "code",  Martin Bauer committed Oct 21, 2019 312  "execution_count": 13,  Michael Kuron committed Oct 17, 2019 313 314 315 316  "metadata": {}, "outputs": [ { "data": {  Martin Bauer committed Oct 21, 2019 317  "image/png": "iVBORw0KGgoAAAANSUhEUgAAA48AAAF0CAYAAACDowz8AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3dfbBkd33f+fdHI4lxBqwHBFglCSMS2RYkSOApGZdS5skImRCk1MKuFOJMKLlUScDBXu8mwruBXWxXmc2WiVMBY8UoklMGQYQxsy4ZoRWwJLHBGoECegBrkFk0HkVjMRJgGUmrme/+0WdM62rudM+9t/t2f+f9qjrVfR66+3fgSh99z+/8fidVhSRJkiRJR3LcZjdAkiRJkrT4LB4lSZIkSRNZPEqSJEmSJrJ4lCRJkiRNZPEoSZIkSZrI4lGSJEmSNJHFoyRJkiRpIotHSZIkSdJEFo9qL8kzklyd5KEk+5L8/Ga3SZKkY5nZLC0ni0cdC34P+BrwA8BlwP+Z5Ac2t0mSJB3TzGZpCVk8qrUkrwOoqndX1WNV9Sngz4AfGva/MMmBJGduZjslSTpWHCmbk/zHJP9lWD6VZOvmtlbSOItHdfd64OOHVpIcB5wEPDBs+hfAfwDOnX/TJEk6Jh0pm38IeFlVXQh8G3jRprRQ0mFZPKq7HwO+Obb+SuDBqvpqkhcB9wM3YfEoSdK8HDabgT8FjquqJ5I8HTgd+JNNaJ+kVVg8qq0kJwDnAG9IsjXJC4H3MeptBPh54N3AXVg8SpI0cxOy+VzgB5J8BrgHeG9VPbxpjZX0FMdvdgOkGToX+DpwB6NbYfYBv1xVNyQ5H7gQ+CCwZVgkSdJsHSmb3wT8RlW9I8lJwOeA3960lkp6CotHdfYi4O6q+pfAv1yx7yrgx6rqIYAkfzzvxkmSdAw6Ujb/TeDW4f0pwLfm2TBJk1k8qrPzgLtXbkzyo8B3DxWOg0eTPLOqvrnyeEmStGEOm82DFwIXJXkLcBD4x3NrlaSpOOZRnb0I+MrKjVV1W1W9ecW2n7BwlI4syVlJPp3k7iR3JnnbYY5Jkn+TZHeSLyV5ydi+HUnuGZYd8229pAVx2GwGqKrXV9WPVtWrqurVVXX7nNsmLZ15Z3OqaqPPQZLUUJLTgdOr6gtJngHcBlxaVXeNHfNa4GeB1zKaUfHXq+rHkpwK7AK2AzV89kdX3AEgSZKOwryz2Z5HSdJUqur+qvrC8P47jG49O2PFYZcAv10jnwNOHoLtNcDNVbV/CKWbgYvn2HxJktqZdzZbPEqSjlqS5wEvBj6/YtcZwH1j63uGbattlyRJG2Ae2TzXCXO2bNtWJ5x86jx/UpIW0v/38H4OPPJINvI7X/OKbfXN/QfW/PnbvvTYncCjY5uurqqrVx43PLz7o8DPVdW3V+4+zFfXEbZrk5143Nb6vuOesdnNkKRN992D3+Hxg4+azUcw1+LxhJNP5ax/+vPz/ElJWkj3ve89G/6dD+4/wOdvOnPNnz/h9K89WlXbj3TM8IDvjwK/U1W/e5hD9gBnja2fCewdtr98xfbPrLmx2jDfd9wz+PGT/t5mN0OSNt0ffetjG/6d3bLZ21YlSVNJEuADjJ7R9murHLYT+IfDzG4vBb5VVfcDNzGagv+UJKcAFw3bJEnSGs07m33OoyS1URyog7P8gQuBnwa+nOTQFPq/CDwXoKreD9zIaDa33cBfAm8e9u1P8kt87wHg76qq/bNsrCRJm69XNls8SlITBRyc4TDCqvrPHH58xPgxBbxllX3XANfMoGmSJC2kbtls8ShJjRxkplc3JUnSUeqUzY55lCRJkiRNZM+jJDVRFAfKp19IkrQoumWzxaMkNTLLcRWSJOnodcpmi0dJaqKAA40CSpKkZdctmy0eJamRTlc3JUnqoFM2O2GOJEmSJGkiex4lqYmCVoPyJUladt2y2eJRkhrp8yQpSZJ66JTNU922muTkJDck+UqSu5P8eJJTk9yc5J7h9ZRZN1aStLqiOLCORcvFbJakxdctm6cd8/jrwCeq6keA84C7gauAW6rqHOCWYV2StFkKDqxj0dIxmyVp0TXL5onFY5LvB34C+ABAVT1eVQ8DlwDXDYddB1w6q0ZKkqTvMZslSZthmp7H5wN/Dvz7JF9M8ltJtgHPqar7AYbXZx/uw0muTLIrya4DjzyyYQ2XJD1ZMRpXsdZFS2XDsvnxenR+rZakY0y3bJ6meDweeAnwG1X1YuARjuI2mKq6uqq2V9X2Ldu2rbGZkqTJwoF1LFoqG5bNJ2brrNooSWqWzdMUj3uAPVX1+WH9BkaB9UCS0wGG132zaaIkaRoFHKy1L1oqZrMkLYFu2TyxeKyq/wbcl+SHh02vAu4CdgI7hm07gI/PpIWSJOlJzGZJ0maY9jmPPwv8TpITgXuBNzMqPD+S5ArgG8AbZ9NESdK0FvEWF82M2SxJS6BTNk9VPFbV7cD2w+x61cY2R5K0VkWvgNKRmc2StPi6ZfO0PY+SpCVwsPoElCRJHXTKZotHSWqi29VNSZKWXbdsnma2VUmSJEnSMc6eR0lqoggHvCYoSdLC6JbNFo+S1EincRWSJHXQKZstHiWpiW7jKiRJWnbdstniUZLaCAeqz60xkiQtv17Z3OdMJEmSJEkzY8+jJDVRwEGvCUqStDC6ZbPFoyQ10mlchSRJHXTKZotHSWqiqte4CkmSll23bO5zJpIkSZKkmbHnUZIaOdjo1hhJkjrolM0Wj5LUxOhZUt5QIknSouiWzRaPktTGbMdVJLkGeB2wr6r+5mH2/8/Am4bV44FzgWdV1f4kXwe+AxwAnqiq7TNrqCRJC6NXNvcpgyXpGHdoOvC1LlO4Frh41d+v+ldVdX5VnQ+8Hfh/qmr/2CGvGPZbOEqSjgndstniUZI0lar6LLB/4oEjlwMfmmFzJEk65s07my0eJamRA5U1L8BpSXaNLVeupQ1J/hqjq6AfHdtcwCeT3LbW75UkaRl1ymbHPEpSE0XWOyj/wQ26pfTvAv9lxW0xF1bV3iTPBm5O8pXhaqkkSW11y2aLR0lq5OBiPIj4MlbcFlNVe4fXfUk+BlwAWDxKktrrlM0LcSaSpPU7NB34WpeNkOQk4GXAx8e2bUvyjEPvgYuAOzbkByVJWmDdstmeR0nSVJJ8CHg5o/EXe4B3AicAVNX7h8P+HvDJqnpk7KPPAT6WBEa588Gq+sS82i1JUlfzzmaLR0lqovirwfWz+f6qy6c45lpG04aPb7sXOG82rZIkaXF1y2aLR0lqZMpnQkmSpDnplM0Wj5LURBUcWIxB+ZIkiX7Z3OdMJEmSJEkzY8+jJLURDjK7cRWSJOlo9cpmi0dJaqLodWuMJEnLrls2WzxKUiMb9UwoSZK0MTpls8WjJDVRhIMznA5ckiQdnW7Z3KcMliRJkiTNjD2PktRIp1tjJEnqoFM2WzxKUhMFHGw0KF+SpGXXLZstHiWpjXCg0XTgkiQtv17ZbPEoSU10u7opSdKy65bNfc5EkiRJkjQz9jxKUiOdbo2RJKmDTtls8ShJTVSl1a0xkiQtu27ZbPEoSY0caBRQkiR10Cmb+5yJJEmSJGlm7HmUpCYKONhoXIUkScuuWzZbPEpSG2l1a4wkScuvVzZPVTwm+TrwHeAA8ERVbU9yKvBh4HnA14H/vqoemk0zJUmTjJ4l1efqpo7MbJakxdctm4+mDH5FVZ1fVduH9auAW6rqHOCWYV2StIkOcNyaFy0ls1mSFlynbF5Piy4BrhveXwdcuv7mSJKkdTCbJUkzM23xWMAnk9yW5Mph23Oq6n6A4fXZh/tgkiuT7Eqy68Ajj6y/xZKkwyrCwVr7oqWzIdn8eD06p+ZK0rGnWzZPO2HOhVW1N8mzgZuTfGXaH6iqq4GrAbaecVatoY2SpCkdXMBbXDQzG5LNJx3/LLNZkmaoUzZPVTxW1d7hdV+SjwEXAA8kOb2q7k9yOrBvhu2UJE1QBQcW8CqlZsNslqTF1y2bJ5bBSbYlecah98BFwB3ATmDHcNgO4OOzaqQkaTqdbo3R6sxmSVoenbJ5mp7H5wAfS3Lo+A9W1SeS3Ap8JMkVwDeAN86umZIkaYzZLEmau4nFY1XdC5x3mO3fBF41i0ZJko7eaFB+n3EVWp3ZLEnLoVs2TzthjiRpCRxg8W5xkSTpWNYpmy0eJamJgoUcHyFJ0rGqWzb36UOVJEmSJM2MPY+S1EavcRWSJC2/Xtnc50wkSRwka14mSXJNkn1J7lhl/8uTfCvJ7cPyjrF9Fyf5apLdSa7awFOWJGmhdcpmex4lqYk5PIj4WuDfAr99hGP+U1W9bnxDki3Ae4FXA3uAW5PsrKq7ZtVQSZIWQbdstniUpEZmeWtMVX02yfPW8NELgN3D4yVIcj1wCWDxKElqr1M2e9uqJOmQ05LsGluuXMN3/HiS/5rkD5K8cNh2BnDf2DF7hm2SJOnIFiqb7XmUpCZGDyJe160xD1bV9nV8/gvAD1bVXyR5LfB7wDlw2EEbtY7fkSRpKXTLZnseJamRWQ7Kn6Sqvl1VfzG8vxE4IclpjK5mnjV26JnA3nX/oCRJS6BTNtvzKElNbPaDiJP8APBAVVWSCxhdoPwm8DBwTpKzgT8DLgP+/qY1VJKkOemWzRaPkqSpJPkQ8HJG4y/2AO8ETgCoqvcDbwD+SZIngO8Cl1VVAU8keStwE7AFuKaq7tyEU5AkqZV5Z7PFoyQ1MuMZ3S6fsP/fMpou/HD7bgRunEW7JElaZJ2y2eJRkrqodQ/KlyRJG6lZNls8SlITBRsyuF6SJG2Mbtls8ShJjXS6uilJUgedstlHdUiSJEmSJrLnUZKa2OzpwCVJ0pN1y2aLR0lqpFNASZLUQadstniUpCaKXjO6SZK07Lpls8WjJDXSaUY3SZI66JTNTpgjSZIkSZrInkdJ6qJ6jauQJGnpNctmi0dJaqLbjG6SJC27btls8ShJjXQKKEmSOuiUzY55lCRJkiRNZM+jJDXRbTpwSZKWXbdstniUpEaqUUBJktRBp2y2eJSkRjo9S0qSpA46ZbPFoyQ1Uc2mA5ckadl1y2YnzJEkSZIkTWTPoyQ10mlchSRJHXTKZotHSWqj14xukiQtv17ZbPEoSY10uropSVIHnbLZ4lGSmih6DcqXJGnZdctmJ8yRJEmSJE1kz6MkdVGjKcElSdKCaJbNFo+S1EinBxFLktRBp2y2eJSkJopeg/IlSVp23bLZMY+SJEmSpInseZSkNno9S0qSpOXXK5stHiWpkU6D8iVJ6qBTNk9922qSLUm+mOT3h/Wzk3w+yT1JPpzkxNk1U5I0jaqsedHyMZslafF1yuajGfP4NuDusfV3A++pqnOAh4ArNrJhkqSjU9UroDQVs1mSFli3bJ6qeExyJvB3gN8a1gO8ErhhOOQ64NJZNFCSJD2V2SxJmrdpex7/NfDPgYPD+jOBh6vqiWF9D3DG4T6Y5Moku5LsOvDII+tqrCTpyA5W1rxMkuSaJPuS3LHK/jcl+dKw/GGS88b2fT3Jl5PcnmTXBp7ysWxDsvnxenT2LZWkY1inbJ5YPCZ5HbCvqm4b33yYQw87FLSqrq6q7VW1fcu2bdO0SZK0RqPbY9a2TOFa4OIj7P9T4GVV9SLgl4CrV+x/RVWdX1Xb13Ju+p6NzOYTs3UmbZQkjXTK5mlmW70QeH2S1wJbge9ndLXz5CTHD1c4zwT2TvODkqTZmeX4iKr6bJLnHWH/H46tfo5RNmg2zGZJWhKdsnliz2NVvb2qzqyq5wGXAZ+qqjcBnwbeMBy2A/j4ehoiSVqfYu0D8odgO+3QrYzDcuU6mnMF8AdPah58Mslt6/xeYTZL0rLols3rec7jvwCuT/LLwBeBD6zjuyRJm+/BjbilNMkrGAXU3x7bfGFV7U3ybODmJF+pqs+u97f0FGazJPWyUNl8VMVjVX0G+Mzw/l7ggqP5vCRptjb7OcRJXsRo9s+fqqpvHtpeVXuH131JPsYoPyweN4DZLEmLrVM2H81zHiVJi2yTnyWV5LnA7wI/XVV/MrZ9W5JnHHoPXAQcdlY4SZJaaZbN67ltVZK0aGZ4eTPJh4CXMxp/sQd4J3ACQFW9H3gHo8dFvG/0yEGeGG61eQ7wsWHb8cAHq+oTs2upJEkLpFE2WzxKkqZSVZdP2P8zwM8cZvu9wHlP/YQkSVqPeWezxaMkNTLL6cAlSdLR65TNFo+S1MiUDxSWJElz0imbLR4lqYmi19VNSZKWXbdstniUpC4KaBRQkiQtvWbZ7KM6JEmSJEkT2fMoSY10GlchSVIHnbLZ4lGSOmkUUJIktdAomy0eJamNtBqUL0nS8uuVzRaPktRJo6ubkiS10CibnTBHkiRJkjSRPY+S1EX1epaUJElLr1k2WzxKUieNbo2RJKmFRtls8ShJrfS5uilJUg99stkxj5IkSZKkiex5lKROGt0aI0lSC42y2eJRkjppFFCSJLXQKJstHiWpiwIazegmSdLSa5bNFo+S1Eg1uropSVIHnbLZCXMkSZIkSRPZ8yhJnTS6uilJUguNstniUZI6aTSuQpKkFhpls8WjJDWSRlc3JUnqoFM2WzxKUhdFq1tjJElaes2y2QlzJEmSJEkT2fMoSW2k1bgKSZKWX69stniUpE4a3RojSVILjbLZ4lGSOmkUUJIktdAomx3zKEmSJEmayJ5HSeqk0dVNSZJaaJTNFo+S1EXRalC+JElLr1k2WzxKUiOdHkQsSVIHnbLZMY+S1EmtY5kgyTVJ9iW5Y5X9SfJvkuxO8qUkLxnbtyPJPcOyYz2nKEnSUmmUzRaPkqRpXQtcfIT9PwWcMyxXAr8BkORU4J3AjwEXAO9McspMWypJ0rHhWuaYzRaPkqSpVNVngf1HOOQS4Ldr5HPAyUlOB14D3FxV+6vqIeBmjhx0kiRpCvPOZsc8SlIjmzyu4gzgvrH1PcO21bZLktRep2yea/H4tAce5a//2lfm+ZOStJD2fevR2Xzx+mZ0Oy3JrrH1q6vq6qP4/OF+vI6wXQvgseds5Wv/9Ec2uxmStOkee9/W2Xxxo2y251GSuphycP0RPFhV29fx+T3AWWPrZwJ7h+0vX7H9M+v4HUmSlkOzbHbMoyRpo+wE/uEws9tLgW9V1f3ATcBFSU4ZBuNfNGyTJEmztaHZbM+jJHUyw5tBk3yI0VXK05LsYTRL2wkAVfV+4EbgtcBu4C+BNw/79if5JeDW4aveVVVHGtwvSVIfjbJ5YvGYZCvwWeBpw/E3VNU7k5wNXA+cCnwB+Omqenz6U5UkbbRZDsqvqssn7C/gLavsuwa4ZhbtOhaZzZK0PDpl8zS3rT4GvLKqzgPOBy4eujzfDbynqs4BHgKuOJofliTNwAwfRKyFYjZL0rJolM0Ti8fhmSB/MayeMCwFvBK4Ydh+HXDpTFooSZKexGyWJG2GqSbMSbIlye3APkYPkPwa8HBVPTEcsupzQZJcmWRXkl2P14ymppckjTS6uqkj26hsPvDII/NpsCQdqxpl81TFY1UdqKrzGU3hegFw7uEOW+WzV1fV9qrafmJm9OwUSRKp9S1aLhuVzVu2bZtlMyXpmNYtm49qttWqejjJZ4CXAicnOX64wnnoeSGSpM20vgcRawmZzZK04Bpl88SexyTPSnLy8P77gJ8E7gY+DbxhOGwH8PFZNVKSNKVGt8ZodWazJC2RRtk8Tc/j6cB1SbYwKjY/UlW/n+Qu4Pokvwx8EfjADNspSZK+x2yWJM3dxOKxqr4EvPgw2+9lNMZCkrQgFnF8hDae2SxJy6NTNh/VmEdJ0oJrFFCSJLXQKJstHiWpiwWdmU2SpGNWs2ye6lEdkiRJkqRjmz2PktRJo6ubkiS10CibLR4lqZNGASVJUguNstniUZIa6TSuQpKkDjpls2MeJUmSJEkTWTxKkiRJkibytlVJ6qTRrTGSJLXQKJstHiWpi2bPkpIkaek1y2aLR0nqpFFASZLUQqNstniUpE4aBZQkSS00ymYnzJEkSZIkTWTPoyQ1EXqNq5Akadl1y2aLR0nqpFFASZLUQqNstniUpC6azegmSdLSa5bNjnmUJEmSJE1kz6MkddLo6qYkSS00ymaLR0nqpFFASZLUQqNstniUpEY6jauQJKmDTtls8ShJnTQKKEmSWmiUzU6YI0mSJEmayJ5HSeqiaHV1U5Kkpdcsmy0eJamRTuMqJEnqoFM2e9uqJHVS61imkOTiJF9NsjvJVYfZ/54ktw/LnyR5eGzfgbF9O9dzmpIkLY1G2WzPoyQ1Msurm0m2AO8FXg3sAW5NsrOq7jp0TFX9/NjxPwu8eOwrvltV58+uhZIkLZ5O2WzPoyRpWhcAu6vq3qp6HLgeuOQIx18OfGguLZMk6dg012y2eJSkTtZ3a8xpSXaNLVeu+PYzgPvG1vcM254iyQ8CZwOfGtu8dfjezyW5dD2nKUnS0miUzd62KkldrH9GtweravsR9meVXz2cy4AbqurA2LbnVtXeJM8HPpXky1X1tbU2VpKkhdcsm+15lKQmss5lCnuAs8bWzwT2rnLsZay4Laaq9g6v9wKf4cljLiRJaqdbNls8SpKmdStwTpKzk5zIKISeMjNbkh8GTgH+aGzbKUmeNrw/DbgQuGvlZyVJ0lGZazZ726okdTLDGd2q6okkbwVuArYA11TVnUneBeyqqkNhdTlwfVWNt+Zc4DeTHGR04fJXx2eCkySprUbZbPEoSY3M+kHEVXUjcOOKbe9Ysf6/HeZzfwj8rZk2TpKkBdQpmy0eJamTGQeUJEk6So2y2eJRkjppFFCSJLXQKJudMEeSJEmSNJE9j5LURc1+XIUkSToKzbLZ4lGSOmkUUJIktdAomy0eJamRTlc3JUnqoFM2WzxKUieNAkqSpBYaZbMT5kiSJEmSJrLnUZIa6XRrjCRJHXTK5ok9j0nOSvLpJHcnuTPJ24btpya5Ock9w+sps2+uJGlVtc5FS8NslqQl0Sybp7lt9QngF6rqXOClwFuSvAC4Crilqs4BbhnWJUmbqVFA6YjMZklaFo2yeWLxWFX3V9UXhvffAe4GzgAuAa4bDrsOuHRWjZQkSd9jNkuSNsNRjXlM8jzgxcDngedU1f0wCrEkz17lM1cCVwJsPe7p62mrJOkIQq9xFZrOerP5+JO8s1WSZqVbNk8922qSpwMfBX6uqr497eeq6uqq2l5V20/M1rW0UZI0rUa3xmiyjcjmLdu2za6BkqRW2TxVz2OSExiF0+9U1e8Omx9IcvpwZfN0YN+sGilJmk5qAZNGM2E2S9Jy6JTN08y2GuADwN1V9Wtju3YCO4b3O4CPb3zzJElTazajm1ZnNkvSkmiWzdP0PF4I/DTw5SS3D9t+EfhV4CNJrgC+AbxxNk2UJEkrmM2SpLmbWDxW1X9mNNbzcF61sc2RJK1Hp0H5Wp3ZLEnLo1M2H9Vsq5KkBdcooCRJaqFRNls8SlIjna5uSpLUQadstniUpE4aBZQkSS00yuapn/MoSZIkSTp22fMoSV1Ur1tjJElaes2y2eJRkjppFFCSJLXQKJstHiWpidDr6qYkScuuWzY75lGSJEmSNJE9j5LUSTW6vClJUgeNstniUZIa6XRrjCRJHXTKZotHSeqiaDUoX5Kkpdcsmy0eJamRHNzsFkiSpHGdstkJcyRJkiRJE9nzKEmdNLo1RpKkFhplsz2PktRIau3LVN+fXJzkq0l2J7nqMPv/UZI/T3L7sPzM2L4dSe4Zlh0bd9aSJC2uTtlsz6MkdVHMdDrwJFuA9wKvBvYAtybZWVV3rTj0w1X11hWfPRV4J7B9aOltw2cfmlmDJUnabM2y2Z5HSWpkxlc3LwB2V9W9VfU4cD1wyZRNew1wc1XtH0LpZuDitZyjJEnLpFM2WzxKkg45LcmuseXKFfvPAO4bW98zbFvpv0vypSQ3JDnrKD8rSZK+Z6Gy2dtWJamT9d0Z82BVbT/C/kzxi/8X8KGqeizJPwauA1455WclSeqnUTbb8yhJTYSZ3xqzBzhrbP1MYO/4AVX1zap6bFj9d8CPTvtZSZK66ZbNFo+S1EXV+pbJbgXOSXJ2khOBy4Cd4wckOX1s9fXA3cP7m4CLkpyS5BTgomGbJEl9Nctmb1uVJE2lqp5I8lZGwbIFuKaq7kzyLmBXVe0E/lmS1wNPAPuBfzR8dn+SX2IUcgDvqqr9cz8JSZIamXc2WzxKUiPTPhNqrarqRuDGFdveMfb+7cDbV/nsNcA1M22gJEkLplM2WzxKUidOQSNJ0mJplM0Wj5LUyKyvbkqSpKPTKZstHiWpiwIONkooSZKWXbNsdrZVSZIkSdJE9jxKUid9Lm5KktRDo2y2eJSkRjqNq5AkqYNO2WzxKEmdTPdAYUmSNC+NstniUZIa6XR1U5KkDjplsxPmSJIkSZImsudRkrooWg3KlyRp6TXLZotHSWoiQBqNq5Akadl1y2aLR0nq5OBmN0CSJD1Jo2x2zKMkSZIkaSJ7HiWpkU63xkiS1EGnbLZ4lKQumg3KlyRp6TXLZotHSWqjWj2IWJKk5dcrmy0eJamRTg8iliSpg07Z7IQ5kiRJkqSJ7HmUpE4a3RojSVILjbJ5Ys9jkmuS7Etyx9i2U5PcnOSe4fWU2TZTkjRRQQ6ufdHyMJslaUk0y+Zpblu9Frh4xbargFuq6hzglmFdkrTZqta+aJlci9ksScuhUTZPLB6r6rPA/hWbLwGuG95fB1y6we2SJEmrMJslSZthrWMen1NV9wNU1f1Jnr3agUmuBK4E2Hrc09f4c5KkqSzeRUrNz5qy+fiTvLtVkmaqUTbPfMKcqroauBrgpOOf1eh/OklaPFnAW1y0eMazeesZZ/lHI0kz1Cmb1/qojgeSnA4wvO7buCZJktas0bgKHTWzWZIWUaNsXmvxuBPYMbzfAXx8Y5ojSVqzAg6uY9GyM5sladE0y+ZpHtXxIeCPgB9OsifJFcCvAq9Ocg/w6mFdkiTNgdksSdoME8c8VtXlq+x61Qa3RZK0DqFajavQ6sxmSVoO3bJ55hPmSJLmqFFASZLUQqNstniUpE4aBZQkSS00ymaLR0nq4tCgfEmStBiaZfNaZ1uVJA3eBvsAAAqhSURBVEmSJB1DLB4lqZFUrXmZ6vuTi5N8NcnuJFcdZv//mOSuJF9KckuSHxzbdyDJ7cOycwNPW5KkhdUpm71tVZI6meG4iiRbgPcyegzEHuDWJDur6q6xw74IbK+qv0zyT4D/A/gfhn3frarzZ9ZASZIWUaNstudRktqoUUCtdZnsAmB3Vd1bVY8D1wOXPKkFVZ+uqr8cVj8HnLmhpyhJ0lLplc0Wj5KkQ05LsmtsuXLF/jOA+8bW9wzbVnMF8Adj61uH7/1ckks3qM2SJHW2UNnsbauS1EWx3ltjHqyq7UfYn1V+9akHJv8A2A68bGzzc6tqb5LnA59K8uWq+tramytJ0oJrls0Wj5LUyWynA98DnDW2fiawd+VBSX4S+F+Al1XVY4e2V9Xe4fXeJJ8BXgxYPEqSemuUzd62KkmNzHhGt1uBc5KcneRE4DLgSTOzJXkx8JvA66tq39j2U5I8bXh/GnAhMD6YX5Kkljplsz2PktTJDGd0q6onkrwVuAnYAlxTVXcmeRewq6p2Av8KeDrwH5MAfKOqXg+cC/xmkoOMLlz+6oqZ4CRJ6qlRNls8SpKmVlU3Ajeu2PaOsfc/ucrn/hD4W7NtnSRJx555ZrPFoyR1UcDB2V3dlCRJR6lZNls8SlIbUz8TSpIkzUWvbLZ4lKROGgWUJEktNMpmi0dJ6qRRQEmS1EKjbPZRHZIkSZKkiex5lKQumg3KlyRp6TXLZotHSWqjoA5udiMkSdJf6ZXNFo+S1EmjcRWSJLXQKJsd8yhJkiRJmsieR0nqotm4CkmSll6zbLZ4lKROGt0aI0lSC42y2eJRkjppFFCSJLXQKJstHiWpjWoVUJIkLb9e2eyEOZIkSZKkiex5lKQuCjjY51lSkiQtvWbZbPEoSZ00ujVGkqQWGmWzxaMkddIooCRJaqFRNls8SlIb1epZUpIkLb9e2eyEOZIkSZKkiex5lKQuCqr6DMqXJGnpNctmi0dJ6qTRrTGSJLXQKJstHiWpk0aD8iVJaqFRNjvmUZIkSZI0kT2PktRFVasHEUuStPSaZbPFoyR10ujWGEmSWmiUzRaPktRINbq6KUlSB52y2eJRktqoVlc3JUlafr2y2QlzJEmSJEkT2fMoSV0UrZ4lJUnS0muWzevqeUxycZKvJtmd5KqNapQkaY3q4NoXtWA2S9KCaZTNa+55TLIFeC/wamAPcGuSnVV110Y1TpI0vQKq0dVNHT2zWZIWS7dsXk/P4wXA7qq6t6oeB64HLtmYZkmSjlrVzK9uTurVSvK0JB8e9n8+yfPG9r192P7VJK/ZsPPWOLNZkhZJs2xeT/F4BnDf2PqeYdvKxl6ZZFeSXY/Xo+v4OUnSZhrr1fop4AXA5UlesOKwK4CHqupvAO8B3j189gXAZcALgYuB9w3fp4111Nl84JFH5tY4SdLGmnc2r6d4zGG2PaVPtqqurqrtVbX9xGxdx89Jkiapg7XmZQrT9GpdAlw3vL8BeFWSDNuvr6rHqupPgd3D92ljHXU2b9m2bQ7NkqRjV6dsXk/xuAc4a2z9TGDvOr5PkrRes701Zpperb86pqqeAL4FPHPKz2r9zGZJWjSNsnk9j+q4FTgnydnAnzHq8vz7R/rAtw88+OBN+//d/wucBjy4jt9eJJ7LYupyLl3OAzyXlX5wIxoy7js8dNP/XTecto6v2Jpk19j61VV19dj6NL1aqx0zVY+Y1u2os/mxvXse3P2//oLZvLg8l8XT5TzAc1nJbJ6QzWsuHqvqiSRvBW4CtgDXVNWdEz7zLIAku6pq+1p/e5F4Loupy7l0OQ/wXOahqi6e8U9M06t16Jg9SY4HTgL2T/lZrZPZPOK5LKYu59LlPMBzmYdu2byu5zxW1Y1V9UNV9der6lfW812SpIX3V71aSU5k1Ku1c8UxO4Edw/s3AJ+qqhq2XzbM+HY2cA7wx3Nq9zHFbJakY8pcs3k9t61Kko4hq/VqJXkXsKuqdgIfAP5Dkt2MrmpeNnz2ziQfAe4CngDeUlUHNuVEJElqYt7ZvFnF49WTD1kansti6nIuXc4DPJcWqupG4MYV294x9v5R4I2rfPZXAHvCFlenv2vPZTF1OZcu5wGeSwvzzOaMeiwlSZIkSVrdusY8SpIkSZKODXMvHpNcnOSrSXYnuWrev78eSa5Jsi/JHWPbTk1yc5J7htdTNrON00hyVpJPJ7k7yZ1J3jZsX8Zz2Zrkj5P81+Fc/vdh+9lJPj+cy4eHAcQLL8mWJF9M8vvD+rKex9eTfDnJ7Yeml17Gvy+AJCcnuSHJV4Z/Zn58Wc9FWo3ZvPnM5sVlNi8es3nzzLV4TLIFeC/wU8ALgMuTvGCebVina4GV0+1eBdxSVecAtwzri+4J4Beq6lzgpcBbhv8flvFcHgNeWVXnAecDFyd5KfBu4D3DuTwEXLGJbTwabwPuHltf1vMAeEVVnT82bfYy/n0B/Drwiar6EeA8Rv//LOu5SE9hNi8Ms3lxmc2Lx2zeJPPuebwA2F1V91bV48D1wCVzbsOaVdVnGc1QNO4S4Lrh/XXApXNt1BpU1f1V9YXh/XcY/QN3Bst5LlVVfzGsnjAsBbwSuGHYvhTnkuRM4O8AvzWshyU8jyNYur+vJN8P/ASjWcqoqser6mGW8FykIzCbF4DZvJjM5sVjNm+ueRePZwD3ja3vGbYts+dU1f0w+hc/8OxNbs9RSfI84MXA51nScxluJ7kd2AfcDHwNeLiqnhgOWZa/s38N/HPg4LD+TJbzPGD0HwmfTHJbkiuHbcv49/V84M+Bfz/csvRbSbaxnOcircZsXjBm80IxmxeP2byJ5l085jDbnO51kyR5OvBR4Oeq6tub3Z61qqoDVXU+cCajK+jnHu6w+bbq6CR5HbCvqm4b33yYQxf6PMZcWFUvYXQb3FuS/MRmN2iNjgdeAvxGVb0YeARvg1E/y/zvmnbM5sVhNi8ss3kTzbt43AOcNbZ+JrB3zm3YaA8kOR1geN23ye2ZSpITGIXT71TV7w6bl/JcDhluWfgMo7EiJyc59BzTZfg7uxB4fZKvM7pl7JWMrnYu23kAUFV7h9d9wMcY/YfDMv597QH2VNXnh/UbGAXWMp6LtBqzeUGYzQvHbF5MZvMmmnfxeCtwzjBL1YnAZcDOObdho+0EdgzvdwAf38S2TGW4X/8DwN1V9Wtju5bxXJ6V5OTh/fcBP8lonMingTcMhy38uVTV26vqzKp6HqN/Lj5VVW9iyc4DIMm2JM849B64CLiDJfz7qqr/BtyX5IeHTa8C7mIJz0U6ArN5AZjNi8dsXkxm8+ZK1Xx72pO8ltFVmy3ANVX1K3NtwDok+RDwcuA04AHgncDvAR8Bngt8A3hjVa0cuL9Qkvxt4D8BX+Z79/D/IqOxFct2Li9iNCh6C6OLIR+pqncleT6jq4SnAl8E/kFVPbZ5LZ1ekpcD/1NVvW4Zz2No88eG1eOBD1bVryR5Jkv29wWQ5HxGEyWcCNwLvJnhb40lOxdpNWbz5jObF5vZvFjM5s0z9+JRkiRJkrR85n3bqiRJkiRpCVk8SpIkSZImsniUJEmSJE1k8ShJkiRJmsjiUZIkSZI0kcWjJEmSJGkii0dJkiRJ0kQWj5IkSZKkif5/tWv+V9+zaRMAAAAASUVORK5CYII=\n",  Michael Kuron committed Oct 17, 2019 318  "text/plain": [  Martin Bauer committed Oct 21, 2019 319  "