{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Transport equation with source term\n", "\n", "$$\n", "\\renewcommand{\\DdQq}[2]{{\\mathrm D}_{#1}{\\mathrm Q}_{#2}}\n", "\\renewcommand{\\drondt}{\\partial_t}\n", "\\renewcommand{\\drondx}{\\partial_x}\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 propose to add a source term in the advection equation. The problem reads\n", "$$\\drondt u + c \\drondx u = S(t, x, u), \\quad t>0, , \\quad x\\in(0, 1),$$\n", "\n", "where $c$ is a constant scalar (typically $c=1$).\n", "Additional boundary and initial conditions will be given in the following.\n", "$S$ is the source term that can depend on the time $t$, the space $x$ and the solution $u$.\n", "\n", "In order to simulate this problem, we use the $\\DdQq{1}{2}$ scheme and we add an additional `key:value` in the dictionary for the source term. We deal with two examples." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## A friction term\n", "\n", "In this example, we takes $S(t, x, u) = -\\alpha u$ where $\\alpha$ is a positive constant. \n", "The dictionary of the simulation then reads:" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "%matplotlib inline\n", "import sympy as sp\n", "import numpy as np\n", "import pylbm\n", "import math" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYQAAAD8CAYAAAB3u9PLAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAvtUlEQVR4nO3de3RU9b3//+d7JpMLSbiFcA13EQSEgCGK9QJqrXgB29pW5bS2VZHVY89lfc9Z0tavv+/XSy9n+fut1rZnWazWb0+1ytcq2lMrWrloAVuCQIAgELkmYAjhFgK5zvv3x56dhBDIZfbMnknej7VmZWbPnr0/s5PMaz6X/dmiqhhjjDEBvwtgjDEmMVggGGOMASwQjDHGRFggGGOMASwQjDHGRFggGGOMATwKBBG5RUR2ikipiCy5wDpzRGSziGwXkTVe7NcYY4x3JNrzEEQkCOwCPg+UARuAe1S1pNU6/YF1wC2qekBEBqvqkah2bIwxxlNe1BAKgVJV3aOq9cArwII269wLvK6qBwAsDIwxJvGkeLCNEcDBVo/LgCvbrHMpEBKR1UA28DNV/W17GxORRcAigMzMzCsmTZrkQRGNMaZ32Lhx41FVze3Oa70IBGlnWdt2qBTgCuBGIANYLyIfqequ816ouhRYClBQUKBFRUUeFNEYY3oHEdnf3dd6EQhlwMhWj/OAQ+2sc1RVa4AaEfkAmI7T92CMMSYBeNGHsAGYICJjRSQVuBt4q806bwLXikiKiPTBaVLa4cG+jTHGeCTqGoKqNorIw8AKIAi8oKrbRWRx5PlnVXWHiLwDFANh4Nequi3afRtjjPFO1MNOY6m9PoSGhgbKysqora31qVT+SE9PJy8vj1Ao5HdRjDEJTEQ2qmpBd17rRR9CXJWVlZGdnc2YMWMQaa8/u+dRVaqqqigrK2Ps2LF+F8cY00Ml3dQVtbW15OTk9JowABARcnJyel2tyBgTX0kXCECvCgNXb3zPxpj4SspAMMYY4z0LhG64+uqrO1zngQceoKTEmc7phz/8YZdfn5WV1b3CGWNMNyXdKKMdO3Zw2WWX+VSi7snKyuL06dNRvyYZ37sxJr6iGWVkNYRucL+9r169mjlz5nDXXXcxadIkFi5ciBuwc+bMoaioiCVLlnD27Fny8/NZuHDhOa8/ffo0N954IzNnzuTyyy/nzTff9OcNGWMMSTjstLVYdbR2pda0adMmtm/fzvDhw/nc5z7H2rVrueaaa5qf//GPf8wvfvELNm/efN5r09PTeeONN+jbty9Hjx7lqquuYv78+daBbIzxhdUQolRYWEheXh6BQID8/Hz27dvX6deqKt///veZNm0aN910E+Xl5VRUVMSusMYYcxFJXUNIhP6PtLS05vvBYJDGxsZOv/all16isrKSjRs3EgqFGDNmjJ1rYIzxjdUQ4iAUCtHQ0HDe8pMnTzJ48GBCoRCrVq1i//5uz1prjDFRs0CIg0WLFjFt2rTmTmXXwoULKSoqoqCggJdeegm7GJAxxk827DSJ9Ob3bozpHBt2aowxJmoWCMYYYwALBGOMMREWCMYYYwALBGOMMREWCMYYYwALBGOMMREWCMYYYwALhG771a9+xbBhw8jPz2f69Ol85StfYe/evRd9zR/+8AeuvPJKpk+fTkFBAStWrIhTaY0xpmOeBIKI3CIiO0WkVESWtPP8HBE5KSKbI7fHvNivn4qLi3n88cfZvHkzW7Zs4cYbb+RLX/rSBSfce/nll3n66ad588032bJlC7///e+57777KCsri3PJjTGmfVEHgogEgV8C84DJwD0iMrmdVT9U1fzI7fFo9+u3rVu3MnXq1ObHixcv5rPPPuPgwYPnrVtTU8OSJUtYtmwZQ4cOBWDChAnMmTOH999/P25lNsaYi/GihlAIlKrqHlWtB14BFniw3Q6JxObWGdu2bWPKlCnnLMvIyOD48ePnrfvKK68wc+ZMRo4cec7ytLQ0zpw50+333xupKjU1NX4Xw5geyYtAGAG0/lpcFlnW1mwR2SIifxaRKe08D4CILBKRIhEpqqys9KB43jt48CDZ2dn07du3eVlDQwOHDx8G4P777+euu+5qfm7btm1Mnz79vO1s2bKFSZMmsXz5ch588EEWLFjAu+++G/s3kKR27tzJrFmzGDJkCFu3bvW7OMb0OF4EQnvfqds2pH8MjFbV6cDPgeUX2piqLlXVAlUtyM3NveiOVWNz60hxcfF5tYPf/OY33HDDDUyfPp3nn3/+nOf69u1LfX39OcvWr19PTU0N119/PXfeeSfPPfccL774Iq+++mrHBeiFnn/+eWbOnMnGjRupqanhueee87tIxvQ4XgRCGdC6LSQPONR6BVU9paqnI/ffBkIiMsiDffuibf/Bu+++y49+9COefvrpdte/7bbbWLZsGW6NZ9euXTzwwAO88MILBAItv4Inn3ySf/zHf4xt4ZPQypUreeCBBzhz5gw33XQTAK+++mqXrk5njOmYF5fQ3ABMEJGxQDlwN3Bv6xVEZChQoaoqIoU4QVTlwb59sXXrVlavXs3777+PqnLZZZfxzjvvMHHixHbXLyws5NFHH+Wmm26irq6OpqYmfvvb3zJ79mzAaRdfsmQJ8+bNY+bMmfF8K0nhT3/6EwDf/e53+dnPfsakSZPYtWsXK1eu5Oabb/a5dMb0HFHXEFS1EXgYWAHsAJap6nYRWSwiiyOr3QVsE5EtwDPA3ZrIV+bpwEsvvUR5eTkbN27k448/5qWXXmoOg6qqKhYvXsymTZv40Y9+1Pyab33rW2zZsoU1a9aQmppKZmZm83M///nP+ctf/sJrr73Gs88+G/f3k+g++OADAO644w5EhHvvdb5vvPzyy34Wy5gex66YlkR643s/ffo0/fv3B+DEiRNkZWWxa9cuJk6cSHZ2NhUVFWRkZPhbSGMSiF0xzfRY69evp6mpiZkzZ5KVlQXApZdeSkFBAdXV1fz3f/+3zyU0puewQDAJzW0uuu66685ZvnDhQsBpvjPGeMMCwSS0Dz/8EIBrr732nOXueR4rV6684HQhxpiusUAwCauuro6PPvoIgGuuueac5/Ly8hg0aBDV1dU2H5QxHrFAMAmrqKiIuro6pkyZQk5OznnPuycHbt++Pd5FM6ZHskAwCetCzUUuCwRjvGWBYBKW26F8oUCYPNmZVNcCwRhvWCAkkOXLl1NSUuJ3MRKCqrJ+/Xqg4xqCHTNjvGGBkEAsEFpUVlZy4sQJ+vfvT15eXrvrtA4EG2lkTPQsELrhd7/7HYWFheTn5/PQQw/xt7/9jWnTplFbW0tNTQ1Tpkxh27ZtnD59mhtvvJGZM2dy+eWX8+abbzZv47e//S3Tpk1j+vTpfP3rX2fdunW89dZb/Pu//zv5+fl8+umnPr5D/5WWlgJwySWXIBe4SEVubi65ublUV1e3e2EiY0zXeDG5nX86ezWbrrrIt80dO3bw6quvsnbtWkKhEN/5znfYuXMn8+fP59FHH+Xs2bP8wz/8A1OnTqWxsZE33niDvn37cvToUa666irmz59PSUkJTz31FGvXrmXQoEEcO3aMgQMHMn/+fG6//fZzrqXQW+3evRtwAuFiJk+ezJo1a9i+fTujRo2KR9GM6bGSOxB88P7777Nx40ZmzZoFwNmzZxk8eDCPPfYYs2bNIj09nWeeeQZw2sG///3v88EHHxAIBCgvL6eiooKVK1dy1113MWiQMwP4wIEDfXs/icqtIUyYMOGi602ZMoU1a9ZQUlLCvHnz4lE0Y3qs5A4EH9qNVZX77rvvnJlMAT777DNOnz5NQ0MDtbW1ZGZm8tJLL1FZWcnGjRsJhUKMGTOG2tpaVPWCzSDG0dkagg09NcY71ofQRTfeeCOvvfYaR44cAeDYsWPs37+fRYsW8cQTT7Bw4UIeeeQRAE6ePMngwYMJhUKsWrWK/fv3N29j2bJlVFVVNW8DIDs7m+rqah/eVeLpSg0BLBCM8UJy1xB8MHnyZJ588kluvvlmwuEwoVCIBQsWkJKSwr333ktTUxNXX301K1euZOHChdxxxx0UFBSQn5/PpEmTAOdD7Ac/+AHXX389wWCQGTNm8OKLL3L33Xfz4IMP8swzz/Daa68xfvx4n9+tP1T1nE7li3HPRXBHGlnNy5jus+shJJHe8t4rKysZPHgwffv25cSJEx1+yA8ePJjKykr2799vHcum17PrIZgepXVzUWe+8VuzkTHesEAwCaezHcouCwRjvJGUgZDIzVyx0pvec2c7lF3uenv27IlZmYzpDZIuENLT06mqqupVH5CqSlVVFenp6X4XJS4626HsGj16NEDzKC5jTPck3SijvLw8ysrKqKys9LsocZWenn7BOX16mq42GVkgGOONpAuEUCjE2LFj/S6GiRFVbQ6EzjYZtQ4EG3pqTPd50mQkIreIyE4RKRWRJRdZb5aINImITdZj2lVVVcXJkyfJzs4mNze3U68ZMGAAWVlZnD59muPHj8e4hMb0XFEHgogEgV8C84DJwD0iMvkC6/0EWBHtPk3P1ZlZTtsSkeZawoEDB2JWNmN6Oi9qCIVAqaruUdV64BVgQTvrfRf4A3DEg32aHqqrI4xc1o9gTPS8CIQRQOvJ6Msiy5qJyAjgi8CzHW1MRBaJSJGIFPW2jmPT9Q5llwWCMdHzIhDaq9e3HRP6U+ARVW3qaGOqulRVC1S1oLNtyKbn2LdvHwDjxo3r0uvcKSssEIzpPi9GGZUBI1s9zgMOtVmnAHgl0iY8CLhVRBpVdbkH+zc9iHvls5EjR3aw5rmshmBM9LwIhA3ABBEZC5QDdwP3tl5BVZvHiYrIi8B/WxiY9pSVlQF0+ZwLCwRjohd1IKhqo4g8jDN6KAi8oKrbRWRx5PkO+w2MAeccBKshGOMfT05MU9W3gbfbLGs3CFT1m17s0/Q8VVVV1NbW0q9fP7Kzs7v02mHDhhEKhaisrOTMmTP06dMnRqU0pudKurmMTM/lNhd1tXYAEAgEml9n5yIY0z0WCCZhuM1F3Z2zyU5OMyY6FggmYXS3/8Bl/QjGRMcCwSQMCwRj/GWBYBJGd4ecuuzkNGOiY4FgEobVEIzxlwWCSRgWCMb4ywLBJIRwOBx1k5EbJOXl5TQ2NnpWNmN6CwsEkxCOHj1KfX09AwYMIDMzs1vbSEtLY9iwYTQ1NVFeXu5xCY3p+SwQTEKItrnI5dYuLBCM6ToLBJMQvAqEESOcS3FYIBjTdRYIJiFYIBjjPwsEkxCi7VB2WSAY030WCCYhWA3BGP9ZIJiEYIFgjP8sEExCsEAwxn8WCMZ34XC4+QPc/UDvrtaBoKpRl82Y3sQCwfiuoqKCxsZGBg0aREZGRlTbys7OJjs7m9raWo4fP+5RCY3pHSwQjO+8GmHksmYjY7rHAsH4zqvmIpedrWxM91ggGN8dOnQI8C4Q3O24NQ9jTOdYIBjfeV1DsCYjY7rHk0AQkVtEZKeIlIrIknaeXyAixSKyWUSKROQaL/Zregb3g3v48OGebM8CwZjuiToQRCQI/BKYB0wG7hGRyW1Wex+Yrqr5wLeBX0e7X9NzWA0hSTU1wRtvQKTJzyQ/L2oIhUCpqu5R1XrgFWBB6xVU9bS2DArPBGyAuGkWqz4EC4TYaP5PXrIEvvQlyM+Hv//dzyIZj3gRCCOAg60el0WWnUNEviginwB/wqklGANYk1Eyqa+H66+Hhwa/AU8/7SysrIQ5c2D5cj+LZjzgRSBIO8vOqwGo6huqOgm4E3jighsTWRTpZyiqrKz0oHgmkdXU1HDy5EnS0tLIycnxZJuDBw8mGAxy9OhR6urqPNmmcfzwh3D4w938R+U3Afjd5T+hduH9cPasU1v4+GN/C2ii4kUglAGtJ6DJAy7YqKiqHwDjRWTQBZ5fqqoFqlqQm5vrQfFMInObi4YPH45Ie98tui4YDDJs2LBztm+it3UrPPUU/I5/oB+neCN4F1/f+u8sPPMcfPvbTlvSiy/6XUwTBS8CYQMwQUTGikgqcDfwVusVROQSify3i8hMIBWo8mDfJsl53VzksmYjbzU2Op/54xp3ciV/hwEDmLbheVJShDffEiq/+h1nxWXLnM5mk5SiDgRVbQQeBlYAO4BlqrpdRBaLyOLIal8GtonIZpwRSV9Tm3nM4P0II5cFgreWLoWiIvhGv8h3vdtvZ/yMvtx+u/P5/5vNM2H8eKiogDVr/C2s6TZPzkNQ1bdV9VJVHa+qT0WWPauqz0bu/0RVp6hqvqrOVtW/erFfk/y8HmHkskDw1rJlzs8Hh0QCYf58AO6/33n4/AuCfu1u58Err8S5dMYrdqay8ZU1GSW+6mpYuxaGBCrJLV0HqanwhS8AcMstMGwY7NoFmyZGAuEPf4CGBh9LbLrLAsH4ypqMEt/KlU4fwsPj3kbCYZg7F7KzAUhJgfvuc9b7+aqpMGUKHDsGf/mLjyU23WWBYHxlTUaJ7513nJ9fTj23ucj17chZRcuWQe0XrdkomVkgGF/FuoZgM55GR9UJhDRquXTvCmfhHXecs86ECXDNNXDmDLw/8KvOwj/+EcLhOJfWRMsCwfgmHA6fcx6Cl9xAOHTokF1KMwq7d8O+fbAgexXBszUwYwa0c93refOcn3/aNQGGD4fjx52OBZNULBCMb6qqqmhoaGDAgAFRXzqzrczMTPr37099fT1Hjx71dNu9idtc9M1hkTttageuOXOcn6vXCMye7Tz46KPYFs54zgLB+CZWzUUu60eI3opIK1FB49+cO9dd1+56BQXQpw/s2AHVU69yFq5fH4cSGi9ZIBjfxGrIqcsCITq1tbBqFYSoZ1D5ZmfhFVe0u25qKlx9tXN/QzASCFZDSDoWCMY3sRph5LJAiM7Gjc6cdXeO34bU1cGll0L//hdc3202eqvsCmc86rZtzkkMJmlYIBjfWJNRYnMnLr19SORaB4WFF13fDYT3/prhXCMhHHbmuzBJwwLB+MaajBLbpk3Oz4LwBufOrFkXXX/WLMjIgJISODPN+hGSkQWC8Y01GSU2t4YwurJzgdC6H2FrlvUjJCMLBOMbazJKXHV1sH07ZFJDn73bnT6B/PwOX+c2G719rNXQUzsPJGlYIBjfuGcR5+XlxWT7Fgjdt22bM3/RnaM+duYvuvxypz2oA24gvL5pLOTmOpfX3Ls3toU1nrFAML44e/YsVVVVhEIhYnVlvNzcXEKhEMeOHePs2bMx2UdP5TYX3ZLTueYi16xZEArB9hKhocD6EZKNBYLxRevmokAgNn+GgUDALqXZTW4gFGhkhFEnAyEtDaZPd1qJDg670lm4YUMMSmhiwQLB+CLWzUUuazbqHneE0eiKyId5B0NOW3NX/Tic79wpLvauYCamLBCMLywQEldjI2zZAgOpIuPwHqfvYPLkTr/eDYT3KqY5d4qLrWM5SVggGF/EOxBsGuzO27nTmbZi3tDNzoL8fGeUUSe5rUvvbMuDAQOgqgoOH/a8nMZ7FgjGF/EKBHf7VkPoPLf/4IbBW507l1/epddPnOhcUO3AQaF+Uqtagkl4FgjGF24gjGxnbn0vWZNR17mBMDO0zbnTxUAIBlvmwCvPsUBIJhYIxhfWh5C43A7lMdXdqyFASz/CFiKBsGWLByUzseZJIIjILSKyU0RKRWRJO88vFJHiyG2diEz3Yr8meR08eBCwQEhE27eDEKZv2XZnwdSpXd6G24+wstJqCMkk6kAQkSDwS2AeMBm4R0TaDknYC1yvqtOAJ4Cl0e7XJK+6ujqOHDlCMBhkyJAhMd2XO3HeoUOHCNs1fjtUWQlHj8LUzH0EztTAsGGQk9Pl7bg1hNd3TUVF4JNPnPkwTELzooZQCJSq6h5VrQdeARa0XkFV16nq8cjDj4DYfi00Ca31dZSDwWBM95WRkcHAgQNpbGyksrIypvvqCUpKnJ9fGN795iJwLrs8eDCUH+9Dw5gJzljWTz7xqJQmVrwIhBHAwVaPyyLLLuR+4M8e7NckqXj1H7is2ajz3ECYnRVdIIi01BI+G2zNRsnCi0CQdpa1exaKiMzFCYRHLrgxkUUiUiQiRfaNrmeyQEhcbiBMCUcCoRv9By63H2F7wAIhWXgRCGVA67GDecB5E8eIyDTg18ACVa260MZUdamqFqhqQawmPTP+incguPtxO7LNhW2P9CMPP969IaetuTWED07YSKNk4UUgbAAmiMhYEUkF7gbear2CiIwCXge+rqq7PNinSWLxDoTRo0cDsH///rjsL5mVlEAqdWSV73TafbowZUVbBQXOzzf2RAYVWg0h4UUdCKraCDwMrAB2AMtUdbuILBaRxZHVHgNygP8Ukc0iYhda7cX8CoQDBw7EZX/JqqoKKipgRvonSFMTXHJJp66BcCGDBsG4cbCzbjRNmdnOxisqPCyx8VrnJyi5CFV9G3i7zbJnW91/AHjAi32Z5BfvQBg1ahRgNYSO7Njh/Pz8sG3OQPEomotchYWwZ49wdOjlDPl0ndMmFeOhxqb77ExlE3dWQ0hMbofylX2iG2HUmtuxvDMU6Zzeti3qbZrYsUAwcdXQ0MDhw4cRkeaL18SaexGeQ4cOUV9fH5d9JiM3EC5r9C4Q3I7l9dUWCMnAAsHE1WeffYaqMnToUEKhUFz2GQqFGD58OKpq02BfhBsIw6uiH3LqmjEDAgF477AFQjKwQDBxFa9ZTtuyZqOOlZRAX06ScfQgpKc7ncpRysx0cmVLuFUg2MVyEpYFgomrePcfuGzo6cWdPAnl5XBFauQb/OTJzjzWHpg1C46SS03WYKiuBjsfJGFZIJi4cr+hxzsQbKTRxbkjjG7s5kVxLsbtR9ibac1Gic4CwcTVnj17ABg7dmxc92tNRhfn9h/Myoh8WHvQf+ByRxptrLVASHQWCCau9u7dC/gXCFZDaJ8bCBMbvK8hTJ3qdEmsPTnFWWCBkLAsEExcuYEwbty4uO7XmowuzgkEZdhR7wMhFHJGG23DagiJzgLBxE04HPa9hnDgwAHURrmcp6QEhnGY1NPHYeBA58I4HioshO1MadlZU5On2zfesEAwcfPZZ59RV1fHoEGDyMrKiuu+s7KyGDhwYPPV2kyL06dh/36YmdLq/ANpb1b77ps1C07Rj8qMkc6V0z791NPtG29YIJi48at24LJ+hPa5I4zmDvK+ucjljjQqbrJmo0RmgWDixh1hFO/+A5fbj2Ajjc7ldigXpMcuEC65BPr3h431FgiJzALBxI3VEBKTGwgT6rwfcuoScZqNmjuWt271fB8mehYIJm4sEBJTSQkEaGLI0UgyxCAQwGk22kqk9mGBkJAsEEzc+HVSmsuajNpXUgKXUEqwoRZGjYJ+/WKyn1mzoITJNEoK7NoFNTUx2Y/pPgsEEzd+nYPgshrC+c6cgb17YXog+msod6SwEOpJ4xO5zJngzvoREo4FgomL+vp6ysrKCAQCzd/U480C4Xw7dzqfzXMHbHYWxDAQhg2DESNgUzhyjeUtW2K2L9M9FggmLvbv34+qkpeXF7frILSVm5tLVlYWx48f59ixY76UIdG4HcqFoU3OnRkzYrq/wkLYQiQQNm+O6b5M11kgmLjwu7kIQES49NJLAdi5c6dv5UgkzSOMauITCLNmwWbynQdWQ0g4FggmLvweYeSaOHEiAJ988omv5UgUJSWQyxH6Vh+C7GwYPz6m+zunhrBlC4TDMd2f6RoLBBMXfo8wcrmBYDUER0kJzCBSO5g+3bneZQzNng0nQ7mUM9wZZRT5uzCJwQLBxEWi1BAmTZoEWCAA1NZCaSnMkM3Oghg3FwH06QNXXmn9CInKk0AQkVtEZKeIlIrIknaenyQi60WkTkT+zYt9muSSCH0IYDWE1nbscFpsrsuK1BDy8+Oy37lzrR8hUUUdCCISBH4JzAMmA/eIyOQ2qx0D/gl4Otr9meSUKE1GEyZMAKC0tJTGxkZfy+K34mLnZz7x6VB2zZ1rNYRE5UUNoRAoVdU9qloPvAIsaL2Cqh5R1Q1Agwf7M0nGHeaZkZHB0KFDfS1LZmYmI0eOpKGhobnW0lsVF0Mmpxl2erdzFZspU+Ky39mzoSSUD0DTJqshJBIvAmEEcLDV47LIsm4RkUUiUiQiRZWVlVEXzvhva2TemsmTJyMez7PfHdZs5CguhulsQVSdMEhNjct+09NhyOcu4QwZBMsPgp0TkjC8CIT2/sO7fUkqVV2qqgWqWpCbmxtFsUyicAPh8hieBdsVFgiO4uJWI4zi1Fzkuv6GYMtEd9aPkDC8CIQyYGSrx3nAIQ+2a3qIRAsEG2kEFRVw5AhcGaczlNuaOxc+ZqbzYMOGuO7bXJgXgbABmCAiY0UkFbgbeMuD7ZoeItECwWoILR3KV6b6EwiFhfBx6mwAalevj+u+zYWlRLsBVW0UkYeBFUAQeEFVt4vI4sjzz4rIUKAI6AuEReRfgMmqeira/ZvEpqpsi8xqaYGQOIqLIZU6xp7d7iyYNi2u+09NhcZZs2Et6Np1zgx7CdC/1NtFHQgAqvo28HabZc+2uv8ZTlOS6WUOHjzIqVOnGDRoEEOGDPG7OADk5eWRkZFBRUUFJ06coH///n4XKe6Ki6GAIkLheqdDuW/fuJdh+pcvoXLtIHJPHXHm4Pb5HBVjZyqbGGvdXJQII4wAAoFAr5/krrgYruVD58G11/pShi9+SVjH1QDUr17nSxnMuSwQTEwlWv+Bqzd3LDc0OHMYNQfCNdf4Uo7Ro+HACKcfofw1C4REYIFgYipRA6E3z3q6axc01Ie5NrDWWeBTDQGg/zynhhD4m3UsJwILBBNTiRoIbnk+/vhjn0sSf8XFMJVt9A2fdK6h7NMV7AAKv1NAAynkHSum4Vi1b+UwDgsEEzMNDQ3N38CnxGlahM66+mrnm+n69esJ97I5+ROh/8A1cUYfPknPJ0iYrc//3deyGAsEE0M7d+6koaGBcePGkZWV5XdxzjF8+HBGjx7NqVOnKHEvG9ZLFBUlTiAAVF/uhHP5a9Zs5DcLBBMzidpc5HJrCevW9Z4OzaYm+NtHmlCBMOxLzu+hz+Z19PIJaH1ngWBiJtEDYfZsZ4RLbwqEkhIYdHovIzgEOTlw2WV+F4kx9zi/h5n163nrdUsEP1kgmJj56KOPAMiP04VXuqo31hDWr28z3DQBzg2RUSM5njuBAZxg9VNr/S5Or2aBYGLi9OnT/PWvfyUQCDB37ly/i9OuadOm0adPH3bv3k1vmWp9/XqYw2rngU/nH5xHhD733AnA2OLlzfMsmfizQDAxsXr1ahoaGpg1axYDBw70uzjtCoVCFBYWAs5oo96gaF09C3jTeTBvnr+FaSXta3cCcCfL+fkz3Z4930TJAsHExIoVKwD4whe+4HNJLq43NRtVVcHIXX9hIMcJT54StyukdcqVV9I4aAhj2ceW/yqmqsrvAvVOFggmJiwQEs9HH8HdvAJA4O6v+VyaNoJBUr44H4B59ct59tkO1jcxYYFgPLd37152795Nv379mptkEtVVV10FwIYNG6ivr/e5NLG14cNa7mS58+BrCRYIAF/8IuA0G/34x1Be7nN5eiELBOM5t3Zw0003kZLiyQzrMZOTk8OkSZOora3lr3/9q9/Fianw2+/Ql2pOjJ0BkdleE8oNN0BWFjPYTM7pffyP/+F3gXofCwTjOTcQbr75Zp9L0jlf+cpXAHjuued8LknsNDXB5TteBSB4bwLWDgDS0uDWWwH4augNXn0V3n/f5zL1MqKauD36BQUFWlRU5HcxTBc0NDSQk5NDdXU1+/btY/To0X4XqUMHDhxgzJgxpKSkUF5eTm5urt9F8tyW9WcYf/VgsqiBPXtg7Fi/i9S+116Dr3yFEwPGMPj4TsZPSmXTJkhP97tgyUNENqpqQXdeazUE46l3332X6upqJk6cmBRhADBq1CjmzZtHQ0MDL774ot/FiYkjjz9LFjV8mntl4oYBOP0IkybR//g+vjfkN3zyCSxc6NRwTOxZIBjPhMNhHn30UQDuv/9+n0vTNQ899BAAS5cuJZFrzd1SUcFV7/5vAKq+85jPhelAMAj/2ynro/oEg/vW8vrr8E//5Fx22cSWBYLxzMsvv8zmzZvJy8vj4Ycf9rs4XXLrrbcyYsQISktLWbVqld/F8VT1Pz9KdvgUK4K3Mm3JrX4Xp2N33QXTphE6Us7a+5aSlgb/+Z/wve9ZTSHWLBCMJ+rq6pprB48//jgZGRk+l6hrUlJSeOCBBwB4+umne841Ej7+mKxlz9NACiu+8P8lR1t8IACPPw7AJct+yP/91TFE4Cc/gblzYf9+n8vXg1kgXMSpU7Bpk9PP9ZvfwPPPwwsvwLvvwu7d0MOHrXfJL3/5S/bv38+UKVP4xje+4XdxuuXBBx8kMzOTP//5z/zLv/xL8jcd7d8P992HqPIz/pnZ35zod4k6b/58KCyEigru+H/nsOqVCoYNgw8/hOnTnValigq/C9kDqWrUN+AWYCdQCixp53kBnok8XwzM7Mx2r7jiCo2XcFi1pER16VLV++5THT9e1Wm1vPAtPV31uutUf/AD1VWrVOvq4lbchHH27Fn9t3/7NxURBfSPf/yj30WKynvvvaepqakK6GOPPRa3/Z4+rfrnP6s+9ZTqQw+p3nab83f405+qrl3r/H12yXvvqebkqILu4hLNTT2h1dWxKHkMlZWpTprk/LNNmKBVfy/VBQta/v9SU1W/+lXV555TLS3txjHqoYAi7eZnedTDTkUkCOwCPg+UARuAe1S1pNU6twLfBW4FrgR+pqpXdrTtS4eO1KUP/CupqWFCoTDBYJsVOlP2Vuu0fq91dUE++yyd/fuzOLA/k/37M6k5EzrnpaGUJgYOrCcnp44+GQ0AhMPCyZOpHDuWyomTaeesn5HWyLhxp7j00lNMmHCK/v0bLlqe9our561z6lSIo0fTqaxM58yZFM6eTaG+LoCqM3txMBgmNTVMamoTaWnOz9TUMGmRn6mpjZHlYcJhoalJaGx0fjY1CuFwy+NwOEBjI6gKgYASEEVECQQUCUBAQCRMTU01x49XsWPHNo4cOYIQYO7cG7jhhhsAaYlNAAUNO+8pEHDK3LItCIed9bRJz4td9zXu6wKiBAKR10TWaXtfcNYJBltu7n6amiK3Rm253+S8NhSClBQo/XQH/3fZ74F6cgcNYMqUieTljaSpKZWmphCNjSk0NKTQ2BikoT5AQ4NT0U5JUVJSlOy+DfTrW0+/fg1kZjbS3gTTjY1w+HAf9u3LZN/eTMrLM2gKX7jCnje8hptuOsyIEWcvuI6okrlvH/23bqXfjh1IOMyO0dfyuf3LmTS7iSeeSL5pREPHjzNjyRKyP/0UFeHk5MkUj57Lqt2XsXF3HmfJQCNHOD2tiZycWgYMrKdPRhPp6U1kRH6mpYUJBsMtf0eBVn/Xkb9H9/fUMiO4tnl84efPn0X83Ne2/dnWBfcp5/y44HZb7/f2n/5zt4edelE7mA2saPX4e8D32qzzK5yQcB/vBIZ1tO0rOvqKbje72a3dWwPo46AB3o0s+oZGPk2S7tYP9PegZxPguCbDjShqCF7MKzACONjqcRlOLaCjdUYAh9tuTEQWAYsALiWT96UANIBeoLtD2/3+1Zl1wgSDdaSknCaUUkMwdJpgoO7i2+ngYiLhplTqG/pRX9+Phoa+KG2rNE5pkDAikeESGkA1eN77a11mkUZSgmcJBmsJBOqRQFPL6yPrajiIEkQ10HI/HCQc+dmyD0VQRMIgzn0VRWh5jLj/iwIq7ivcgwAIIgECgSCBYAqpqSGQgLNG22800nJfRSDyF0urv14h8k1HWl7Q+lCrtlpfIYw0fytr9ZKWbUTW0TCEFTRSA3FqGK1qJwFpvg9ODcK9NYUh3KTUNYRpaGgiHA4TkCZEGiHyU2hCAk1AOPJ7cH6X4aY0msLpNDWlo5py3u/T/a0Fg2dJTT1OauoxUkMnkED7VwtTTaGmZhQ1Z0ahBElJqaZ//2KCwfP/Xo+kpVHcvz/b+/XjYPVVhDd/nmDwLFdffZJQ6IZ2t58MngN+19hIYVUV00+cIKuxkczGRkKtOv81nEJTOI1wONX5ez/nFsD9e3aI+1078jdO8/LOc/8nOrneBZz/8o7K0OZ5bf2MQrj7U7l7EQjtlb7te+zMOs5C1aXAUnDOVL6xaHVkOdTWQk0NnDnj/Ozo/pkzThNAair06QNDhsDgwTBunHNuTiyn2WlogC1bnAuSbNjgdEKXlsLRo7R892klMxNGjoQxY5yrGk6a1PKzB54422ucOgUHDkB1NdTVOWEzciSMGuXM1NAVO3fC7bc7f0fDQvDWW1DQTsPAPTh/f9OnO4+ffDKDJUuWR/tWTLKI4ip4XnwklgEjWz3OAw51Y52LEoGMDOeWDEIh55+17T9sXZ3z4XDqVMt76tMHsrMT4mqGxmN9+8LUqd5sa+JEZwrrL38Z1qyB666D//ov53Fbv/gF7NgBl1wC//qv3uzf9HxeDDvdAEwQkbEikgrcDbzVZp23gG+I4yrgpKqe11zUG6SlwaBBLbWUoUOdDw0LA9MZOTnOsOdvfxvOnnXO4XryyXOHQO/dC//rfzn3f/rTrtdETO8VdSCoaiPwMLAC2AEsU9XtIrJYRBZHVnsb2IMz7PQ54DvR7teY3io1FX79a/iP/3C+SPzP/wmjR8Njj8E3v+nMbH3qFNx2m3MzprNstlNjktif/gSPPALbt7csCwScmsMzzzj9ZqZ3iWa208S+eokx5qJuu825hMCaNfDii9C/P3z3uzB+vN8lM8nIAsGYJCcCc+Y4N2OiYXMZGWOMASwQjDHGRFggGGOMASwQjDHGRFggGGOMASwQjDHGRFggGGOMASwQjDHGRFggGGOMASwQjDHGRFggGGOMASwQjDHGRFggGGOMASwQjDHGRFggGGOMASwQjDHGRFggGGOMASwQjDHGRFggGGOMASwQjDHGREQVCCIyUETeE5HdkZ8DLrDeCyJyRES2RbM/Y4wxsRNtDWEJ8L6qTgDejzxuz4vALVHuyxhjTAxFGwgLgP8Tuf9/gDvbW0lVPwCORbkvY4wxMZQS5euHqOphAFU9LCKDoy2QiCwCFkUe1lkzU7NBwFG/C5EA7Di0sGPRwo5Fi4ndfWGHgSAifwGGtvPUD7q704tR1aXA0si+i1S1IBb7STZ2LBx2HFrYsWhhx6KFiBR197UdBoKq3nSRHVeIyLBI7WAYcKS7BTHGGOOvaPsQ3gLui9y/D3gzyu0ZY4zxSbSB8GPg8yKyG/h85DEiMlxE3nZXEpHfA+uBiSJSJiL3d3L7S6MsX09ix8Jhx6GFHYsWdixadPtYiKp6WRBjjDFJys5UNsYYA1ggGGOMifA9EETkFhHZKSKlInLemc7ieCbyfLGIzPSjnPHQiWOxMHIMikVknYhM96Oc8dDRsWi13iwRaRKRu+JZvnjqzLEQkTkisllEtovImniXMV468T/ST0T+KCJbIsfiW36UM9Y6mg6o25+bqurbDQgCnwLjgFRgCzC5zTq3An8GBLgK+JufZfb5WFwNDIjcn9ebj0Wr9VYCbwN3+V1uH/8u+gMlwKjI48F+l9vHY/F94CeR+7k4MySk+l32GByL64CZwLYLPN+tz02/awiFQKmq7lHVeuAVnOkwWlsA/FYdHwH9I+c89DQdHgtVXaeqxyMPPwLy4lzGeOnM3wXAd4E/0LPPf+nMsbgXeF1VDwCoak89Hp05Fgpki4gAWTiB0BjfYsaedjwdULc+N/0OhBHAwVaPyyLLurpOT9DV93k/zjeAnqjDYyEiI4AvAs/GsVx+6MzfxaXAABFZLSIbReQbcStdfHXmWPwCuAw4BGwF/llVw/EpXkLp1udmtHMZRUvaWdZ2HGxn1ukJOv0+RWQuTiBcE9MS+aczx+KnwCOq2uR8GeyxOnMsUoArgBuBDGC9iHykqrtiXbg468yx+AKwGbgBGA+8JyIfquqpGJct0XTrc9PvQCgDRrZ6nIeT7F1dpyfo1PsUkWnAr4F5qloVp7LFW2eORQHwSiQMBgG3ikijqi6PSwnjp7P/I0dVtQaoEZEPgOlATwuEzhyLbwE/VqchvVRE9gKTgL/Hp4gJo1ufm343GW0AJojIWBFJBe7GmQ6jtbeAb0R6za8CTmpkhtUepsNjISKjgNeBr/fAb3+tdXgsVHWsqo5R1THAa8B3emAYQOf+R94ErhWRFBHpA1wJ7IhzOeOhM8fiAE5NCREZgjPz5564ljIxdOtz09cagqo2isjDwAqcEQQvqOp2EVkcef5ZnBEktwKlwBmcbwA9TiePxWNADvCfkW/GjdoDZ3js5LHoFTpzLFR1h4i8AxQDYeDXqtrjpo3v5N/FE8CLIrIVp9nkEVXtcdNiR6YDmgMMEpEy4P8BQhDd56ZNXWGMMQbwv8nIGGNMgrBAMMYYA1ggGGOMibBAMMYYA1ggGGOMibBAMMYYA1ggGGOMifj/ATxM5EHNwuDpAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "C, ALPHA, X, u, LA = sp.symbols('C, ALPHA, X, u, LA')\n", "c = 0.3\n", "alpha = 0.5\n", "\n", "def init(x):\n", " middle, width, height = 0.4, 0.1, 0.5 \n", " out = np.zeros_like(x)\n", " mask = (np.abs(x%1-middle)<=width)\n", " out[mask] = height/width**10 * (x[mask]%1-middle-width)**5 * \\\n", " (middle-x[mask]%1-width)**5\n", " return out\n", " \n", "def solution(t, x):\n", " return init(x - c*t)*math.exp(-alpha*t)\n", "\n", "dico = {\n", " 'box': {'x': [0., 1.], 'label': -1},\n", " 'space_step': 1./128,\n", " 'scheme_velocity': LA,\n", " 'schemes': [\n", " {\n", " 'velocities': [1,2],\n", " 'conserved_moments': u,\n", " 'polynomials': [1, LA*X],\n", " 'relaxation_parameters': [0., 2.],\n", " 'equilibrium': [u, C*u],\n", " 'source_terms': {u: -ALPHA*u},\n", " },\n", " ],\n", " 'init': {u: init},\n", " 'parameters': {\n", " LA: 1., \n", " C: c, \n", " ALPHA: alpha\n", " },\n", " 'generator': 'numpy',\n", "}\n", "\n", "sol = pylbm.Simulation(dico) # build the simulation\n", "viewer = pylbm.viewer.matplotlib_viewer\n", "fig = viewer.Fig()\n", "ax = fig[0]\n", "ax.axis(0., 1., -.1, .6)\n", "x = sol.domain.x\n", "ax.plot(x, sol.m[u], width=2, color='k', label='initial')\n", "while sol.t < 1:\n", " sol.one_time_step()\n", "ax.plot(x, sol.m[u], width=2, color='b', label=r'$D_1Q_2$')\n", "ax.plot(x, solution(sol.t, x), width=2, color='r', label='exact')\n", "ax.legend()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## A source term depending on time and space\n", "\n", "If the source term $S$ depends explicitely on the time or on the space, we have to specify the corresponding variables in the dictionary through the key *parameters*. The time variable is prescribed by the key *'time'*. Moreover, sympy functions can be used to define the source term like in the following example. This example is just for testing the feature... no physical meaning in mind !" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYQAAAD8CAYAAAB3u9PLAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAqvUlEQVR4nO3de3xcdZ3/8ddnkjRJk96TXtPSFErTFtpa0nIRkUuRi6vF/cH+wK6LLrXWXdjd32/1YX/iIg90V3xYfYjoPrAgVh+CyIoKClqQcm+RttB7aRt7TW80oS1t7pN8f3+cmWaSJs1lzsw5k7yfj8c8zsycM+d852Qy7/l+v+d8jznnEBERiQRdABERCQcFgoiIAAoEERGJUSCIiAigQBARkRgFgoiIAD4Fgpldb2bbzazCzJZ0ssyVZrbezLaY2St+bFdERPxjyZ6HYGZZwA7gWqASWAPc5pzbmrDMUGAVcL1zbp+ZjXTOvZfUhkVExFd+1BDmAhXOuV3OuUbgCWB+u2U+DfzGObcPQGEgIhI+2T6sYxywP+FxJXBxu2XOB3LM7GVgEPCAc+7nHa3MzBYBiwAKCgouKisr86GIIiL9w7p166qcc8W9ea0fgWAdPNe+HSobuAi4BsgHVpvZm865HWe80LllwDKA8vJyt3btWh+KKCLSP5jZ3t6+1o9AqATGJzwuAQ52sEyVc64GqDGzV4GZeH0PIiISAn70IawBJptZqZkNAG4Fnmm3zNPAR8ws28wG4jUpbfNh2yIi4pOkawjOuaiZ3QmsALKAR51zW8xscWz+Q865bWb2J2Aj0AI84pzbnOy2RUTEP0kfdppKHfUhNDU1UVlZSX19fUClCkZeXh4lJSXk5OQEXRQRCTEzW+ecK+/Na/3oQ0iryspKBg0axMSJEzHrqD+773HOUV1dTWVlJaWlpUEXR0T6qIwbuqK+vp4RI0b0mzAAMDNGjBjR72pFIpJeGRcIQL8Kg7j++J5FJL0yMhBERMR/CoReuOyyy7pcZuHChWzd6g3n9F//9V89fn1hYWHvCici0ksZd5TRtm3bmDp1akAl6p3CwkJOnTqV9Gsy8b2LSHolc5SRagi9EP/1/vLLL3PllVdy8803U1ZWxoIFC4gH7JVXXsnatWtZsmQJdXV1zJo1iwULFrR5/alTp7jmmmuYPXs2F154IU8//XQwb0hEhAw87DRRqjpae1Jreuedd9iyZQtjx47lwx/+MG+88QaXX3756fn3338/P/zhD1m/fv0Zr83Ly+O3v/0tgwcPpqqqiksuuYRPfvKT6kAWkUCohpCkuXPnUlJSQiQSYdasWezZs6fbr3XO8dWvfpUZM2Ywb948Dhw4wJEjR1JXWBGRs8joGkIY+j9yc3NP38/KyiIajXb7tY899hhHjx5l3bp15OTkMHHiRJ1rICKBUQ0hDXJycmhqajrj+RMnTjBy5EhycnJ46aWX2Lu316PWiogkTYGQBosWLWLGjBmnO5XjFixYwNq1aykvL+exxx5DFwMSkSDpsNMM0p/fu4h0jw47FRGRpCkQREQEUCCIiEiMAkFERAAFgoiIxCgQREQEUCCIiEiMAkFERAAFQq/9+Mc/ZsyYMcyaNYuZM2dyyy23sHv37rO+5qmnnuLiiy9m5syZlJeXs2LFijSVVkSka74Egpldb2bbzazCzJZ0MP9KMzthZutjt3v82G6QNm7cyH333cf69evZsGED11xzDX/7t3/b6YB7jz/+OEuXLuXpp59mw4YN/PKXv+T222+nsrIyzSUXEelY0oFgZlnAj4AbgGnAbWY2rYNFX3POzYrd7kt2u0HbtGkTF1xwwenHixcv5vDhw+zfv/+MZWtqaliyZAlPPvkko0ePBmDy5MlceeWVvPjii2krs4jI2fhRQ5gLVDjndjnnGoEngPk+rLdLZqm5dcfmzZuZPn16m+fy8/M5duzYGcs+8cQTzJ49m/Hjx7d5Pjc3l9ra2l6///7IOUdNTU3QxRDpk/wIhHFA4s/iythz7V1qZhvM7I9mNr2D+QCY2SIzW2tma48ePepD8fy3f/9+Bg0axODBg08/19TUxKFDhwC44447uPnmm0/P27x5MzNnzjxjPRs2bKCsrIzf/e53fP7zn2f+/Pk8//zzqX8DGWr79u3MmTOHUaNGsWnTpqCLI9Ln+BEIHf2mbt+Q/jZwjnNuJvAg8LvOVuacW+acK3fOlRcXF591w86l5taVjRs3nlE7+OlPf8rVV1/NzJkz+clPftJm3uDBg2lsbGzz3OrVq6mpqeGjH/0oN910Ew8//DDLly/nV7/6VdcF6Id+8pOfMHv2bNatW0dNTQ0PP/xw0EUS6XP8CIRKILEtpAQ4mLiAc+4D59yp2P3ngBwzK/Jh24Fo33/w/PPP861vfYulS5d2uPzHP/5xnnzySeI1nh07drBw4UIeffRRIpHWP8E3v/lN/vmf/zm1hc9AK1euZOHChdTW1jJv3jwAfvWrX/Xo6nQi0jU/LqG5BphsZqXAAeBW4NOJC5jZaOCIc86Z2Vy8IKr2YduB2LRpEy+//DIvvvgizjmmTp3Kn/70J6ZMmdLh8nPnzuVrX/sa8+bNo6GhgebmZn7+859z6aWXAl67+JIlS7jhhhuYPXt2Ot9KRnj22WcBuOuuu3jggQcoKytjx44drFy5ko997GMBl06k70i6huCciwJ3AiuAbcCTzrktZrbYzBbHFrsZ2GxmG4AfALe6MF+ZpwuPPfYYBw4cYN26dbz99ts89thjp8OgurqaxYsX88477/Ctb33r9Gs+97nPsWHDBl555RUGDBhAQUHB6XkPPvggf/7zn/n1r3/NQw89lPb3E3avvvoqAJ/4xCcwMz79ae/3xuOPPx5ksUT6HF0xLYP0x/d+6tQphg4dCsDx48cpLCxkx44dTJkyhUGDBnHkyBHy8/ODLaRIiOiKadJnrV69mubmZmbPnk1hYSEA559/PuXl5Zw8eZI//OEPAZdQpO9QIEioxZuLrrjiijbPL1iwAPCa70TEHwoECbXXXnsNgI985CNtno+f57Fy5cpOhwsRkZ7JyEDoj18A/fE9NzQ08OabbwJw+eWXt5lXUlJCUVERJ0+e1HhQIj7JuEDIy8ujurq6X31BOueorq4mLy8v6KKk1dq1a2loaGD69OmMGDHijPnxkwO3bNmS7qKJ9El+nIeQViUlJVRWVhLWYS1SJS8vj5KSkqCLkVadNRfFTZ8+nVdeeYUtW7Zw/fXXp7NoIn1SxgVCTk4OpaWlQRdD0iDeodxZIEyb5g2qqxqCiD8yrslI+gfnHKtXrwbOXkMA2Lp1a9rKJdKXKRAklI4ePcrx48cZOnRop01liYHQn/qURFJFgSChVFFRAcB5552HdXKRiuLiYoqLizl58mSHFyYSkZ5RIEgo7dy5E/AC4WzUjyDiHwWChFK8hjB58uSzLqd+BBH/KBAklLpbQ9C5CCL+USBIKPW0hqBAEEmeAkFCxznXplP5bOJ9CDrSSCR5CgQJnaqqKk6cOMHgwYMpKjr7lVbjRxqdOnVKRxqJJEmBIKGT2FzU2SGnidRsJOIPBYKETnc7lOMUCCL+UCBI6HS3QzkuvtyuXbtSViaR/kCBIKHT3Q7luHPOOQeAvXv3pqxMIv2BAkFCp6dNRgoEEX8oECRUnHOnA6G7TUaJgaBDT0V6z5dAMLPrzWy7mVWY2ZKzLDfHzJrN7GY/tit9T3V1NSdOnGDQoEEUFxd36zXDhg2jsLCQU6dOcezYsRSXUKTvSjoQzCwL+BFwAzANuM3MpnWy3LeBFcluU/qu7oxy2p6Zna4l7Nu3L2VlE+nr/KghzAUqnHO7nHONwBPA/A6Wuwt4CnjPh21KH9XTI4zi1I8gkjw/AmEckHiKaGXsudPMbBzwKeChrlZmZovMbK2Zre1v102WnncoxykQRJLnRyB0VK9v37P3feArzrnmrlbmnFvmnCt3zpV3tw1Z+o49e/YAMGnSpB69bsKECYACQSQZ2T6soxIYn/C4BDjYbply4IlYm3ARcKOZRZ1zv/Nh+9KHxMcjGj9+fBdLtqUagkjy/AiENcBkMysFDgC3Ap9OXMA5Vxq/b2bLgT8oDKQjlZWVAJ1eR7kzCgSR5CUdCM65qJndiXf0UBbwqHNui5ktjs3vst9ABLxzEFRDEAmOHzUEnHPPAc+1e67DIHDOfdaPbUrfU11dTX19PUOGDGHQoEE9eu2YMWPIycnh6NGj1NbWMnDgwBSVUqTv0pnKEhrx5qKe1g4AIpHI6dfpXASR3lEgSGjEm4t62n8Qp5PTRJKjQJDQ6G3/QZz6EUSSo0CQ0FAgiARLgSCh0dtDTuN0cppIchQIEhqqIYgES4EgoaFAEAmWAkFCoaWlJekmo3iQHDhwgGg06lvZpC1dg6jvUiBIKFRVVdHY2MiwYcMoKCjo1Tpyc3MZM2YMzc3NHDhwwOcSCkBzM8ybB+efD2+8EXRpxG8KBAmFZJuL4uK1CwVCavz0p7ByJezcCR/9KCxdqhpDX6JAkFDwKxDGjfMuxaFA8N+pU3DPPd79667zagtf/jJ8/evBlkv8o0CQUFAghN93vwuHDsGcOfDHP8ITT3jPL1vmhYNkPgWChEKyHcpxCoTUOHQIvvMd7/7SpWAGf/d3UFoKR46oP6GvUCBIKKiGEG4PPgg1NTB/PlxxhfecGdx8s3f/qaeCK5v4R4EgoaBACLeVK73pF7/Y9vnEQGhpSW+ZxH8KBAkFBUJ41dTAunUQicCll7adN2cOjB8PBw7AX/4STPnEPwoECVxLS8vpL/D4F3pvJQaC0/GQvnjrLYhGYdYsGDy47bzEZqNf/zrtRROfKRAkcEeOHCEajVJUVER+fn5S6xo0aBCDBg2ivr6eY8eO+VTC/u2117zp5Zd3PD8xEJTBmU2BIIHz6wijODUb+SseCB/5SMfzL7kExo6Fffu8piXJXAoECZxfzUVxOlvZP9EorF7t3e+shhCJwPXXe/d1+GlmUyBI4A4ePAj4Fwjx9cRrHtJ769d7ncrnnQejR3e+3MUXe9O33kpLsSRFFAgSOL9rCGoy8k9X/Qdxc+d6UwVCZvMlEMzsejPbbmYVZrakg/nzzWyjma03s7Vm1sXHS/qT+Bf32LFjfVmfAsE/r7/uTTvrP4ibPh3y86GiAt5/P/XlktRIOhDMLAv4EXADMA24zcymtVvsRWCmc24W8I/AI8luV/oO1RDCybmuO5TjcnJg9mzv/po1qS2XpI4fNYS5QIVzbpdzrhF4ApifuIBz7pRrPSi8ANDBaXJaqvoQFAjJqaiAo0dh5EivD6Er8WYjnaCWufwIhHHA/oTHlbHn2jCzT5nZu8CzeLUEEUBNRmG1ebM3nT3bOwGtK+pHyHx+BEJHH5UzagDOud8658qAm4BvdLoys0Wxfoa1R48e9aF4EmY1NTWcOHGC3NxcRowY4cs6R44cSVZWFlVVVTQ0NPiyzv5o61ZvOn1695ZPDASdoJaZ/AiESiBxAJoS4GBnCzvnXgXONbOiTuYvc86VO+fKi4uLfSiehFm8uWjs2LFYd36GdkNWVhZjxoxps37puXggTGvfI9iJ0lIYMcJrZtq7N3XlktTxIxDWAJPNrNTMBgC3As8kLmBm51nsv93MZgMDgGofti0Zzu/mojg1GyWvp4FgpmajTJd0IDjnosCdwApgG/Ckc26LmS02s8Wxxf4XsNnM1uMdkfS/nUYeE/w/wihOgZCc5mZ4913v/tSp3X+dAiGzZfuxEufcc8Bz7Z57KOH+t4Fv+7Et6Vv8PsIoToGQnD17oL4exo2DIUO6/zoFQmbTmcoSKDUZhVNPm4vi5szxpm+/rQvmZCIFggRKTUbh1NtAKC72xjyqqfFqGZJZFAgSKDUZhdOWLd60p4EAMGOGN9240b/ySHooECRQqa4haMTT3ultDQHgwgu96aZN/pVH0kOBIIFpaWlpcx6Cn+KBcPDgQV1Ks4daWmDbNu++AqF/USBIYKqrq2lqamLYsGFJXzqzvYKCAoYOHUpjYyNVVVW+rruv27cPamu9voDhw3v+ejUZZS4FggQmVc1FcepH6J1kmovAO28hKwt27oS6Ov/KJamnQJDApOqQ0zgFQu8kGwh5eTB5ctumJ8kMCgQJTKqOMIpTIPROsoEAajbKVAoECYyajMIp/qu+J0NWtKeO5cykQJDAqMkonHbu9KZTpvR+HQqEzKRAkMCoySh8jh+H6moYONA7yqi31GSUmRQIEhg1GYXPX//qTc89t3tXSevMOedAYSEcOeJdH0EygwJBAhM/i7ikpCQl61cg9FxFhTc999zk1hOJwAUXePfVbJQ5FAgSiLq6Oqqrq8nJySFVV8YrLi4mJyeH999/nzodEN8t8RrCeeclvy41G2UeBYIEIrG5KBJJzccwEonoUpo9FK8h+BEI6ljOPAoECUSqm4vi1GzUM4l9CMlSIGQeBYIEQoEQTqmoIWze7F2SU8JPgSCBSHcgaBjsrtXWwsGDkJMD48cnv77hw71LcNbVwa5dya9PUk+BIIFIVyDE168aQtfiX9qlpd7gdH5Qs1FmUSBIIOKBMN6Pn6JnoSaj7vOzuShORxplFgWCBEJ9COHj1zkIiVRDyCy+BIKZXW9m282swsyWdDB/gZltjN1WmdlMP7YrmWv//v2AAiFM/DwHIU6BkFmSDgQzywJ+BNwATANuM7P2A+fuBj7qnJsBfANYlux2JXM1NDTw3nvvkZWVxahRo1K6rfjAeQcPHqSlpSWl28p0qWgyKiuD7Gxv3TU1/q1XUsOPGsJcoMI5t8s51wg8AcxPXMA5t8o5dyz28E0gtT8LJdQSr6Oc5VfvZSfy8/MZPnw40WiUoxpU56xS0WSUm+uNmupc63UWJLz8CIRxwP6Ex5Wx5zpzB/BHH7YrGSpd/QdxajbqWmOjdy3lSAQmTvR33Wo2yhx+BEJHYyK6Dhc0uwovEL7S6crMFpnZWjNbq190fZMCIXz27PEueTlhgver3k860ihz+BEIlUDisYMlwBkDx5jZDOARYL5zrrqzlTnnljnnyp1z5aka9EyCle5AiG8n3pEtZ/JzyIr2VEPIHH4EwhpgspmVmtkA4FbgmcQFzGwC8BvgM865HT5sUzJYugPhnHPOAWDv3r1p2V4mSkWHclw8EDZu9PoSJLyyk12Bcy5qZncCK4As4FHn3BYzWxyb/xBwDzAC+G/zrroRdc6VJ7ttyUxBBcK+ffvSsr1MlIoO5bgJE2DwYKiq8i6Yk8yV2CS1kg4EAOfcc8Bz7Z57KOH+QmChH9uSzJfuQJgwYQKgGsLZpOIchDgzr5bwxhtes5ECIbx0prKknWoI4ZPKJiNQP0KmUCBIWjU1NXHo0CHM7PTFa1ItfhGegwcP0tjYmJZtZpLm5taB7SZNSs02FAiZQYEgaXX48GGcc4wePZqcnJy0bDMnJ4exY8finNMw2B2orISmJhgzBgoKUrMNHXqaGRQIklbpGuW0PTUbdS7VzUUAF1zgTbdu1cVywkyBIGmV7v6DOB162rlUHmEUN3Sod9Gd+vrW7Un4KBAkreK/0NMdCDrSqHOpPMIokZqNwk+BIGm1K9Z7WVpamtbtqsmoc+moIYA6ljOBAkHSavfu3UBwgaAawpnS0YcACoRMoECQtIoHwqRUHd/YCTUZdcy51I5jlEhNRuGnQJC0aWlpCbyGsG/fPpwG1Dnt8GGorYXhw2HYsNRua8oUyMnxznk4dSq125LeUSBI2hw+fJiGhgaKioooLCxM67YLCwsZPnz46au1iSddHcrghUFZmXd/y5bUb096ToEgaRNU7SBO/QhnSlf/QZyajcJNgSBpEz/CKN39B3HxfgQdadQqXUcYxaljOdwUCJI2qiGETzqbjKC1hrBhQ3q2Jz2jQJC0USCET7prCBdd5E3XrdMQFmGkQJC0CeqktDg1GZ0p3X0II0fCxIlQU6OO5TBSIEjaBHUOQpxqCG1VVcHx41BY6H1Rp8vFF3vTv/wlfduU7lEgSFo0NjZSWVlJJBI5/Us93RQIbW3f7k2nTPGuapYuCoTwUiBIWuzduxfnHCUlJWm7DkJ7xcXFFBYWcuzYMd5///1AyhAmiYGQTgqE8FIgSFoE3VwEYGacf/75AGyPfxv2Y/FdENslafOhD0F2tteHcPJkerctZ6dAkLQI+gijuCmxn8PvvvtuoOUIgx07vGm6awj5+TBzpjeO0tq16d22nJ0CQdIi6COM4uKBoBpCcE1GoGajsFIgSFqEpYZQFhtMp78HQjTaeshpupuMQIEQVr4Egpldb2bbzazCzJZ0ML/MzFabWYOZfcmPbUpmCUMfAqiGELdnDzQ1QUkJFBSkf/uJgaDBZ8Mj6UAwsyzgR8ANwDTgNjOb1m6x94F/AZYmuz3JTGFpMpo8eTIAFRUVRKPRQMsSpCCbiwAmT/aus3zoEMQusy0h4EcNYS5Q4Zzb5ZxrBJ4A5icu4Jx7zzm3BmjyYXuSYeKHeebn5zN69OhAy1JQUMD48eNpamo6XWvpj4IOhEgE5s717q9eHUwZ5Ex+BMI4YH/C48rYc71iZovMbK2ZrT169GjShZPgbYoNbTlt2jQsnWdAdULNRsEHAsAVV3jTF18MrgzSlh+B0NF/eK9bBZ1zy5xz5c658uLi4iSKJWERD4QL42MfB0yBEI5AuPZab/rCC8GVQdryIxAqgfEJj0uAgz6sV/qIsAWCjjRqPQchiCOM4i66yLts5+7drcNwS7D8CIQ1wGQzKzWzAcCtwDM+rFf6iLAFQn+vIXzwgdeZm5sLAQ0rBUBWFlx9tXdftYRwSDoQnHNR4E5gBbANeNI5t8XMFpvZYgAzG21mlcD/Bb5mZpVmNjjZbUv4OefYvHkzoEAIi3jtYPJk70s5SGo2CpdsP1binHsOeK7dcw8l3D+M15Qk/cz+/fv54IMPKCoqYtSoUUEXB4CSkhLy8/M5cuQIx48fZ+jQoUEXKa3C0H8QFw+ElSu9C+YEHVD9nc5UlpRKbC4KwxFGAJFIpF8PchemQJg0ybsdP65xjcJAgSApFbb+g7j+3LG8bZs3DUMggJqNwkSBICkV1kDoz6Oerl/vTWfNCrIUrRQI4aFAkJQKayDEy/P2228HXJL0OnnSG9RuwACIVZICd/XV3pnLq1ZBdXXQpenfFAiSMk1NTad/gU+fPj3g0rR12WWXAbB69WpaWloCLk36bNzoTadP90IhDIYN82oJ0Sj8z/8EXZr+TYEgKbN9+3aampqYNGkShYWFQRenjbFjx3LOOefwwQcfsHXr1qCLkzbvvONNw9JcFPf3f+9Nf/GLYMvR3ykQJGXC2lwUF68lrFq1KuCSpE+8/+BDHwq0GGe46SYYOBDeeANiA+NKABQIkjJhD4RLL70U6J+BELYaQmEhfOpT3v3HHgu2LP2ZAkFS5s033wRgVti+fWL6Ww2hqQliJ40zc2awZelIYrORLpoTDAWCpMSpU6d4/fXXiUQiXHXVVUEXp0MzZsxg4MCB7Ny5k/4w1Pq770JDA5x7LgwO4cAx8+bByJHe0Bo6SS0YCgRJiZdffpmmpibmzJnD8OHDgy5Oh3Jycpgbu0rL6n5wlZawNhfFZWfDbbd59x95JNiy9FcKBEmJFStWAHDdddcFXJKz60/NRmEPBIAvfAHM4Kc/hb17gy5N/6NAkJRQIIRPWI8wSjR1qldLaGqC//zPoEvT/5gLce9NeXm5W6vGxIyze/duJk2axJAhQ6iqqiI725dBdVOiurqaoqIi8vLyOHHiBAPCcraWz5yDESPg2DHvovbjen2R29Tbvh2mTfPOXt6xA0pLgy5RZjGzdc658t68VjUE8V28djBv3rxQhwHAiBEjKCsro76+ntdffz3o4qTM/v1eGBQVwdixQZfm7KZM8Y44ikbhm98MujT9iwJBfBcPhI997GMBl6R7brnlFgAefvjhgEuSOvE+89mzvTb6sPuP//CujfCzn7U2dUnqKRDEV01NTbz44otA+PsP4hYuXIiZ8dRTT/XZw0/jI4lec02w5eiu886DxYu9i+bccot32U9JPQWC+Or555/n5MmTTJkyhXPOOSfo4nTLhAkTuOGGG2hqamL58uVBF8d3zsHzz3v340NNZ4LvfAdmzPBGZ/3853WyWjooEMQ3LS0tfO1rXwPgjjvuCLg0PfOFL3wBgGXLlhHmAy16Y8cOrw+huDicZyh3Jj/fG/20sBCefBJ+8IOgS9T3KRDEN48//jjr16+npKSEO++8M+ji9MiNN97IuHHjqKio4KWXXgq6OL5KbC6KZNh//PnnQ7xr59/+Db7+ddUUUinDPh4SVg0NDadrB/fddx/5+fkBl6hnsrOzWbhwIQBLly7tU9dIiDcXZUgf/xluvRUefNALs/vugwULvAv9iP90HsJZnDgBf/2rd/vgA2hp8Y7QGD/eOzZ64sTwXGQkaN/73vf493//d6ZPn86GDRvIysoKukg9duDAAaZMmUJNTQ133XUXDzzwAJYJh+ScRVOTd/7ByZNes1FJSdAl6r1nn/XC4dQp7/DZr34VvvhFyMsLumThksx5CDjnkr4B1wPbgQpgSQfzDfhBbP5GYHZ31nvRRRe5dGlpcW7bNucefti522937txznfMqp53fcnOdu+IK5+6+27mXXnKusTFtxQ2Nuro696UvfcmZmQPc73//+6CLlJQXXnjBDRgwwAHunnvuScs245+9b3/buU99yrlZs5wbNsy5KVOc+8d/dG75cufq6nq37tde8z6rZWX+ljkoGzY4d+mlrf+DI0Z4/6+/+Y1zhw8HXbpwANa6Xn6XJ11DMLMsYAdwLVAJrAFuc85tTVjmRuAu4EbgYuAB59zFXa27tLTM3X//MvLymsnLaybZH52J77WuLsLu3YPYsmUomzcPZevWoZw40fbn/oABzYwdW8uYMXUMHtxEJOKIRiMcPZrLoUMDOXKkbbPIwIFRPvShasrLq5gzp4ri4oakyug9hurqXPbtK6CyspBjxwZQU5NNba13wpcZZGe3kJ/fTH5+lIEDo+TlNTNwYDN5edHYtIn8fO9xS0uEpiajqSlCNBpJmLY+F41GaGmBSMSRleWIRCArq4VIxHvOrIUTJ6o5dKiS119/hT179mCWxe23f5bbb/8sYB1GKBBbR+vNzKt5OedN298H73j0xFsk0jq/ubn1fuJrsrMhJ6d1Gol4v5bb36JRb+oc5OZ6t7VrX+Pee7+Kc3VMnDiGq666jOnTL6SpaSCNjXnU1+fS0DCAuroB1NV5fwszR26u93cYNaqeMWNqGTu2joKC5g7/zs3NxtatQ1i1qohVq4o5cGDgWT8XEybU8OUvb2Xq1J4df7l8+SR+8YtSbrppP3feuaNHrw0r5+AvfxnB8uWTqKhoO2zr0KGNjB9fw4gRjQwf3kBhYZTc3GZyc1tOT7OzWzDzPsvgfQa9W9vH4BLmeY/bzk88p8O1edzR/DPntf1f795rut7eokWX9LqG4EcgXArc65y7Lvb4/wE4576VsMyPgZedc7+MPd4OXOmcO3T2dZc7SGwyagBqEm61ndyPP64FsoABQCEwEhgNTAY6OiTyIPB6wm0j0PE/tGc48GHgo3iVpPbXDX4XWA28CWzFqyAd7mRdQ4AS4DxgKlCWMA3hWMXSTUeAvwLvA/V43XaT8f7OiT8oqoFngRV4v6/24X1GLwe+AEzB+yzeD/wH8S+Err0FzAE+AfwhubcSSlOBT+L91pyJ93/U3/W+yciPcQXGAfsTHlfi1QK6WmYccEYgmNkiYJF3/0Jyct6nuTmP5uY8IDd2S344ZbMmBg6sZOjQrQwZsplhwzaTl3c4IYFHAN0Zx78OsxXACurqRlJdPYeqqjm8//4smpvL8L7QP9dmu9nZtWRl1QHQ0jKA5uZ8mps774TNyTlOQUElBQX7yMurIju7huzsWrwvBaOlJZvm5oFEo956WqcD2zxubs7HrJlIpIlIpAmz6On73i2KmTf11h3Bufgtq800OzuX3NwC8vIGMmLEcLKystr8muroBm1//cdv7WsMiffBqwUk3lpaWmsKXu2l7Tqc8375x3/9R6Pe6xJrDO1vAI2NUF/vXTOgvh7q6x3Hj9dy/Hg90Wjj6f2elXWKSKSGrKxTZGXVkpVVixk0N+cSjRZQXz+aurpx1NWNpaVlFDCqw79rfn4lRUWrKCpaxZAhW4hE4h3Zg4ELYvc30Nz8b+ze/Vn27bsZuJsxYz5EWdn3zviF2d6xYxfyzjtzyM4+xWWXtZCdffVZl89ca4A1OAcNDcXU1pbQ2DicxsbhRKMFsf+xvNNT57JxzvBasi3hPgmP6WDemfPb8h53PP9s8zqan7hMR/M6n19TQ6/5EQgd9bq1/6R2ZxnvSeeWAcsg3qk8PPa8949aUwO1td60o1vivNpa7wtgwADveq2jRnkX4Jg0CSZNyiE7uxQoBT7e6zffmcZG2LDBGzJgzRrYudM7waa6OoempiE0NbX9JVNQ4HVWT5zojfg4dSqUlXnToqKhwFBavyQkPQwoiN16rqUFDh70Dko4eRLq6rzP8bnneodTDhpUAvxd7Na1F16A+fPh0KEbue66G3nkEc7ajBo/Ce3uuwu5995ne/UeJPMkcxyEH4FQCYxPeFyC1/bS02XOysw7miAvzztqIuwGDIA5c7xbovp678vhgw+8X7N5eV5YDR6cGWPMSPdFIt5RPX4d2XPttd6RNn/zN7B8uVfzWb6841BYtQr+/GcYNAj+9V/92b70fX6ch7AGmGxmpWY2ALgVeKbdMs8A/2CeS4ATXfUf9FV5ed4Zo+ee6x26OmYMDBmiMJDuueoq+OMfvRrlL34Bn/mMFwztfeMb3vRf/gWGDUtvGSVzJR0IzrkocCdeb9g24Enn3BYzW2xmi2OLPQfswutVfRj4p2S3K9JfXXEFrFjhDenwy196F5Q5dqx1/lNPwZ/+5M3/P/8nuHJK5tGJaSIZ6s034brrvObHQYO8y09u2uSFBcDdd+t6Av2RLpAj0g9dcgm8+irMm+f1Sy1d6oXBkCHw3e/CvfcGXULJNOG+nJWInNXMmd7RR2+9Bd//vnck3d13e0M7iPSUAkGkD5g7Fx5/POhSSKZTk5GIiAAKBBERiVEgiIgIoEAQEZEYBYKIiAAKBBERiVEgiIgIoEAQEZEYBYKIiAAKBBERiVEgiIgIoEAQEZEYBYKIiAAKBBERiVEgiIgIoEAQEZEYBYKIiAAKBBERiVEgiIgIkGQgmNlwM3vBzHbGpsM6We5RM3vPzDYnsz0REUmdZGsIS4AXnXOTgRdjjzuyHLg+yW2JiEgKJRsI84Gfxe7/DLipo4Wcc68C7ye5LRERSaHsJF8/yjl3CMA5d8jMRiZbIDNbBCyKPWxQM9NpRUBV0IUIAe2HVtoXrbQvWk3p7Qu7DAQz+zMwuoNZd/d2o2fjnFsGLItte61zrjwV28k02hce7YdW2hettC9amdna3r62y0Bwzs07y4aPmNmYWO1gDPBebwsiIiLBSrYP4Rng9tj924Gnk1yfiIgEJNlAuB+41sx2AtfGHmNmY83sufhCZvZLYDUwxcwqzeyObq5/WZLl60u0LzzaD620L1ppX7Tq9b4w55yfBRERkQylM5VFRARQIIiISEzggWBm15vZdjOrMLMzznQ2zw9i8zea2ewgypkO3dgXC2L7YKOZrTKzmUGUMx262hcJy80xs2Yzuzmd5Uun7uwLM7vSzNab2RYzeyXdZUyXbvyPDDGz35vZhti++FwQ5Uy1roYD6vX3pnMusBuQBfwVmAQMADYA09otcyPwR8CAS4C/BFnmgPfFZcCw2P0b+vO+SFhuJfAccHPQ5Q7wczEU2ApMiD0eGXS5A9wXXwW+HbtfjDdCwoCgy56CfXEFMBvY3Mn8Xn1vBl1DmAtUOOd2OecagSfwhsNINB/4ufO8CQyNnfPQ13S5L5xzq5xzx2IP3wRK0lzGdOnO5wLgLuAp+vb5L93ZF58GfuOc2wfgnOur+6M7+8IBg8zMgEK8QIimt5ip57oeDqhX35tBB8I4YH/C48rYcz1dpi/o6fu8A+8XQF/U5b4ws3HAp4CH0liuIHTnc3E+MMzMXjazdWb2D2krXXp1Z1/8EJgKHAQ2Af/qnGtJT/FCpVffm8mOZZQs6+C59sfBdmeZvqDb79PMrsILhMtTWqLgdGdffB/4inOu2fsx2Gd1Z19kAxcB1wD5wGoze9M5tyPVhUuz7uyL64D1wNXAucALZvaac+6DFJctbHr1vRl0IFQC4xMel+Ale0+X6Qu69T7NbAbwCHCDc646TWVLt+7si3LgiVgYFAE3mlnUOfe7tJQwfbr7P1LlnKsBaszsVWAm0NcCoTv74nPA/c5rSK8ws91AGfBWeooYGr363gy6yWgNMNnMSs1sAHAr3nAYiZ4B/iHWa34JcMLFRljtY7rcF2Y2AfgN8Jk++OsvUZf7wjlX6pyb6JybCPwa+Kc+GAbQvf+Rp4GPmFm2mQ0ELga2pbmc6dCdfbEPr6aEmY3CG/lzV1pLGQ69+t4MtIbgnIua2Z3ACrwjCB51zm0xs8Wx+Q/hHUFyI1AB1OL9Auhzurkv7gFGAP8d+2UcdX1whMdu7ot+oTv7wjm3zcz+BGwEWoBHnHN9btj4bn4uvgEsN7NNeM0mX3HO9blhsWPDAV0JFJlZJfB1IAeS+97U0BUiIgIE32QkIiIhoUAQERFAgSAiIjEKBBERARQIIiISo0AQERFAgSAiIjH/H1UURwKiYsOOAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "t, C, X, u, LA = sp.symbols('t, C, X, u, LA')\n", "c = 0.3\n", "\n", "def init(x):\n", " middle, width, height = 0.4, 0.1, 0.5 \n", " out = np.zeros_like(x)\n", " mask = (np.abs(x%1-middle)<=width)\n", " out[mask] = height/width**10 * (x[mask]%1-middle-width)**5 * \\\n", " (middle-x[mask]%1-width)**5\n", " return out\n", " \n", "dico = {\n", " 'box': {'x': [0., 1.], 'label': -1},\n", " 'space_step': 1./128,\n", " 'scheme_velocity': LA,\n", " 'schemes': [\n", " {\n", " 'velocities': [1, 2],\n", " 'conserved_moments': u,\n", " 'polynomials': [1, LA*X],\n", " 'relaxation_parameters': [0., 2.],\n", " 'equilibrium': [u, C*u],\n", " 'source_terms': {u: -sp.Abs(X-t)**2*u},\n", " },\n", " ],\n", " 'init': {u: init},\n", " 'generator': 'cython',\n", " 'parameters': {LA: 1., C: c},\n", "}\n", "\n", "sol = pylbm.Simulation(dico) # build the simulation\n", "viewer = pylbm.viewer.matplotlib_viewer\n", "fig = viewer.Fig()\n", "ax = fig[0]\n", "ax.axis(0., 1., -.1, .6)\n", "x = sol.domain.x\n", "ax.plot(x, sol.m[u], width=2, color='k', label='initial')\n", "while sol.t < 1:\n", " sol.one_time_step()\n", "ax.plot(x, sol.m[u], width=2, color='b', label=r'$D_1Q_2$')\n", "ax.legend()" ] }, { "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.12" }, "widgets": { "application/vnd.jupyter.widget-state+json": { "state": {}, "version_major": 2, "version_minor": 0 } } }, "nbformat": 4, "nbformat_minor": 4 }