{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Von Karman vortex street\n",
    "\n",
    "$$\n",
    "\\renewcommand{\\DdQq}[2]{{\\mathrm D}_{#1}{\\mathrm Q}_{#2}}\n",
    "\\renewcommand{\\drondt}{\\partial_t}\n",
    "\\renewcommand{\\drondx}{\\partial_x}\n",
    "\\renewcommand{\\drondy}{\\partial_y}\n",
    "\\renewcommand{\\drondtt}{\\partial_{tt}}\n",
    "\\renewcommand{\\drondxx}{\\partial_{xx}}\n",
    "\\renewcommand{\\drondyy}{\\partial_{yy}}\n",
    "\\renewcommand{\\dx}{\\Delta x}\n",
    "\\renewcommand{\\dt}{\\Delta t}\n",
    "\\renewcommand{\\grandO}{{\\mathcal O}}\n",
    "\\renewcommand{\\density}[2]{\\,f_{#1}^{#2}}\n",
    "\\renewcommand{\\fk}[1]{\\density{#1}{\\vphantom{\\star}}}\n",
    "\\renewcommand{\\fks}[1]{\\density{#1}{\\star}}\n",
    "\\renewcommand{\\moment}[2]{\\,m_{#1}^{#2}}\n",
    "\\renewcommand{\\mk}[1]{\\moment{#1}{\\vphantom{\\star}}}\n",
    "\\renewcommand{\\mke}[1]{\\moment{#1}{e}}\n",
    "\\renewcommand{\\mks}[1]{\\moment{#1}{\\star}}\n",
    "$$\n",
    "\n",
    "In this tutorial, we consider the classical $\\DdQq{2}{9}$ to simulate the Von Karman vortex street modeling by the Navier-Stokes equations.\n",
    "\n",
    "In fluid dynamics, a Von Karman vortex street is a repeating pattern of swirling vortices caused by the unsteady separation of flow of a fluid around blunt bodies. It is named after the engineer and fluid dynamicist Theodore von Karman. For the simulation, we propose to simulate the Navier-Stokes equation into a rectangular domain with a circular hole of diameter $d$. \n",
    "\n",
    "The $\\DdQq{2}{9}$ is defined by:\n",
    "\n",
    "* a space step $\\dx$ and a time step $\\dt$ related to the scheme velocity $\\lambda$ by the relation $\\lambda=\\dx/\\dt$,\n",
    "* nine velocities $\\{(0,0), (\\pm1,0), (0,\\pm1), (\\pm1, \\pm1)\\}$, identified in pylbm by the numbers $0$ to $8$,\n",
    "* nine polynomials used to build the moments\n",
    "\n",
    "$$ \n",
    "\\{1, \\lambda X, \\lambda Y, 3E-4, (9E^2-21E+8)/2, 3XE-5X, 3YE-5Y,X^2-Y^2, XY\\},\n",
    "$$\n",
    "\n",
    "where $E = X^2+Y^2$.\n",
    "\n",
    "* three conserved moments $\\rho$, $q_x$, and $q_y$,\n",
    "* nine relaxation parameters (three are $0$ corresponding to conserved moments): $\\{0,0,0,s_\\mu,s_\\mu,s_\\eta,s_\\eta,s_\\eta,s_\\eta\\}$, where $s_\\mu$ and $s_\\eta$ are in $(0,2)$,\n",
    "* equilibrium value of the non conserved moments\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\\mke{3} &= -2\\rho + 3(q_x^2+q_y^2)/(\\rho_0\\lambda^2), \\\\ \\mke{4} &= \\rho-3(q_x^2+q_y^2)/(\\rho_0\\lambda^2), \\\\ \\mke{5} &= -q_x/\\lambda, \\\\ \\mke{6} &= -q_y/\\lambda, \\\\ \\mke{7} &= (q_x^2-q_y^2)/(\\rho_0\\lambda^2), \\\\ \\mke{8} &= q_xq_y/(\\rho_0\\lambda^2),\\end{aligned}\n",
    "$$\n",
    "\n",
    "where $\\rho_0$ is a given scalar.\n",
    "\n",
    "This scheme is consistant at second order with the following equations (taken $\\rho_0=1$)\n",
    "\n",
    "$$\n",
    "\\begin{gathered} \\drondt\\rho + \\drondx q_x + \\drondy q_y = 0,\\\\ \\drondt q_x + \\drondx (q_x^2+p) + \\drondy (q_xq_y) = \\mu \\drondx (\\drondx q_x + \\drondy q_y ) + \\eta (\\drondxx+\\drondyy)q_x, \\\\ \\drondt q_y + \\drondx (q_xq_y) + \\drondy (q_y^2+p) = \\mu \\drondy (\\drondx q_x + \\drondy q_y ) + \\eta (\\drondxx+\\drondyy)q_y,\\end{gathered}\n",
    "$$\n",
    "\n",
    "with $p=\\rho\\lambda^2/3$.\n",
    "\n",
    "We write a dictionary for a simulation of the Navier-Stokes equations on $(0,1)^2$.\n",
    "\n",
    "In order to impose the boundary conditions, we use the bounce-back conditions to fix $q_x=q_y=\\rho v_0$ at south, east, and north where the velocity $v_0$ could be $v_0=\\lambda/20$. At west, we impose the simple output condition of Neumann by repeating the second to last cells into the last cells. \n",
    "\n",
    "The solution is governed by the Reynolds number $Re = \\rho_0v_0d / \\eta$, where $d$ is the diameter of the circle. Fix the relaxation parameters to have $Re=500$. The relaxation parameters related to the bulk viscosity $\\mu$ should be large enough to ensure the stability (for instance $\\mu=10^{-3}$).\n",
    "\n",
    "We compute the stationary solution of the problem obtained for large enough final time. We plot the vorticity of the solution with the function imshow of matplotlib. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Reynolds number:  5.000e+02\n",
      "Bulk viscosity :  1.000e-03\n",
      "Shear viscosity:  1.000e-05\n",
      "relaxation parameters: [0.0, 0.0, 0.0, 1.4450867052023122, 1.4450867052023122, 1.9923493783869939, 1.9923493783869939, 1.9923493783869939, 1.9923493783869939]\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAACcCAYAAABvA1AIAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJztnXuUZVV95z+/ulXd1U+6G2gEmqdIMsZERDBmgRkNOioxYmZiQiQIidGYiUbzmIBxliGuGDUziTEzPoJoREHFZ2AmyYquKD6SEQSCgkFFHh2gG7p5dNPd1Y967PnjnN+pX/3uPvuee+tW3VvV+7tWrVv37H32/p19zt3nu3+vLSEEMjIyMjKWPkYGLUBGRkZGRn+QJ/SMjIyMZYI8oWdkZGQsE+QJPSMjI2OZIE/oGRkZGcsEeULPyMjIWCbIE3pGRkbGMsFhPaGLyD+KyNsjxy8QkYdFZLTP/d0vIi803y8UkSdE5D/2s5+lBhG5QkSuGUC/QUROW8D253VdInKpiHyjQ50bReTXe+0j0e5FIrLX/E2U4/XssvwKEZl0dU5NtPcqEdkqIvtE5G9FZJMp2yQiXyjLtorIq4bh3KWIw3pCBz4KXCwi4o5fDFwbQphaqI5F5BLgfcDPhhC+2sP5rf5Ltfjo90uznxhm2RYaIYRrQwhr9Q/4r8C9wG2m2nW2Tgjh3lhbIvJjwF9T/K6OASaA95sq7wMOlWUXAR8ozxnYuUsWIYTD9g9YBewGftoc2wgcAJ5Zfj8C+BiwE9gK/HdgpCy7FPgG8D+BJ4D7gJcm+rsfeCHwOuBR4CxX/hng4VKmrwE/Zso+CnwA+HtgX9nORyke0H8A9gL/DDwF+MtSnu8BzzJtXA7cA+wB/g34eVPW+FrKdj7rjr0X+Kvy/+OAG4DHgR8CrzX1rgA+C1wDPAm8geJHNVlew7fNuH8Y2A48BPwJ0CrLPmD7B94N/BMgEVlPA75ajumjFJMQ5fiGciz3Ar8EPB94ELisvA8fL+u+DLgd2AX8C/ATpv3jgM+Vz8d9wG+Xx18Su66asWy7J8B/oHgOp8vzd0XOfUdZfqCs878X8LfyFeCP3H28puG5fwp8wnx/ajk264A15f+nm/KPA+8a5LlL9W/gAgz6D/gQcJX5/hvA7eb7x4Dry4fgZOAHwGvKskvLH+xrgRbwm8C22MRS1r+//PE/QvnCcOW/VvazkmJStnJ8lGJSOodiZTVeHnsUeHb5/cvlpPLqUp4/Ab5i2nglxQQ0QjGB7QOO7fZagJMo2M768nuLYuJ9bvn9qxQvmnHgDIrJ7ryy7Iqyn1eUcqwiMjkAf0vBrtYAm4Gbgd8oy1aX9+FS4HnlGGypGfNPAm81Y3auKQvAaeb784EpihfEylK2M4EdwE+W13lJeR9Xlm3eCrwNWAGcSsFiX2yuNTnpNbgn3+hw/o3Ar3eosyvxd3mD38hJFC+OU8yxKyiex8eB7wK/mTj/euAyd2wvxXP7LGC/K/t94P8M8tyl+ne4q1wArgZeKSKryu+vLo+pWuOXgLeEEPaEEO4H/pxiCafYGkL4UAhhujzvWIolXB1eBHwTuMMXhBA+UvZzkOIH80wROcJUuT6E8M8hhJkQwoHy2BdCCLeW378AHAghfKyU5zqKB1fb/0wIYVt5/nXA3cBzur2WEMJWiqX3K8pDPwNMhBC+KSInAOdS/JAOhBBuB65yY/b/Qgh/W8qx37cvIscALwXeHELYF0LYAbwHuLDsfwL4FeAvKJj+G0MID/p2SkxSTEjHlfIkddLADAUTPVjK9lrgr0MIN4UQpkMIVwMHgecCZwNHhxDeHkI4FAqVw4dUziZocE/mjRDChsTfuxo08Wrg6yGE+8yxT1OsIo6mGKO3icgv15y/lmLyt9hNQV5SZYM8d0nisJ/Qyx/4TuCC0qhzNvCJsvgoCua11ZyyFTjefH/YtDVR/rs20eXrgdOBq6zuXkRaIvIuEblHRJ6kYIEqg+KBSHuPmP/3R75XsojIq0XkdhHZJSK7gGe49ru5lk8A+gN+FbNjdhzweAhhj6nrxyx2HRYnAWPAdiPrX1MwdZXvZgo2LBSTSx3+oKxzs4h8V0R+rUPfO83LUmX5PZWjlOUEius8CTjOlf0h6Rf6HDS4J8OAiuQoQgj/Vr6IpkMI/0KhcvuFmvP3AuvdsfUUaqZU2SDPXZI47Cf0Eh+jeGgvBr4YQtBJ8VFmGZ7iRAqdbq/YAZxHoSqwBppXARdQ6MaPoFDvQDEZKXpOjSkiJ1GwxzcAR4YQNgB3uva7wWeA54vIFuDnmZ3QtwGbRMQyHT9m/jr89wcoWPBRhkmuDyFUBisR+S0Ktcc2ikk7ihDCwyGE14YQjqNQp72/g2dLTJZ3OFa7OoTwybLsPle2LoRwfk1bc9DgnjS53x3rOE8U//eHHc49h+Ll9dkGctQ9S98FnmnaPJXi3v2g/BsVkaeZ+s8szxnkuUsSeUIv8DGKifS1GCZSqh4+DbxDRNaVP8DfpVjm94wQwjYKNcVLROQ95eF1FJPYYxQ64j+dTx8RrKH40e0EEJFfpWCDPSGEsJNCf/s3FJPaXeXxBygMh+8UkXER+QngNcC1ieYeAU4WkZGyje3AF4E/F5H1IjIiIk+V0r1TRE6nsA/8CsVL+A9E5IxYwyLyyvKlA4WxN1Dog7XfWle7Eh8CXi8iPykF1ojIz5YvrJuBJ0XkMhFZVa6yniEiZ8euK4JO9+QRYIuIrEjI1/EawlxPFP/X6Tm7BPicW3Gpa+/GckyeA/w2hc46hmuBnxOR54nIGuDtwOdL9eI+4PPA28uxPYeC2Hx8wOcuTfRbKb9U/ygmpyeAle74RooJfCcFI3sbzsvF1Z9jaHNl9wMvNN9PKdt8J4Vq43qKJd9WihVD1RaFAfRPXHtzjgG/Dtxovp8GTJnv76AwYj1KoX/+KqVBrdtrKcsvLuv8N3d8C/B/y77uAV5vyq6g3QB6JIWHzRPAbeWxIyi8WR6k0G3+K4VuepRiIr3cnP+bFDaJlREZ/4xidbC3lOV1puz1FMbcXcAvUnq5RNp4CfCtst52itXJurLsOArD68Ol/N/Uexy7rkjbqXuyAvg7La85/6co2OYTlF5GffxNjJfXfF6k7JMU5GMvhTfVb7vyvcDzzPdXAf9OYfS9HthkyjZRGMH3lXVe5doayLlL8U/KC8vIyMjIWOLIKpeMjIyMZYI8oWdkZGQsE+QJPSMjI2OZIE/oGRkZGcsEjSZ0EdkgIp8Vke+JyF0i8lNlprIvicjd5efGhRY2IyMjI6MejbxcRORqitDfq0qf2NUUEXGPhxDeJSKXAxtDCJel2jnqqKPCySef3Aexe0MIgfbEihmLidTz1rSs1za6Qeo58WWxut2c3y2G2TNtGH9fwzxeTXHbbbc9GkI4ulO9jhO6iKwHvg2cGkxlEfk+8PwQwnYROZbC//lHUm2dddZZ4ZZbbqm+xyZY7SL2YMTKUvVj52s9f16qbUWnfuvqN3mgmj50TfqdL7r9UXZTP/VS7XWcUpO9fje+x9EyxczMTEcZYhgZKRa7em32Gv0x+1n3coidn9EdFmIiH9TLodVq3RpCOKtTvSYql1Mpgmr+RkT+VUSuKqOujglFRB/l5+bYySLyOhG5RURu2blzZxeXkJGRkZHRDZok8B+lSCH6xhDCTSLyXooczo0QQrgSuBIKhm7L+rFUjTHsFLPupe351m/KsBaSifV7RdErlMmm+l+ocai7rhjDj7H32LG6tlMrTPt9PuqbhVxJLTV0+/vtZz/DhCYT+oMU4dA3ld8/SzGhPyIixxqVy46FElIRm7S7WaKmbkaq7bp6/vh8bvZCPygLOVHE0I0qLNbvQqiRurneJpN2NxN8037nW9ZNneWEhSAhS3EMO6pcQggPAw+IiOrHz6PYWeUGisQ9lJ91iXkyMjIyMhYBTfdMfCNwbenhci/wqxQvg0+LyGsoEtu8cmFEnEVTQ9F839a9sKNO8szX+NukbFjQreFzMVQ9qX5jaLKi6ZWZd4te7nWvqsdhxLCM6VJAowk9FLvOxCys5/VXnIyMjIyMXjGQXc3tGzfFepqUeT3myMhII1fEFIP0ZerK1sSgF+ujCVtKXXednHX9NpGpU1u9INWOdQfsxfjZLcPutc1u2m5ql5gPwxxmu8xiYblcx2Igh/5nZGRkLBMsOkOvY6LzZZypt3iMYeux2Pm+bXteTJY62Hb6FUzTRI/b1FukU99NPXdS9873ERvLbjxKUsf66brWzxXBMDLMlLeWR8rGESsbxus9XLDoE3o3RsVel+R1E7KdoHo1xDX5ITRxafQRiSl1TlPMV6YmL1BblvLrbvLCnZ6ebjvmkXqpdjPpp16qTbEUjIxNVGudjnVTp1+G4WEcy6WIrHLJyMjIWCYYiFFU0cm9rYmrX4o1pYxzPveGVcF0Y4ytk73uWIq51qkprEyxdrpdEjeRt1NZp2uMMWooWLlfnXgWPjIyUtWZmpoCZu/PzMxM2z3Uz1arRavVmnPMnuf768U4azEMaoZeIlSHEU1lzEw+jczQMzIyMpYJFp2he6ZUx1ztMY8Ye1d9rDI025e202q1ahlritk00UumDK5N+puZmWlry8qaYu8pA2RqZVA3vvMN9ImNhd6fEEJ1j/RzcnISgIMHDwIFK9djytD1c3p6ek5bVt5Wq8XY2BgAK1asmPM5Pj7O6OjoHFlSq7JumeAg9OtNjNrLDcvxmvqJzNAzMjIylgkGokNvGnATCxpSeEZuPVl8/ZgHjP+MMfuYl0yqLHYNWlbnXWPP88zcnt+EvafOa+Ke1qQsBdtHjAUDjI6OcuDAAQAOHToEzGXf+qkyeB349PR0VV9ZvGXsth+gYuzj4+OsXr0aoPpcuXJldV6na54PY+83W+/VVpKx/JEZekZGRsYywaIzdO9vXedpkGK1IlIx6hgD9V4Ulg3HGLmtE5MlhRjTTuljPfNUxFIWxDwzFNY+0EvQUuq8bhGTW+VTpqxsemJiovrfy6f3ZGxsrG3lZb1WtG1l+Jbxa9v+c//+/ezduxeAVatWAbB27VqgYOzj4+NzZPLP1XzGq5OtIla3qbdKv5l5Pxl+9khZfAw0l0vKHTDlthhzi9Mf7tjYWJvLWhPYun6y1yW9iNS23TTqtc7lzpbFVE11bnixF0ETl8ZO8jaBb9sGCqmqQydb/QwhVJO8Vx/pOdb9UOHdF4FqEt60aVPVthpW9+/fP+fz4MGD1Qtgz5491TGtYyd3mDWm2v56Hac6zDeoJ6aybCKj/W0tpIqmk+wZ/UdWuWRkZGQsEywqQw8hMDU1Nceg6NmdZW0xt0P97usru7Psw7PZsbGxtmMxFUzKGNuNMTFmjI0xc71+b0BUmexYxHLQ1LktdpK3F1gGWKc+Gh0dZWJiAphVi6i8K1asmHNdsU+YZeR+LMbGxuYYT2GWadt7v27dOmCWvU9OTlYyqepFVw0TExNt7pFr1qwBZg2ni80om7DbmDG3nyH/C4mU+iijd2SGnpGRkbFMsKgMXXXQMUZWx1x9Pf1eZySEWX266kFj/aT063WZGK3bopezkz2gLtWA1TvbEHYLy95jRsI6l8aUUbWT0Tl2Xb7Mr3JUN24Nn3pMmW5qJWRXa/6YvX4/TqpLn5mZaWPaukIIIVQybNiwAZjVpe/evbti63rMX9v4+HjURXUQSNmZUuy9KZow5YUYgyarjYw0MkPPyMjIWCYYSPpcyyDr3BZnZmaiembfhtfBj42NVczcY2pqqqpng0q0P/3u+40x3W506bGAqJjeOSavnlPn9RGTt4nrWyr0PwZ/n+zKQGVXljs9PT0n5N72a+0mdX2kXDgt89RjuhoQkbbQf2Xok5OTla5d5VTZNm/ezGOPPQbAvn375nza69b6gw7iid3flMurr2P/Tz0jTa4vtspbSBfKzNrTGMgGF03UHfb/mE+wRyxKUevp5Dc2NlZN5HX9xl4WMRfBmB95nTFVRNomcm8YtOdZNQEUE2ad333suq1MXk2QMqZ2o0ayLxkf8blixYpqQtVjMcNy3RI7VteOpb8vqUyMer/tWHh1zPj4eGVEtX7rMGtAbbVabdGni43Uvfe/lVi92LMW+0yVdZLFjnPseZrvZL+UXCEH8cLPKpeMjIyMZYKBui2Ojo7Wui22Wq02pm3VDt7opm3abHxa39bR87QfrauMMraC0PMtw0ixWV82NTXVdg2e2ViVgFcbHDhwoI2BWSaq9ZSN2tVAXf53i7oALlvm71PMxVDHyeZrSRm769w8Y0bkmLyxaN+UMVXl8/fAqrt0LNXFUZn66OhoVRZTjzUJ6OkVdW3H1E8xtZVvx/9vEVuhxsa0bpUUY+g26tevGPvJYgdtrB4GZIaekZGRsUyw6G6L1mBp3a68XjKEUJtvJYQwJ5CoE5RBrly5so0ZxPSh82VbMfbiw919PpIQQhXMsnv3bgCeeOKJqq7WtznDFd41UNtZu3ZtNYax/DZ1hjF7/d5mYPX6Krtek2W3PpeL7SOVR0cRyykf++7LvB7e6s09Y9Qx3LdvX7Vy8+Ok17h///5qXNU4aldu/lr6mR6gri3L0GOrljo9t22zbpVk4cfUrlB8psypqam2MbT57+ty9PTTmHo4IzP0jIyMjGWCgexYZFlqXbCG9aKw+nFFHdMdHx9v0/NqO/v376/qazKm2N6VVp8Oc93ivCyxNAaeoYyNjbUxXW3T7tqzffv2OWWaFXB8fLzNW8Tu8uND2p988snqGtevXz/nemMBSYqUC1uM4WuZrgy0/4MHD9YyQKtL71cAS8zl1K9MWq1W1Zbqxe246fOj8MzzwIED1eqom6Rv80Unt0LPnmMrBO9h1S0bjtX146wr3cnJyShr1+/ermUZe53uPqM5Gk/oItICbgEeCiG8TEROAT4FbAJuAy4OIRxKtRFCmLNBgd0M2LvcjY6OzvG11mMK/QHqA2GjQvWY/gD1c8OGDdWP2RthY6qBlFHJH4tNkLEluTfmqjx79uypfhRHH3101abKqst8/4OYnJysfKY1ylEn1l27drVtJnHEEUcAxSRcl2Y45dqo/dsfp89weOjQoTa1ine7bIqUq2rsheRf9Hr+1NRUNS7e13zfvn1tz5N/cdr//bOakq2fRlJ/L+z1+v5szEXKJz32WRdBHHNt9FkwR0dH28bO1vHG9ViGzSa/u4w4ulG5vAm4y3x/N/CeEMLTgCeA1/RTsIyMjIyM7tCIoYvIFuBngXcAvyvFK/RngFeVVa4GrgA+0Kkty+xarVati6BlH7GAE32zW0YCBTtUVqpGLH3D79mzp+3t71UZIYQ2Zh1z44u5S8ZWElpHWaxnL8pqjzjiiIolqupF21m1alVbLhTLvnQsjj/+eGCWoe/cubP63zOpjRs3tkXLNmFC1sioqyJdGagKY3p6us0IbI3Y/WZcllHGXEahGBM1MqvRWRn6oUOH2oy4/hkYGRlpc6PtJ1LqrxgzV5n8eTbIq85ldWRkpG0lkgpWShmm/SpgamqqbdWsmJycrFVX2efCy+brZdSjKUP/S+APAL0bRwK7QgjqavEgcHzsRBF5nYjcIiK37Ny5c17CZmRkZGTUoyNDF5GXATtCCLeKyPP1cKRq9NUZQrgSuBLgrLPOCitWrIjqAE1/QPHGtwZDPaZ1PKOyrNwbHi0rUJamTEYNj7GNhhXax4oVK9p0tNZg6kPLbfCSD/lX5qyZ/7Zv387jjz8OzLrFHXnkkUCR21uvT3OOaDsbNmyo+tu1axcwu+vOMcccU8nnswi2Wi02btw45/oUdpw9K7TfddWhTFcZ+szMzJyQe5hl6FaP2i+jqDV+eyOutSfo+ChD19WRdWn0KQusW23KfbBfuvJU8FDMAKplNmeNbysW4FPnOhrboSmVVsCutqEYNx07PWbTbaicfrWTWrnZ603ZVDKaqVzOAV4uIucD48B6Csa+QURGS5a+Bdi2cGJmZGRkZHRCxwk9hPAW4C0AJUP//RDCRSLyGeAXKDxdLgGub9DWHP2q1belXJZ8mL4NBlKWpax0//79FcNVRmaZ71FHHTWnfW1T9e2WISkrtl4cdvcje55eD8wyGpXj0KFDlXzKsPW8Bx54AID77ruvalNdDJXZ2CAeH1I/Pj5e9aPsTOts2LChSjqlrow6JqOjo9V5WkcRY0He7bLValWMXD+1XwufJGtqaqq6zl7yi8d0y7FyHS/ryqmrFKs712uKhfPDXFuNd7VLrTRSYxiTv64s5ZoYwuweqt7t0iYT8yw6lkguhrp+baBbzPbgc9LHmLrKbZl6J5ms3EuNqaee1X5iPoFFl1EYSH9IoVP/cH9EysjIyMjoBV0FFoUQbgRuLP+/F3hOtx1OT0/P2VHIh7TrW2t8fHxO2DXMMsmZmZmKeSnzffTRR4FCJ/7II48AhQ5Z6wPceOON3H///cCs7lp1pirTypUrK19t7U9Z/QknnFAxa2WlynhFpGLWqpdXRgizTFGZjDL1O++8sypXeb1vfAihLRe3Hae6MHsRmePna2V64okn2uSN2TO8p4Rl6Hpf9NMydG8zsIFfvbCqTiH/vl/vGz8xMTHHTx5oi3Ow//uUFDYlsNc/xzBf5hhbofo2Dx48OGc/VSub9ev2XiMpfXXMHhB7LmLeLdqXXYHDXL99HUOvl7deYrYfjzq/+4Vm6in7SZNjTdpLxcA0xaJnW/TRfH7CsEtcvSiNdrQ3UY+pOkUn3/3791cP+bvf/W4Arrzyyur8H/3RHwVmJ3sfJGJVGDrRaV/HHnssp556KgCnn346MBsEtGvXrmop7zch3rhxYyWTvoBuvvlmoFC1QPGg64bG3v1xcnJyjmHWjtfMzExbsJLK/eSTT1YvHq+OsS9FO9nWwf9wpqamqra8akrbh/bldizXRwqpSMmYq53PdW5z39RFMFp59dNnsFy9enWb62jTSaQbVUts83J/vTZgzgfI2d9WLF9KE9Spkqwh1LsL+4kd5j6/UIy7nqdjaa/XtxHrvy7oyRrw+4XUpBojAd20Za+lU/vdIOdyycjIyFgmWPRsi2NjY3PyoMSy4sHcgBlVE6iKQN3OYJaNqsrluOOO473vfS8A73znO+f0/xtnn83FBw+yZs8eJiYm+MGP/zjbTjoJRNi8eTNQMHyvSlC1yo4dOyqVzcMPPwzAmWeeCcBTn/rUiqnqqsEaLvWat27dCsA999wDFME/ULB463IGzFm61mVrtIxIx0uZ96OPPsqOHTvmXIN11dMVhTI+HcsYg4wFZPnUCjZzpGfo1nDb7wCdmMrFGjx9XzEXPYXPXKnP3Jo1a6pjvS7zU/Xrgocs89R7reM9OTlZ65oYcwlugk5GZ9+2jpfeU7sfgQ/Cs6skfVZtcJs+//56Y7LFnCnmq37phjGnDLipNAr2/FSQVROZYlj05FyDwOkh8KmpKU676SZWzMwwEgIzIjzzppvYt24dN1x0EZQTekZGRsZSxUCyLaZcuvTNbXXbqndW1rVu3brqXGWjGvb+9a9/nXvvvReAF7/4xRy/dy//6+abaU1PM2KzNYbAikOHGHvsMS6+6ir2/d3fMX3GGUxPT1csWtm4/b5tW+Fur2zcGt3OOOMMYHYFoexu79691f+qM3/wwQer86DQ03v9omU4KQbnw9R1hbBt27bKQKxsXMd11apVbUZC75qYMtxMTU21BbNYhuWDq2ySMM/Qu2FWKSNUSo86MjJSPVs6Xpa9e/uDGs31c+3atV3l4Fc0Neb6NmN6bx1vu/qo2wGrn/rZFAuNrQK8jcIa5r1LozWSembuE/fF9M72+vtpiK47luojlkYhFZxVh/kYR5e9Dv2yW2+ltX8/I3XLR2BkYoI1F10ES8SnNSMjIyOGgTB0q1+N5feGghFpPWW86jK4b9++6q2t7n+qQ//85z/PCSecAMDHfud34HnPQ5pM1Lt3w5e/zNiLXlR5zKj7onUHVKb50EMPzSk7+uij2bJlC9CeJOvAgQPV/8r2NdDHJjLSlYiyRMtwfP50C2VAqidXVr5t27ZKh67nx9hS3W49KTYxMzPT5jViPUX0XF0F6KrjwIEDbf01cQNsgpg7nvUO8isDmzxK5dN7r2kR9BlYtWpVzymAm+jOFTEvDu+5Y90tF3KPzrpgP9tvjKl7VmqfY59czroL+378c+U90iys542Xo1vE7C2p++THIuWK6c+H+H6/vco+ELfF2LLQG1DsQ+Nd9uwxVSXoQ3POOedwyimnALDnwx9mXSR6MQbZtw8+/nEef/azqwdIc6nYKNI6A9WTTz5ZGUpVJm1ny5YtfOtb3wJmJ1s930Zr6iSiE7t9QJr436qqRfvYsWNHNclrPb0Wm1+mTp0TK1PMzMxEN07wZX5iP3jwYDSHdqyPTojJ658j/RwfH29T69kXio+aVVdVHa/YVoXdyGiRcsWM+Xz7F681FnojuY57XeRrJzmbRL1ag15sYvd+6PY348fRqmdsJCm0v8BiRsSYv35qA/dujqWcA+x1e5dXK2edKszK6M9LxQl0wrJWubR27kQa+jwLMFqy/YyMjIyliEV3W/SuPKk3l8KzRBuQpMfUeLV58+bZ7dY2byaINFO5ABx5JFu2bKlcCZWtPf3pTwcK90mvGtIl5OOPP14ZPI899lgATjzxRKBwTVQDq7J3ZTHKBDdv3lwFKanKxb65feZJO07apjJ066JYp6ppErCTgmURMTZh3dhUFv1cqK3c7PLXb41mDZ++3xUrVlQMXVdH+mmzRi6EWqNuKa+wqi3/GTuvblPwWF+x81Muc7GVW4yp+3tgV3Iqu1+VWXk9U9d+7aoyBq8Ss4y97to7GTnrGLZdETVx74ytpOoMpVaL0S2WNUM/9PKXQzk5dsLMmjXse8UrFliijIyMjIXDohtFR0ZG5ujW/JtP3+Y2F4XC7j6jUOanTPQb3/hGxXDf+IY3wPHHE+6+O5rA3SKsWsXkC17AzMxMZQhT1mbZmxpclcVrv2vXrq36Pe6446prALj99tsr46S+lZWZqwH1xBNPrEL/Y/lTvEHI6n/9nqKWlftwdbVDjI6O1u7haNlaSq9Y52oX2zvSBiHFdMHzgV1ZeLuL39nVg+UGAAAXwElEQVTK/m9ZvI6Tjo/P2zIfvWYT1DEye+9jRmyvp7Z7zvZidG7i/teJqfvxVVhW6l0aY7sZefddG5SmsO3567XsOKb79tcQuwfeNhFbpdVpGeyeyd4O0Gm11MT9NIZlzdARYeoTn2BqfDy++0aJmfFxHvvgB6FP3hYZGRkZg8CiM/QQQnLPQJsr2wbB2LJNmzZVZcqetexFL3oR1113HQCXXnopAG96z3v4iT/6I8YnJmDv3oqtT61axcyqVdz1x3/M+DOeAfv2sXXr1soFUj0cVCd+8skn84xnPAOYZegaGLR27drKK0YDm77zne8ARUCSzw6pycGe9rSnAUUmRy2LwbMPy3CUoft85KtXr27b19Ey0bp9JRVNWYQ9pu15XaO1PXi3xX4yX687t14Vnr3bT7tysZ8LoTdXNA2k8izPBnLFsmBCcd0+02bTAC5fL+X54q9lZmam1qXRXkvMK8d7TSnsKtMm+rJ17arQr3BTuvCUq2DMThRj415Pbl0y657x2IrRjluvDP2wCP3fc9pp/PM113Beq8We97+f1qOPsnt0lIdf8AL2nXsujIxQn2cwIyMjY2lgoDp0mwZAYffoVGahDNTux+l1hhqoc9ZZZ/G1r30NoGLqX/7ylwE499xzOfvsswHzhv/61wEq/fWePXsqfa/q0tX7ZN26dRWLVk8aZXIPPfRQ5V2isqhOe2JiolplaG51ZeaahnfTpk1tzCLmARAbL7Uj+B2dxsbG2hiC3TvTs9IUYkzFMwsrv/c4sMExPnFWt6gLeLGy2Dzmvsyz8LGxsbbdfRYSTbxNYt9jqx0/vvb+em+RbtE0bYFFykfdjm0sWCjlv67n+yR+3pvKnhfzC1c0SVyW8sKzq1sfJJXagzXm4RdLGeCvpSkWfUKfnp6OGg30IqwxSgdI3Qe1zrp166r/dfLUoJy9e/dy4YUXArM/6muuuQaA73//+5WK5LTTTgNmgxfuuOOO6hyfw1tVL0cccURlzPQPhJ0U/QO9cePGyq3ypJNOAqiCn1RNMz093TYJ2rGpW7JNTk62TZB267q6KLuYmiH2QPmlvL22uh+sv8f2PGu8mq/bVmxi94FFPojIymQ3gKiL9muCTvJ24w7qEZtUbEZJP5noda9evXrORjILhSbXbes22ZzaZ3C0z47fzNsa3+ueq6b3MjbB1qmH7T1IRZOmDKZ1aqscWJSRkZGRsfgM3YbKWlbq366WNfl86CGE6piyUVXLrFixomJlb37zmwG45JJLgEL1opsyK9NW2K28fFY9+5b2b3/LCNVAq+xB1SybNm2qDKtqDNXzlF3HGL4fNwtrgPGyWMOYlvlc6+Pj43OYfKyPGOyY1C0V7UbBqbDxJnkuUvBt2wAsbxy1DC62EkmpEOrQhHnH3Dy7YfMi0qYismzRu8PZ3PR1z2o/UGcktGWptAepjacVnh3HMjLa9rzqsEkovZU7tkrrxu0wZTD1ZZ0C83pFZugZGRkZywQDMYrGEh15F6sQZhNSKZu27F316spwVYd+8OBBnvKUpwCzofCqZ3/Oc57D+eefD8y6+Kmeziat8vpIlengwYNtBkiFldcz9NWrV7dlzPMBLPb6Um/zGCtVeLepGNO39gHP0D1i7NIaV2ObD/t+vF2gHwE6dewwNk5Wh65j3y90oz9uWt/Lb38ver/00yY688zV2ioWgqF7uXtl6gobnl9nyI+l/Iit/GJ9NBl7P5ZNmH3MuJkKXup3uguPzNAzMjIylgkWPTlXq9Wa417n37R+F3OY1QtatqUsWvXWypzHx8crJubzix84cKAKGlKWo2XW/clbtG06Al0ZxKzfdecdPHiwYsG6ElB0E7xhEbM1eGYQ06/b1YOOXZ2rXoyh25S3fgxiO8v4BE0i7fmje/EssUh55Vg3Rh++be9rLyy2KRutOy9W1x+bmZlpY+g2X763kcRWK4uBmDurL+vkxlfHXi0r9896jE3Xffr/U8cUdc+qZeUL4a3SKxY9H7qNbtNj9jOWZ9hGjyq8K5Z+t1uj+a3rRkdHKz9whf6orcFV69u8GBDPqWIfLH8NVpUR89W2/es1x5CaZMbGxtrcD/1mu/ZarBrI53nxsOd798eUwSc2adsfolfH9Dqhx5b7/sVj74VVVVhYY1sT1P2AO9Xvxc84hNnIan0Bx/Khqztuyi2uX4i50Vp5Y9epZTHXPkVskvaoU0vG5pRUv7G6Xt6Y22Hs/IVWo3SDrHLJyMjIWCZYdKOoNWxY9YZVE+h3vzS2xlQfGWeNYMpolLXYDIdqDNXzvMql1Wq1rQgsE0zloalbbVgZ/GfMuJIKmPEYHR2tVhe6IrERa57dad01a9bUBpzElrjebSvFAGPHrepjPlGZKXZo+44ZGb0x3q5k6qIbFalleycXtm7cQf13ayz0kZ+xyEuF/R0sBFPvZlXV1HjchA37/2OqsjoX324RW6EOOzpeqYicICJfEZG7ROS7IvKm8vgmEfmSiNxdfm5ceHEzMjIyMurQhKFPAb8XQrhNRNYBt4rIl4BLgX8KIbxLRC4HLgcuSzWkbMUH7kC7PtWyy1gQgs+Gp3rvlStXVm9TZaz6fc+ePVX7ylRjzNnLZ4MYUuH5nhnEVhuKOuOO/b+JcaXValWGVk0vYF0kffoCve7x8fG2QBW/SoLZVZEanWPXG2OzdS5zNsdIL3uJxlYrTQxyVobYvpb2XkFcv9/EDS+l0/WyWZtM3RjYIC3P1EMI1YrT7qyk37sJGIuhia0gxtTrrrfT6qqpLaLTOU1WUMsVHRl6CGF7COG28v89wF3A8cAFwNVltauBvN1PRkZGxgDRlQ5dRE4GngXcBBwTQtgOxaQvIpubtGFdxKx+MKX3SnkJWEYCcz0WPIteu3ZtpTdVZhNL4uTD8W34uHerjIUle52dZZV1DN8ips+tC6iYmpqq9OOaHVIxMTFR9euzRFoGp/DeMa1WK7pLjr9ee5366eXVvsbHxyu7Rbc703eDOl26LbO7GflgLM+0U+zSrgxS997bhGJ7ZPrnynptedns/36HJbsC69aLqBdvq5QLZ0zefrHmhfLkWapobC0QkbXA54A3hxCe7OK814nILSJyi24KkZGRkZHRfzSiSCIyRjGZXxtC+Hx5+BERObZk58cCO2LnhhCuBK4EePaznx0mJyert6plaD7ZjvWGscEsWsfmSLZ1XN/AXBbt+46FGcfC1b1MHna1UbfqqJPTnmPPi53vVwjT09MV41U9ucq/fv36apyUwdld7XUM1C/bj7PdrzSWbMuzWnv93itGV1B279Vu97qcT1nKDjE6OtoWFxDzWPLPWuw++z5s0qnYLjs+ZkFXDTp+Nue59w7SewOz9ze2I1UT9JPppjyN6upk9AcdJ3QpRv7DwF0hhL8wRTcAlwDvKj+vb9DWnBzddoKNTaJ+GWsnw7q8yqOjo23GPRs96Jf5PhgntTmDjcr0LxvrxpfKKdEk418s+MK3rT/8AwcOVPLqZG1zX/ggHv3BT05OVtvi+Xtg2/abcKt6Z3Jysm3SVth8OCqLqnrWr1/ftvnvYriEWZWJV2nZZy2WP93XSRmyYzm56zZaGB0dbXORVfWKzY7pjch63w4dOtT2QrBbMsbk1LFYDDTpJ/UyzOgeTRj6OcDFwB0icnt57A8pJvJPi8hrgH8HXrkwImZkZGRkNEHHCT2E8A2g7lV7Xrcd2mCVGCwDjhmYmsCzLcu0fbZDhV1OewZnGXsqI2Idg4uV+WuKsXCr6qnLGy8ilYFXWbB1O/SBUxpsZQ1yOk5aV+vs27evYujeZW5iYqItuMu6kFpDNMy6VK5bty7qtuqxkMEwitTqqlOQS6cye5+9C6RV7/it1BQ2f4veX80aavO2+Dwv+n10dHTJMd4mbpJLHQu9Osqh/xkZGRnLBIse+g/N9GYxFh8LpY+5D8aMqBBPHeBZkw3k8P12YujeCGtZny+LufzF2vRl3ni3cuXKSl7VrapxdM+ePW3HVM8eC9DRXaDU2GYNcqqb1fO1LrRnVJyamqpWQJqnXvdOXb16ddKQvRiIMfU6e00qcKUJQ7cup/oZS56mz6qOmxqqd+/eXa0q7YoN5iZmi+0nsNTRhM0Ow3UOk+tkZugZGRkZywSLztCtu1ssuZBlSjFdsn5PuUTV7Rpu86/XpepUTxx7nl0N1LmuxXZcsddSV2ZZvF91pNzblOVNTU1V+lNl4/q5fv16du/eDcCuXbvmtB0bH7+DPMwyc901ygYYeW8iXSmMjo5WzHzz5iLeTIOeVq5c2VOI90IgtlL0z2YqWCoVtGQ//b2zKQ/8mOv90u/T09NtK0yrN+81eGi5YJjY8TBgINkWFfbhi0VQerVLzG3R/xCtW6L3ObeGJ+/2aPNl1G3cYGVKqUWa+KN7xCaXlGrKqjl85kid0B977LFqQw4937oh6lLeG+Tskt76rcPsi2FycrKadPyWfWvXrq02wz766KOBuRt7DNukk5rY7UvZG07t8+hf2KlrtH3U+agr7KbnPu99jPTErinj8EFWuWRkZGQsEwxE5aKIsYiUu1jMUOpZuFXH+DZSLNie38RtMWb4rJMzZnSLRSLWMfsYE7N96OpCmZy6Cu7evRtNt6BlyrTXrVvXFszio2+twVWZuaoEdu/e3ZYPR1cDRx11VMXM/ZZ9MUPoMLDJWEZP+92qPmLPhb93scCimNqtbvtAq57xWSl9OxbDMJYZg0Nm6BkZGRnLBANh6KkcJ032G62rr6hjKTEWbNMC6Ll1IeK2zLPoFPtP6cJjxqxYgEWdfn5kZKRicNYoCbBp06aKWateXYNTRNp3X/LyTk1NVfpxZfHqTmc3vlaDpxpCN27cWOnMU+H9w8Qm/T2LGb9jRlBFHXu39VPpBFIrOG/jiKWLyMiAzNAzMjIylg0GEljUbfa1WFlKF62IJcmqC96J6TVjrmx1gSed+u3E7K2ePeVBEyuz+3XC3MAV1WV7D5j9+/e3se/Y7jlebmWeloWrnlx19zYBWwz9ZpVNVzmp+r5ejDE38Vyp+27bStXzbNzKUvd9odBLYE9OtjVYLPqE7lUqdUbQmGom5QYYe5BSbo+pzYBTy+Bu8rU0nZD1nLolfep6Y0ZgdW+bnp6uVCR67ClPeQpQZFLUMu9+aCcVb7TTdsbHxysfdUUsi2bqGnpBr7lVUi9Ve6xuEuqH/KkXSF0/Cz15z7f9Xu9HDIf7C6Af9zqrXDIyMjKWCQaiclF0+3ZPGR47MSyIG1xjdeoMrjYvRyzAqE6d0pR51NW37DJ2Tp3aaWRkpIoitRtHQ6E60ejPujGMrTqsWiYWfGO/zwcLwUxTgV+d+ut0D5usJlPqxLrv/cSwG1B7kW+xWf2wj2Fm6BkZGRnLBEMTWNSE0fQaQl+nU687z+vAm5R1clvsdC0po13KnmBl8mHkVi/vg1NmZmbaGHas7TqXO2s8bmL/aILFNvalbByKXhlgpxXmYmDY2WQ/cDhcYzfIDD0jIyNjmWDRGXqdPrhbL5du9NMx3blv27otxlwKfTud2puPvIqUfj6GWJZGfy02cKbOtTDlBtjE5a8T+sWC54sm49vLfRskMmM9vJEZekZGRsYywcB3LOqmTtNjdWyrqQdMqizVdhOf+pQnS10fKR/qGLuMBcP4POx1ctb169HJrztW3/dTJ8dio5MtR78PSxDNoMcrY3gxULfF+aBpkFHd8Tq3R/vD7eaHMx8XzPnUj02stizlqrcQctYhFcU5TOjmepuo1Jq6LfbyzGVkeGSVS0ZGRsYywdC4LXaLXgI5UuqGpqqIXlQmneTU43Xj0dRA7MtSaqAYUgbAbssOB6bZzxXNQo7XUjDoHg7Py2IgM/SMjIyMZYKBuC32iqb63yZG0SYMP4V+M7GF0MGn0hFoeVOk6i6E/nfQrLKbZ8eW92sMFuL6+9XmQrDpXl1eM+YiM/SMjIyMZYIl5eXSlIX3iwWnEPOS6UYH3k3b9lhM3iaudrE6vYzTfGFXDbGy2P/zQa/3vNcVSb+Z9aBXKha9PheZVS8e5jWhi8hLgPcCLeCqEMK7+iJVs75rjzeZ7LtZIjdReXSaODr118lXvVdDrS+L+aM3mTR8ndjkm4oe7TWytInxutOx2Pe6Yyk55ju59kvF1S0W21CbelY61c2YH3pWuYhIC3gf8FLg6cAvi8jT+yVYRkZGRkZ3mA9Dfw7wwxDCvQAi8ingAuDf+iFYDDGW6xlnqix2LLWVnGcPsW3mUsbGlDompgpJBQjFsh36unVs3P5vP31/fgu6mZmZtvFpsmlyJ7WXltVtTj0yMlJb1mq1kiy87jwvX+p7DE1dMpuw2m5WK90y2H6tOFMyNJW/11XZMLL2paI2mo9R9HjgAfP9wfLYHIjI60TkFhG5ZefOnfPoLiMjIyMjhfkw9Ngrq+3VGkK4ErgSQER2isg+4NF59LvYOIos70Iiy7uwyPIuLBZL3pOaVJrPhP4gcIL5vgXYljohhHC0iNwSQjhrHv0uKrK8C4ss78Iiy7uwGDZ556Ny+RbwNBE5RURWABcCN/RHrIyMjIyMbtEzQw8hTInIG4B/pHBb/EgI4bt9kywjIyMjoyvMyw89hPD3wN93edqV8+lzAMjyLiyyvAuLLO/CYqjklWF0EcrIyMjI6B45l0tGRkbGMsGiTegi8hIR+b6I/FBELl+sfptCRE4Qka+IyF0i8l0ReVN5/AoReUhEbi//zh+0rAoRuV9E7ijluqU8tklEviQid5efGwctJ4CI/IgZw9tF5EkRefOwja+IfEREdojIneZYdEylwF+Vz/R3ROTMIZH3f4jI90qZviAiG8rjJ4vIfjPWHxwSeWufARF5Szm+3xeRFw+JvNcZWe8XkdvL4wMf3yribyH/KIym9wCnAiuAbwNPX4y+u5DxWODM8v91wA8oUhpcAfz+oOWrkfl+4Ch37M+Ay8v/LwfePWg5a56Hhyl8a4dqfIGfBs4E7uw0psD5wD9QxGQ8F7hpSOT9T8Bo+f+7jbwn23pDNL7RZ6D8/X0bWAmcUs4hrUHL68r/HHjbsIzvYjH0Kk1ACOEQoGkChgYhhO0hhNvK//cAdxGJfF0CuAC4uvz/auAVA5SlDucB94QQtg5aEI8QwteAx93hujG9APhYKPBNYIOIHLs4khaIyRtC+GIIYar8+k2KGJGhQM341uEC4FMhhIMhhPuAH1LMJYuGlLxS5AP4ReCTiylTCos1oTdKEzAsEJGTgWcBN5WH3lAuXz8yLCqMEgH4oojcKiKvK48dE0LYDsVLCtg8MOnqcSFzfwTDOr6KujFdCs/1r1GsIhSniMi/ishXReR5gxIqgtgzMOzj+zzgkRDC3ebYQMd3sSb0RmkChgEishb4HPDmEMKTwAeApwJnANsplljDgnNCCGdSZLz8LRH56UEL1AllENrLgc+Uh4Z5fDthqJ9rEXkrMAVcWx7aDpwYQngW8LvAJ0Rk/aDkM6h7BoZ6fIFfZi4xGfj4LtaE3nWagEFARMYoJvNrQwifBwghPBJCmA4hzAAfYpGXfCmEELaVnzuAL1DI9ogu+8vPHYOTMIqXAreFEB6B4R5fg7oxHdrnWkQuAV4GXBRKBW+punis/P9WCp306YOTskDiGRjm8R0F/jNwnR4bhvFdrAl96NMElPqwDwN3hRD+why3OtGfB+705w4CIrJGRNbp/xSGsDspxvWSstolwPWDkbAWc1jNsI6vQ92Y3gC8uvR2eS6wW1Uzg4QUG89cBrw8hDBhjh8txT4GiMipwNOAewcj5SwSz8ANwIUislJETqGQ9+bFlq8GLwS+F0J4UA8MxfguorX4fArPkXuAtw7SElwj37kUy7nvALeXf+cDHwfuKI/fABw7aFlLeU+l8AD4NvBdHVPgSOCfgLvLz02DltXIvBp4DDjCHBuq8aV42WwHJikY4mvqxpRCJfC+8pm+AzhrSOT9IYXuWZ/jD5Z1/0v5rHwbuA34uSGRt/YZAN5aju/3gZcOg7zl8Y8Cr3d1Bz6+OVI0IyMjY5kgR4pmZGRkLBPkCT0jIyNjmSBP6BkZGRnLBHlCz8jIyFgmyBN6RkZGxjJBntAzMjIylgnyhJ6RkZGxTJAn9IyMjIxlgv8PlzellU6dZncAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np\n",
    "import sympy as sp\n",
    "import pylbm\n",
    "\n",
    "X, Y, LA = sp.symbols('X, Y, LA')\n",
    "rho, qx, qy = sp.symbols('rho, qx, qy')\n",
    "\n",
    "def bc_in(f, m, x, y):\n",
    "    m[qx] = rhoo * v0\n",
    "\n",
    "def vorticity(sol):\n",
    "    ux = sol.m[qx] / sol.m[rho]\n",
    "    uy = sol.m[qy] / sol.m[rho]\n",
    "    V = np.abs(uy[2:,1:-1] - uy[0:-2,1:-1] - ux[1:-1,2:] + ux[1:-1,0:-2])/(2*sol.domain.dx)\n",
    "    return -V\n",
    "    \n",
    "# parameters\n",
    "rayon = 0.05\n",
    "Re = 500\n",
    "dx = 1./64   # spatial step\n",
    "la = 1.      # velocity of the scheme\n",
    "Tf = 75      # final time of the simulation\n",
    "v0 = la/20   # maximal velocity obtained in the middle of the channel\n",
    "rhoo = 1.    # mean value of the density\n",
    "mu = 1.e-3   # bulk viscosity\n",
    "eta = rhoo*v0*2*rayon/Re  # shear viscosity\n",
    "# initialization\n",
    "xmin, xmax, ymin, ymax = 0., 3., 0., 1.\n",
    "dummy = 3.0/(la*rhoo*dx)\n",
    "s_mu = 1.0/(0.5+mu*dummy)\n",
    "s_eta = 1.0/(0.5+eta*dummy)\n",
    "s_q = s_eta\n",
    "s_es = s_mu\n",
    "s  = [0.,0.,0.,s_mu,s_es,s_q,s_q,s_eta,s_eta]\n",
    "dummy = 1./(LA**2*rhoo)\n",
    "qx2 = dummy*qx**2\n",
    "qy2 = dummy*qy**2\n",
    "q2  = qx2+qy2\n",
    "qxy = dummy*qx*qy\n",
    "\n",
    "print(\"Reynolds number: {0:10.3e}\".format(Re))\n",
    "print(\"Bulk viscosity : {0:10.3e}\".format(mu))\n",
    "print(\"Shear viscosity: {0:10.3e}\".format(eta))\n",
    "print(\"relaxation parameters: {0}\".format(s))\n",
    "\n",
    "dico = {\n",
    "    'box': {'x': [xmin, xmax], \n",
    "            'y': [ymin, ymax], \n",
    "            'label': [0, 2, 0, 0]\n",
    "           },\n",
    "    'elements': [pylbm.Circle([.3, 0.5*(ymin+ymax)+dx], rayon, label=1)],\n",
    "    'space_step': dx,\n",
    "    'scheme_velocity': la,\n",
    "    'parameters': {LA: la},\n",
    "    'schemes': [\n",
    "        {\n",
    "            'velocities': list(range(9)),\n",
    "            'conserved_moments': [rho, qx, qy],\n",
    "            'polynomials': [\n",
    "                1, LA*X, LA*Y,\n",
    "                3*(X**2+Y**2)-4,\n",
    "                (9*(X**2+Y**2)**2-21*(X**2+Y**2)+8)/2,\n",
    "                3*X*(X**2+Y**2)-5*X, 3*Y*(X**2+Y**2)-5*Y,\n",
    "                X**2-Y**2, X*Y\n",
    "            ],\n",
    "            'relaxation_parameters': s,\n",
    "            'equilibrium': [\n",
    "                rho, qx, qy,\n",
    "                -2*rho + 3*q2,\n",
    "                rho-3*q2,\n",
    "                -qx/LA, -qy/LA,\n",
    "                qx2-qy2, qxy\n",
    "            ],\n",
    "        },\n",
    "    ],\n",
    "    'init': {rho:rhoo, \n",
    "             qx:0.,\n",
    "             qy:0.\n",
    "    },\n",
    "    'boundary_conditions': {\n",
    "        0: {'method': {0: pylbm.bc.BouzidiBounceBack}, 'value': bc_in},\n",
    "        1: {'method': {0: pylbm.bc.BouzidiBounceBack}},\n",
    "        2: {'method': {0: pylbm.bc.NeumannX}},\n",
    "    },\n",
    "    'generator': 'cython',\n",
    "}\n",
    "\n",
    "sol = pylbm.Simulation(dico)\n",
    "while sol.t < Tf:\n",
    "    sol.one_time_step()\n",
    "    \n",
    "viewer = pylbm.viewer.matplotlib_viewer\n",
    "fig = viewer.Fig()\n",
    "ax = fig[0]\n",
    "im = ax.image(vorticity(sol).transpose(), clim = [-3., 0])\n",
    "ax.ellipse([.3/dx, 0.5*(ymin+ymax)/dx], [rayon/dx,rayon/dx], 'r')\n",
    "ax.title = 'Von Karman vortex street at t = {0:f}'.format(sol.t)\n",
    "fig.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}