{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Shan-Chen Two-Phase Single-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", "from lbmpy.macroscopic_value_kernels import macroscopic_values_getter, macroscopic_values_setter\n", "from lbmpy.maxwellian_equilibrium import get_weights" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This is based on section 9.3.2 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": [ "N = 64\n", "omega_a = 1.\n", "g_aa = -4.7\n", "rho0 = 1.\n", "\n", "stencil = get_stencil(\"D2Q9\")\n", "weights = get_weights(stencil, c_s_sq=sp.Rational(1,3))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Data structures" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "dim = len(stencil[0])\n", "\n", "dh = ps.create_data_handling((N,)*dim, periodicity=True, default_target='cpu')\n", "\n", "src = dh.add_array('src', values_per_cell=len(stencil))\n", "dst = dh.add_array_like('dst', 'src')\n", "\n", "ρ = dh.add_array('rho')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Force & combined velocity" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The force on the fluid is\n", "$\\vec{F}_A(\\vec{x})=-\\psi(\\rho_A(\\vec{x}))g_{AA}\\sum\\limits_{i=1}^{q}w_i\\psi(\\rho_A(\\vec{x}+\\vec{c}_i))\\vec{c}_i$\n", "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 = sum((psi(ρ[d]) * w_d * sp.Matrix(d)\n", " for d, w_d in zip(stencil, weights)), zero_vec) * psi(ρ.center) * -1 * g_aa" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Kernels" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [], "source": [ "collision = create_lb_update_rule(stencil=stencil,\n", " relaxation_rate=omega_a, \n", " compressible=True,\n", " force_model='guo', \n", " force=force,\n", " kernel_type='collide_only',\n", " optimization={'symbolic_field': src})\n", "\n", "stream = create_stream_pull_with_output_kernel(collision.method, src, dst, {'density': ρ})\n", "\n", "\n", "opts = {'cpu_openmp': False, \n", " 'target': dh.default_target}\n", "\n", "stream_kernel = ps.create_kernel(stream, **opts).compile()\n", "collision_kernel = ps.create_kernel(collision, **opts).compile()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Initialization" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [], "source": [ "method_without_force = create_lb_method(stencil=stencil, relaxation_rate=omega_a, compressible=True)\n", "init_assignments = macroscopic_values_setter(method_without_force, velocity=(0, 0), \n", " pdfs=src.center_vector, density=ρ.center)\n", "\n", "\n", "init_kernel = ps.create_kernel(init_assignments, ghost_layers=0).compile()" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [], "source": [ "def init():\n", " for x in range(N):\n", " for y in range(N):\n", " if (x-N/2)**2 + (y-N/2)**2 <= 15**2:\n", " dh.fill(ρ.name, 2.1, slice_obj=[x,y])\n", " else:\n", " dh.fill(ρ.name, 0.15, slice_obj=[x,y])\n", "\n", " dh.run_kernel(init_kernel)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Timeloop" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [], "source": [ "sync_pdfs = dh.synchronization_function([src.name])\n", "sync_ρs = dh.synchronization_function([ρ.name])\n", "\n", "def time_loop(steps):\n", " dh.all_to_gpu()\n", " for i in range(steps):\n", " sync_ρs()\n", " dh.run_kernel(collision_kernel)\n", " \n", " sync_pdfs()\n", " dh.run_kernel(stream_kernel)\n", " \n", " dh.swap(src.name, dst.name)\n", " dh.all_to_cpu()" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [], "source": [ "def plot_ρs():\n", " plt.title(\"$\\\\rho$\")\n", " plt.scalar_field(dh.gather_array(ρ.name), vmin=0, vmax=2.5)\n", " plt.colorbar()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Run the simulation\n", "### Initial state" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAA0AAAAF0CAYAAAAKF1nQAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAeoUlEQVR4nO3de7Bld1Un8O8yCaI8DJgGevKwsSajRkuE6oowTEkkMgakkkwJM2F8RCZWahwZwdGSgFVYWk4V1EyBWghMD0HCFMNjApIeK4AxImiVRDohPJKGSYwKbWLS4a0oVLrX/HFO8Nrc7nvPvvfcc07vz6dq1z17n332WXV376RXr99v/aq7AwAAMAbfsOgAAAAAdooECAAAGA0JEAAAMBoSIAAAYDQkQAAAwGhIgAAAgNGQAAEAAKMhAQIAAEZDAgQwMlX1iKraV1Wfq6r7qurnFx0TAOwUCRDA+LwryZ8neVySy5L896p63GJDAoCdIQECGJGqenaSdPcruvsr3f2HSf46yb9YbGQAsDMkQADjcnGS6x7cqapvSPItSe5dWEQAsIMkQADj8v1JPrNm/+lJ7u/uTy4oHgDYURIggJGoqtOSnJvkOVX10Kr67iSvSfLixUYGADvn1EUHAMCO+a4kf5nk45kMebsvya9397WLDAoAdlJ196JjAGAHVNWPJ/k33f2ji44FABbFEDiA8XhCkoOLDgIAFkkCBDAe35vkE4sOAgA2o6rOrqr3VdXBqrqtql64zjkXVNUXqurW6fayDa9rCBwAALBsqmp3kt3dfUtVPSLJzUku7e7b15xzQZJf7O5nb/a6KkAAAMDS6e57uvuW6esvZTKM+8ytXlcCBAAALLWq2pPkiUluWuftp1TVR6rq3dMlHk5oR9tgn3HGGb1nz56d/EpYOnd85FOLDgFg6Zz7hHMWHQIs1M0333x/d+9adByz+OEffFh/5rNHBn/+5o9+5bYk/7Dm0L7u3nfseVX18CTvSPKi7v7iMW/fkuTbuvtvq+pZSd6VyZp3x7WjCdCePXty4MCBnfxKWDoXPfY/LToEgKXzngOvWXQIsFBV9VeLjmFW93/2SG5671mDP3/a7j//h+7ee6Jzpot4vyPJm7v7nce+vzYh6u7rq+o1VXVGd99/vGsaAgcAACydqqokVyc52N2vPM45j5uel6o6P5P85jMnuu6OVoAAAICTRedIH53nFzw1yU8k+VhV3To99tIk5yRJd78uyXOS/ExVPZDk75Nc1hu0uZYAAQAAM+skRzO/JXW6+0+S1AbnvDrJq2e5rgQIAAAY5GjmWgGaC3OAAACA0VABAgAAZtbpHDnxdJulJAECAAAGmeccoHmRAAEAADPrJEckQAAAwFisYgVIEwQAAGA0VIAAAICZdaIJAgAAMB6rtwrQJofAVdXpVXVtVX2iqg5W1VOq6tFVdUNV3TH9+ah5BwsAACyHTufIFrZF2ewcoN9M8p7u/s4kT0hyMMlVSW7s7nOT3DjdBwAAxqCTI1vYFmXDBKiqHpnkB5JcnSTd/dXu/nySS5JcMz3tmiSXzitIAACA7bCZCtC3Jzmc5Heq6sNV9fqqeliSx3b3PUky/fmY9T5cVVdW1YGqOnD48OFtCxwAAFiczmQO0NBtUTaTAJ2a5ElJXtvdT0zyd5lhuFt37+vuvd29d9euXQPDBAAAlkvlyBa2RdlMAnQoyaHuvmm6f20mCdG9VbU7SaY/75tPiAAAwLLpJEd7+LYoGyZA3f03ST5dVd8xPXRhktuT7E9y+fTY5Umum0uEAAAA22Sz6wD95yRvrqqHJLkryfMzSZ7eXlVXJPlUkufOJ0QAAGAZLXIo21CbSoC6+9Yke9d568LtDQcAAFgFnZM4AQIAADjW0ZYAAQAAI7CqFaDNdIEDAAA4KagAAQAAM+tUjqxgPUUCBAAADGIOEAAAMAqrOgdIAgQAAAxQOdKrNwRu9SIGAAAYSAUIAACYWSc5uoL1FAkQAAAwiDlAAADAKHSbAwQAALDUVIAAAIBBjhoCBwAAjMFkHaDVG1AmAQIAAAZYzTlAEiAAAGBmq9oGe/UiBgAAGEgFCAAAGORIa4IAAACMQKc0QQAAAMbjqCYIAADAGKxqG+zVixgAAGAgFSAAAGBmndIEAQAAGI9VXAdIAgQAAMysOzmygk0QVi9iAACAgVSAAACAASpHYw4QAAAwAp3VHAInAQIAAAZZxXWAJEAAAMDMOpWjK9gGe/VSNgAAgIFUgAAAgEEMgQMAAEahkxzVBAEAABiHyhFtsAEAgDFY1QrQ6kUMAAAwkAoQAAAwiCFwAADAKHSXIXAAAMB4HOlvGLxtpKrOrqr3VdXBqrqtql64zjlVVb9VVXdW1Uer6kkbXVcFCAAAWEYPJPmF7r6lqh6R5OaquqG7b19zzjOTnDvdvj/Ja6c/j0sFCAAAmFknOZoavG14/e57uvuW6esvJTmY5MxjTrskyZt64oNJTq+q3Se6rgoQAAAwQG1qKNu2fFPVniRPTHLTMW+dmeTTa/YPTY/dc7xrbSoBqqq/TPKlJEeSPNDde6vq0UnelmRPkr9M8m+7+3ObuR4AALDaJusAbakL3BlVdWDN/r7u3nfsSVX18CTvSPKi7v7isW8fJ7TjmqUC9IPdff+a/auS3NjdL6+qq6b7L57hegAAwAo7srUZNfd3994TnVBVp2WS/Ly5u9+5zimHkpy9Zv+sJHef6JpbifiSJNdMX1+T5NItXAsAAOBrqqqSXJ3kYHe/8jin7U/yk9NucE9O8oXuPu7wt2TzFaBO8vtV1Un+x7Q09dgHL97d91TVY44T+JVJrkySc845Z5NfBwAALLNObXUI3EaemuQnknysqm6dHntpknOSpLtfl+T6JM9KcmeSLyd5/kYX3WwC9NTuvnua5NxQVZ/YbNTTZGlfkuzdu/eE4/EAAIDVcXSOTaW7+0+y/hyfted0kp+d5bqbSoC6++7pz/uq6neTnJ/k3qraPa3+7E5y3yxfDAAArK7u5Mh8K0BzsWHKVlUPmy48lKp6WJJ/neTjmYy3u3x62uVJrptXkAAAwPI52jV4W5TNVIAem+R3J3OQcmqS/93d76mqDyV5e1VdkeRTSZ47vzABAAC2bsMEqLvvSvKEdY5/JsmF8wgKAABYbpMmCDuzEOp2mmUdIAAAgK85cuIeBUtJAgQAAMysk4XO5Rlq9WpWAAAAA6kAAQAAA5gDBAAAjMhRc4AAAIAxWNWFUCVAAADAIIbAAbAyLn7/7ese3/+083Y4EgDYORIgAABgZpOFUA2BAwAARkITBAAAYBQshAoAALDkVIAAltDxGhScLN+t0QLAyUEXOAAAYBxaEwQAAGAkOpogAAAAI7KKFaDVG7QHAAAwkAoQwA5aZHODZbLZ34NmCQDLa1XbYEuAAACAQSRAAADAKHR0gQMAAEZkFbvAaYIAAACMhgoQwDbQ3GA+Zvm9apgAsMPaHCAAAGAkdIEDAABGZRUTIHOAAACA0VABAgAAZqYNNsAIaHawvNa7NxojAMxXS4AAAICxWMV1gCRAAADAzHpF22BrggAAAIyGChAAADCIOUAAJxEND1afxggA86QLHAAAMCIqQAAAwCh0NEEAAABYaipAAADA7HrSCnvVSIAAouHBmGiMALB9LIQKAACMQmc1myCYAwQAAIyGChAAADCAdYAAAIARWcUmCJseAldVp1TVh6vq96b7j6+qm6rqjqp6W1U9ZH5hAgAAy6a7Bm+LMsscoBcmObhm/xVJXtXd5yb5XJIrtjMwAABgeXWfxAlQVZ2V5EeSvH66X0menuTa6SnXJLl0HgECAABsl83OAfqNJL+U5BHT/W9N8vnufmC6fyjJmet9sKquTHJlkpxzzjnDIwUAAJbKKjZB2LACVFXPTnJfd9+89vA6p647Baq793X33u7eu2vXroFhAgAAy2YyDG7YtiibqQA9NcnFVfWsJA9N8shMKkKnV9Wp0yrQWUnunl+YAADAslnFhVA3TIC6+yVJXpIkVXVBkl/s7h+rqv+T5DlJ3prk8iTXzTFOgG1z8ftvX3QILJn1/kzsf9p5C4gEYHV0FtvMYKhZusAd68VJ/ktV3ZnJnKCrtyckAACA+ZhpIdTu/qMkfzR9fVeS87c/JAAAYBWs4DqosyVAAAAASZI+SecAAQAArGsFS0BbmQMEAAAwF1X1hqq6r6o+fpz3L6iqL1TVrdPtZZu5rgoQAAAwyJyHwL0xyauTvOkE5/xxdz97lotKgAAAgEHmuaBpd3+gqvZs93UNgQMAAGbWmVSAhm5JzqiqA2u2KweE8ZSq+khVvbuqvnszH1ABAgAAZtdJtjYE7v7u3ruFz9+S5Nu6+2+r6llJ3pXk3I0+JAECTmoXv//2RYfAilrvz87+p523gEgAWE93f3HN6+ur6jVVdUZ333+iz0mAAACAQeY5B2gjVfW4JPd2d1fV+ZlM7/nMRp+TAAEAAMPMMQGqqrckuSCTuUKHkvxKktOSpLtfl+Q5SX6mqh5I8vdJLuveOCWTAAEAAAPUXNtgd/fzNnj/1Zm0yZ6JBAgAABhmgUPghtIGGwAAGA0VIAAAYHaduQ6BmxcJEAAAMMwKDoGTAAEAAAOtXgXIHCAAAGA0VIAAAIBhDIEDAABGQwIEAACMQifRBQ4AABiLXsEKkCYIAADAaKgAAQAAw6xgBUgCBAAADGMOEAAAMBalAgQAAIxCZyWHwGmCAAAAjIYKEAAAMECZAwQAAIzICg6BkwABAADDrGACZA4QAAAwGipAAADAMCtYAZIAASeNi99/+6JD4CS33p+x/U87bwGRACyBjiYIAADAeFgIFQAAGI8VTIA0QQAAAEZDAgQAAIyGIXAAAMAg5gABLNB63bh0hmM76fgGcAxd4AAAgFHoaIIAAACwzFSAAACAYU7GClBVPbSq/qyqPlJVt1XVr06PP76qbqqqO6rqbVX1kPmHCwAALIvq4duibGYI3FeSPL27n5Dk+5JcVFVPTvKKJK/q7nOTfC7JFfMLEwAAWDq9hW1BNkyAeuJvp7unTbdO8vQk106PX5Pk0rlECAAAsE021QShqk6pqluT3JfkhiR/nuTz3f3A9JRDSc48zmevrKoDVXXg8OHD2xEzAACwDE7GClCSdPeR7v6+JGclOT/Jd6132nE+u6+793b33l27dg2PFAAAWBpbmf+zyDlAM3WB6+7PV9UfJXlyktOr6tRpFeisJHfPIT4AAGBZreBCqJvpArerqk6fvv6mJD+U5GCS9yV5zvS0y5NcN68gAQCAJbSCQ+A2UwHaneSaqjolk4Tp7d39e1V1e5K3VtWvJ/lwkqvnGCcAAMCWbZgAdfdHkzxxneN3ZTIfCAAAGKFFzuUZaqY5QAAAAF8jAQIAAEZhwd3chtpUG2wAAICTgQoQAAAwzApWgCRAAADAMBIgAABgLMwBAgAAWGISIAAAYDQMgQMAAIZZwSFwEiAAAGB2K7oOkAQIAAAYRgIEAACMhgQIYLnsf9p5X3fs4vffvoBIWDXr/dkBYPVJgAAAgJlVzAECAADGRAIEAACMwop2gbMQKgAAsHSq6g1VdV9Vffw471dV/VZV3VlVH62qJ23muhIgAABgmN7CtrE3JrnoBO8/M8m50+3KJK/dzEUlQAAAwDBzTIC6+wNJPnuCUy5J8qae+GCS06tq90bXNQcIAAAYZItzgM6oqgNr9vd1974ZPn9mkk+v2T80PXbPiT4kAQIAAIbZWgJ0f3fv3cLna51jG0ZkCBwAALCKDiU5e83+WUnu3uhDKkDA6Ox/2nlfd+zi99++gEhYFuv9mQBgA5tvZjAv+5O8oKremuT7k3yhu084/C2RAAEAAAPNcx2gqnpLkgsymSt0KMmvJDktSbr7dUmuT/KsJHcm+XKS52/muhIgAABgmDkmQN39vA3e7yQ/O+t1JUAAAMAg86wAzYsmCAAAwGioAAEAAMOsYAVIAgQAAMxu8V3gBpEAAQAAM6usvxLpsjMHCAAAGA0VIAAAYBhD4ABW0/6nnfd1xy5+/+0LiIR5W+9eAzDMKrbBlgABAADDSIAAAIDRWMEESBMEAABgNFSAAACA2bU5QAAAwJhIgABOHjrDrT4d3wDmSwUIAAAYjxVMgDRBAAAARkMFCAAAGGQVh8BtWAGqqrOr6n1VdbCqbquqF06PP7qqbqiqO6Y/HzX/cAEAgKXQW9wWZDMVoAeS/EJ331JVj0hyc1XdkOSnktzY3S+vqquSXJXkxfMLFWDxjjepXnOExdPwAGABTsYKUHff0923TF9/KcnBJGcmuSTJNdPTrkly6byCBAAA2A4zzQGqqj1JnpjkpiSP7e57kkmSVFWPOc5nrkxyZZKcc845W4kVAABYEpWTdA7Qg6rq4UnekeRF3f3FzX6uu/d1997u3rtr164hMQIAAMvoJJ0DlKo6LZPk583d/c7p4Xurave0+rM7yX3zChIAAFg+1atXAtowAaqqSnJ1koPd/co1b+1PcnmSl09/XjeXCAFWwGYn4GuWMBuNDQCW2IIrOUNtpgL01CQ/keRjVXXr9NhLM0l83l5VVyT5VJLnzidEAACA7bFhAtTdf5LJHKf1XLi94QAAAKtiFZsgzNQFDgAA4GskQAAAwFioAAFwQpolTGhuAHCSWMEEaNPrAAEAAKw6FSAAAGB2bQgcAAAwJhIgAABgDCoqQABsk51oEnC8RgsaFABwMpMAAQAAw/TqlYAkQAAAwCCGwAEAAOPQ0QQBAAAYjzq66AhmJwECGCnNDgAYIwkQAAAwjCFwAADAWGiCAAAAjENHG2wAAGA8VrEC9A2LDgAAAGCnqAABAADDrGAFSAIEAADMrLKaQ+AkQAAAwOy6V7IJgjlAAADAaKgAAQAAgxgCBwAAjIcECAAAGAsVIAAAYBw6ydHVy4A0QQAAAEZDBQgAABhm9QpAEiAAAGAYc4AAAIDxsBAqAAAwFtXDt01dv+qiqvpkVd1ZVVet8/5PVdXhqrp1uv30RtdUAQIAAJZOVZ2S5LeTPCPJoSQfqqr93X37Mae+rbtfsNnrqgABAACz6y1uGzs/yZ3dfVd3fzXJW5NcstWwJUAAAMDMKkl1D9424cwkn16zf2h67Fg/WlUfraprq+rsjS4qAQIAAIY5uoUtOaOqDqzZrjzm6rXONx6bOf3fJHu6+3uT/EGSazYK2RwgAABgEe7v7r0neP9QkrUVnbOS3L32hO7+zJrd/5nkFRt9qQoQAAAwyJyHwH0oyblV9fiqekiSy5Ls/yffX7V7ze7FSQ5udFEVIAAAYHabb2Yw7PLdD1TVC5K8N8kpSd7Q3bdV1a8lOdDd+5P8XFVdnOSBJJ9N8lMbXVcCBAAADNBzXwi1u69Pcv0xx1625vVLkrxklmtKgAAAgEE2u6DpMjEHCAAAGA0VIAAAYJg5D4Gbhw0rQFX1hqq6r6o+vubYo6vqhqq6Y/rzUfMNEwAAWCqd1NHh26JsZgjcG5NcdMyxq5Lc2N3nJrlxug8AAIxJ9/BtQTZMgLr7A5m0lFvrkvzjKqvXJLl0m+MCAADYdkObIDy2u+9JkunPxxzvxKq6sqoOVNWBw4cPD/w6AABg6fQWtgWZexe47t7X3Xu7e++uXbvm/XUAAMAOqe7B26IMTYDurardSTL9ed/2hQQAAKyEk3EO0HHsT3L59PXlSa7bnnAAAICV0EmObmFbkM20wX5Lkj9N8h1Vdaiqrkjy8iTPqKo7kjxjug8AALDUNlwItbufd5y3LtzmWAAAgBVRWexcnqE2TIAAAADWJQECAABGQwIEAACMwoNNEFbM3NcBAgAAWBYqQAAAwCCaIAAAAOMhAQIAAMahVzIBMgcIAAAYDRUgAABgdp2VrABJgAAAgGFWsA22BAgAABhEFzgAAGA8VjAB0gQBAAAYDRUgAABgdp3k6OpVgCRAAADAAKu5DpAECAAAGEYCBAAAjMYKJkCaIAAAAKOhAgQAAMxOEwQAAGA8Oumjiw5iZhIgAABgGHOAAAAAlpcKEAAAMDtzgAAAgFFZwSFwEiAAAGAYCRAAADAOvZIJkCYIAADAaKgAAQAAs+skR60DBAAAjMUKDoGTAAEAAMNIgAAAgHHolVwHSBMEAABgNFSAAACA2XXSrQkCAAAwFis4BE4CBAAADLOCTRDMAQIAAEZDBQgAAJhdt4VQAQCAEVnBIXASIAAAYJBWAQIAAMahV7ICpAkCAAAwGipAAADA7DoruQ7QlipAVXVRVX2yqu6sqqu2KygAAGAF9NHh24IMrgBV1SlJfjvJM5IcSvKhqtrf3bdvV3AAAMBy6iQ9sgrQ+Unu7O67uvurSd6a5JLtCQsAAFhq3XOvAG004qyqvrGq3jZ9/6aq2rPRNbeSAJ2Z5NNr9g9Njx0b1JVVdaCqDhw+fHgLXwcAAIzFmhFnz0xyXpLnVdV5x5x2RZLPdfc/T/KqJK/Y6LpbSYBqnWNfVwPr7n3dvbe79+7atWsLXwcAACyTPtqDt03YzIizS5JcM319bZILq2q9POVrtpIAHUpy9pr9s5LcvYXrAQAAq2S+Q+A2M+Lsa+d09wNJvpDkW0900a20wf5QknOr6vFJ/jrJZUn+/Yk+cPPNN99fVX+1he9kmDOS3L/oIPg67stycl+Wk/uynLbtvlS9djsuw4TnZTltdF++bacC2S5fyufe+wd97RlbuMRDq+rAmv193b1vzf5mRpxtalTaWoMToO5+oKpekOS9SU5J8obuvm2DzxgDtwBVdaC79y46Dv4p92U5uS/LyX1ZTu7LcnJfltPJeF+6+6I5f8VmRpw9eM6hqjo1ybck+eyJLrqlhVC7+/ok12/lGgAAAOvYzIiz/UkuT/KnSZ6T5A+7ez4VIAAAgHk53oizqvq1JAe6e3+Sq5P8r6q6M5PKz2UbXVcCNA77Nj6FBXBflpP7spzcl+Xkviwn92U5uS8DrDfirLtftub1PyR57izXrA0qRAAAACeNrbTBBgAAWCkSoJNYVf23qvpEVX20qn63qk5f895LqurOqvpkVf3wIuMco6q6aPq7v7Oqrlp0PGNVVWdX1fuq6mBV3VZVL5wef3RV3VBVd0x/PmrRsY5NVZ1SVR+uqt+b7j++qm6a3pO3VdVDFh3jGFXV6VV17fT/LQer6imel8Wqqp+f/vfr41X1lqp6qOdlMarqDVV1X1V9fM2xdZ+Pmvit6d8DPlpVT1pc5OMjATq53ZDke7r7e5P8vyQvSZKqOi+TCWLfneSiJK+pqlMWFuXITH/Xv53kmUnOS/K86T1h5z2Q5Be6+7uSPDnJz07vxVVJbuzuc5PcON1nZ70wycE1+69I8qrpPflckisWEhW/meQ93f2dSZ6QyT3yvCxIVZ2Z5OeS7O3u78lkkvhl8bwsyhsz+XvVWsd7Pp6Z5NzpdmUSC2HtIAnQSay7f3+6Im6SfDCT3ulJckmSt3b3V7r7L5LcmeT8RcQ4UucnubO77+ruryZ5ayb3hB3W3fd09y3T11/K5C9zZ2ZyP66ZnnZNkksXE+E4VdVZSX4kyeun+5Xk6UmunZ7inixAVT0yyQ9k0nEp3f3V7v58PC+LdmqSb5quf/LNSe6J52UhuvsD+fr1Z473fFyS5E098cEkp1fV7p2JFAnQePyHJO+evj4zyafXvHdoeoyd4fe/hKpqT5InJrkpyWO7+55kkiQlecziIhul30jyS0mOTve/Ncnn1/yDjmdmMb49yeEkvzMdnvj6qnpYPC8L091/neS/J/lUJonPF5LcHM/LMjne8+HvAgskAVpxVfUH03G/x26XrDnnlzMZ6vPmBw+tcyntAHeO3/+SqaqHJ3lHkhd19xcXHc+YVdWzk9zX3TevPbzOqZ6ZnXdqkicleW13PzHJ38Vwt4Wazie5JMnjk/yzJA/LZGjVsTwvy8d/1xbIOkArrrt/6ETvV9XlSZ6d5MI1q+IeSnL2mtPOSnL3fCJkHX7/S6SqTssk+Xlzd79zevjeqtrd3fdMhyTct7gIR+epSS6uqmcleWiSR2ZSETq9qk6d/qu2Z2YxDiU51N03TfevzSQB8rwszg8l+YvuPpwkVfXOJP8ynpdlcrznw98FFkgF6CRWVRcleXGSi7v7y2ve2p/ksqr6xqp6fCYT8P5sETGO1IeSnDvt0vOQTCas7l9wTKM0nVtydZKD3f3KNW/tT3L59PXlSa7b6djGqrtf0t1ndfeeTJ6NP+zuH0vyviTPmZ7mnixAd/9Nkk9X1XdMD12Y5PZ4XhbpU0meXFXfPP3v2YP3xPOyPI73fOxP8pPTbnBPTvKFB4fKMX8WQj2JVdWdSb4xyWemhz7Y3f9x+t4vZzIv6IFMhv28e/2rMA/Tf93+jUw69ryhu//rgkMapar6V0n+OMnH8o/zTV6ayTygtyc5J5O/YDy3u4+d2MqcVdUFSX6xu59dVd+eScOQRyf5cJIf7+6vLDK+Maqq78ukOcVDktyV5PmZ/GOq52VBqupXk/y7TP5//uEkP53JXBLPyw6rqrckuSDJGUnuTfIrSd6VdZ6PacL66ky6xn05yfO7+8Ai4h4jCRAAADAahsABAACjIQECAABGQwIEAACMhgQIAAAYDQkQAAAwGhIgAABgNCRAAADAaEiAAACA0fj/0wN5PNSa6lIAAAAASUVORK5CYII=\n", "text/plain": [ "