{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Poiseuille flow\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 a Poiseuille flow modeling by the Navier-Stokes equations." ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "%matplotlib inline" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## The $\\DdQq{2}{9}$ **for Navier-Stokes**\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", "\n", "* nine velocities $\\{(0,0), (\\pm1,0), (0,\\pm1), (\\pm1, \\pm1)\\}$, identified in pylbm by the \n", "numbers $0$ to $8$,\n", "\n", "* nine polynomials used to build the moments\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", "where $E = X^2+Y^2$.\n", "\n", "* three conserved moments $\\rho$, $q_x$, and $q_y$,\n", "\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", "\n", "* equilibrium value of the non conserved moments\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", "\\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", "with $p=\\rho\\lambda^2/3$." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Build the simulation with pylbm\n", "In the following, we build the dictionary of the simulation step by step.\n", "\n", "### The geometry\n", "\n", "The simulation is done on a rectangle of length $L$ and width $W$. We can use $L=W=1$.\n", "\n", "We propose a dictionary that build the geometry of the domain. The labels of the bounds can be specified to different values for the moment." ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "+----------------------+\n", "| Geometry information |\n", "+----------------------+\n", " - spatial dimension: 2\n", " - bounds of the box: [0. 1.] x [-0.5 0.5]\n", " - labels: [0, 1, 2, 3]\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAQoAAAD8CAYAAACPd+p5AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAEE5JREFUeJzt3X+MXWWdx/H3t5Riqohgh/KjnRazdWEAI3rp6iZQTVut3Q1FYvgRhJLgNiquka6rJGwEqn+47rrWhSYyYYUqgYqajROFdbWCGmOhg9RCMYWKKN0WKAqutVug9Lt/3Fschjt9rnNP771T3q9kcs85z5PzfOfOzGee+8yZeyIzkaT9mdTtAiT1PoNCUpFBIanIoJBUZFBIKjIoJBUZFJKKDApJRQaFpKLJ3S5gLNOmTcvZs2d3uwzpoHbvvfc+lZl9pX49GxSzZ89meHi422VIB7WI+HUr/XzpIanIoNC4bd68mQsvvJCTTjqJI444gqlTp3LiiSeyfPlytm/f3u3yVKGefemh3rd161a2b9/Oe9/7XmbMmMHkyZO5//77GRwcZM2aNWzYsIGjjz6622WqAgaFxm3+/PnMnz//ZcfPPPNMzj33XG666SY+8YlPdKEyVc2XHqrcrFmzAHj66ae7XImq4oxCbdu9ezc7d+5k9+7dPPjgg3zyk58EYPHixV2uTFVxRqG23XDDDfT19TFz5kze/e5388wzz3DzzTdzxhlndLs0VcQZhdp29tlnc+KJJ7Jz507uu+8+hoaG2LFjR7fLUoUmfFAcc8y/8sQTf+x2GXrRJOBtbNjwcS6//NuAs4pumD791Tz++McrO9+Ef+lhSPSiYxof67tdyCtW1T8XEz4o1Kv2AP/X7SJUEYNCbfjDGMd/BTwJzOhgLTqQJvwahbrpO9TD4gTgddRnEduAB4ApwLu6V5oqZVCoDacAPwc2An8EAjgCqAF/TT08dDAwKNSGUxofOti5RiGpyKCQVGRQSCoyKCQVGRSSigwKSUUGhaQig0JSkUEhqaiSoIiIRRGxOSK2RMQV++n3vojIiKhVMa6kzmg7KCLiEGAV8B5gALggIgaa9Dsc+Chwd7tjSuqsKmYUc4EtmflIZj4HrAGWNOn3aeBzwO4KxpTUQVUExfHAYyP2tzaOvSgiTgNmZua3KxhPUodVERTR5Fi+2BgxCfgC8A/FE0Usi4jhiBj2zVml3lFFUGwFZo7Yn0H93Uv2OZz6/yLfFRGPAm8DhpotaGbmYGbWMrPW11e8E7ukDqkiKNYDcyLihIiYApwPDO1rzMzfZ+a0zJydmbOBdcBZmTlcwdiSOqDtoMjMPcBHgO8CvwBuy8xNEbEiIs5q9/ySuq+Sd7jKzNuB20cd+9QYfd9RxZiSOscrMyUVGRSSigwKSUUGhaQig0JSkUEhqcigkFRkUEgqMigkFRkUkooMCklFBoWkIoNCUpFBIanIoJBUZFBIKjIoJBUZFJKKDApJRQaFpCKDQlKRQSGpyKCQVGRQSCoyKCQVGRSSigwKSUUGhaQig0JSkUEhqcigkFRUSVBExKKI2BwRWyLiiibtyyPiwYjYGBFrI2JWFeNK6oy2gyIiDgFWAe8BBoALImJgVLf7gFpmvgn4BvC5dseV1DlVzCjmAlsy85HMfA5YAywZ2SEz78zMXY3ddcCMCsaV1CFVBMXxwGMj9rc2jo3lUuCOCsaV1CGTKzhHNDmWTTtGvB+oAfPGaF8GLAPo7++voDRJVahiRrEVmDlifwawbXSniFgAXAmclZnPNjtRZg5mZi0za319fRWUJqkKVQTFemBORJwQEVOA84GhkR0i4jTgeuoh8WQFY0rqoLaDIjP3AB8Bvgv8ArgtMzdFxIqIOKvR7V+A1wBfj4gNETE0xukk9aAq1ijIzNuB20cd+9SI7QVVjCOpO7wyU1KRQSGpyKCQVGRQSCoyKCQVGRSSigwKSUUGhaQig0JSkUEhqcigkFRkUEgqMigkFRkUkooMCklFBoWkIoNCUpFBIanIoJBUZFBIKjIoJBUZFJKKDApJRQaFpCKDQlKRQSGpqJJbCnba3r17+eIXv8j1118PbAFeDZwMvBOY0tXapO74MbAd2AY8AxwBXFXZ2SfkjOLyyy9n+fLlDAwMAIuBAeBu4BZgb1drk7pjLfAr4CjgVZWffcLNKDZt2sS1117LOeecwze/+U0irmm0HAncATwAvKl7BUpd8VHqIQGwCniu0rNPuBnFrbfeSmbysY99bFTLW4BDgY1dqErqtqPKXdpQSVBExKKI2BwRWyLiiibth0XE1xrtd0fE7PGOtX79eiZNmsTcuXNHtRwKHEP9NZqkKrUdFBFxCPW5znuoLxZcEBEDo7pdCjydmX8BfAH45/GOt23bNqZNm8Zhhx3WpPVwYBewZ7ynl9REFTOKucCWzHwkM58D1gBLRvVZAqxubH8DmB8RMZ7Bdu3aNUZIwJ+WXJ4fz6kljaGKoDgeeGzE/tbGsaZ9MnMP8Hvg9eMZbOrUqTz77LNjtO6bSRw6nlNLGkMVQdFsZpDj6ENELIuI4YgY3rFjR9PBjjvuOJ566qkxwuIPwFQm4B9zpJ5WRVBsBWaO2J/By1cUX+wTEZOpXw3yu9EnyszBzKxlZq2vr6/pYKeffjp79+7lnnvuGdXyPPA4cNy4PglJY6siKNYDcyLihIiYApwPDI3qMwQsbWy/D/hBZr5sRtGK8847j4hg5cqVo1p+Rj0sTh3PaSXtR9tz9MzcExEfAb4LHAJ8OTM3RcQKYDgzh4D/AL4aEVuozyTOH+94p556KpdddhnXXXcd55xzDvV1iaeoX5k5C4NCr0w/p37pNtT/8vcCn/nMZwCYNWsWF110UVtnj3H+Yj/garVaDg8PN2174YUXWLlyJYODgzz00C+pr0vs+1+Psf4iIh3MbgR+3bRl3rx53HXXXU3bIuLezKyVzj4hg2KkP13CLWmkzPI/hbUaFBPuEm5JnWdQSCoyKCQVGRSSigwKSUUGhaQig0JSkUEhqcigkFRkUEgqMigkFRkUkooMCklFBoWkIoNCUpFBIanIoJBUZFBIKjIoJBUZFJKKDApJRQaFpCKDQlKRQSGpyKCQVGRQSCoyKCQVGRSSigwKSUUGhaSitoIiIo6KiO9FxMONxyOb9HlzRPw0IjZFxMaIOK+dMSV1XrsziiuAtZk5B1jb2B9tF3BxZp4MLAJWRsTr2hxXUge1GxRLgNWN7dXA2aM7ZOZDmflwY3sb8CTQ1+a4kjqo3aCYnpnbARqPR++vc0TMBaYAv2xzXEkdNLnUISK+DxzTpOnKP2egiDgW+CqwNDP3jtFnGbAMoL+//885vaQDqBgUmblgrLaIeCIijs3M7Y0geHKMfq8FvgP8U2au289Yg8AgQK1Wy1Jtkjqj3ZceQ8DSxvZS4FujO0TEFOA/ga9k5tfbHE9SF7QbFJ8FFkbEw8DCxj4RUYuIGxp9zgXOBC6JiA2Njze3Oa6kDiq+9NifzPwtML/J8WHgA43tm4Gb2xlHUnd5ZaakIoNCUpFBIanIoJBUZFBIKjIoJBUZFJKKDApJRQaFpCKDQlKRQSGpyKCQVGRQSCoyKCQVGRSSigwKSUUGhaQig0JSkUEhqcigkFRkUEgqMigkFRkUkooMCklFBoWkIoNCUpFBIanIoJBUZFBIKjIoJBW1FRQRcVREfC8iHm48Hrmfvq+NiP+JiOvaGVNS57U7o7gCWJuZc4C1jf2xfBr4YZvjSeqCdoNiCbC6sb0aOLtZp4h4KzAd+O82x5PUBe0GxfTM3A7QeDx6dIeImAR8HvjHNseS1CWTSx0i4vvAMU2armxxjA8Dt2fmYxFRGmsZsAygv7+/xdNLOtCKQZGZC8Zqi4gnIuLYzNweEccCTzbp9nbgjIj4MPAaYEpE7MzMl61nZOYgMAhQq9Wy1U9C0oFVDIqCIWAp8NnG47dGd8jMC/dtR8QlQK1ZSEjqXe2uUXwWWBgRDwMLG/tERC0ibmi3OEm9oa0ZRWb+Fpjf5Pgw8IEmx28CbmpnTEmd55WZkooMCklFBoWkIoNCUpFBIanIoJBUZFBIKmr3yky9oj0FbAR+CTwN7AGOBE4G3gZM6V5pqpRBoTbcB6wH/hJ4E/UJ6qPAD4BN1K+5O7RbxalCBoXaMACcAbxqxLHTgaOAHwM/A/6qC3Wpaq5RqA3H89KQ2OeUxmOzfybWRGRQ6AD438bja7pahapjUKhie6m/Neok4NQu16KqGBSq2H8BW4F3AtO6XIuqYlCoQj8A7gHeSn2RUweLCR8U06e/utslCIA7gR8Bbwb+tsu1qOqfiwn/59HHH/94t0t4xbvmmmu4+uofcvHFF3PjjTcyadKE//2jUfyKqi0rVqzg6quv5qKLLjIkDmITfkah7lm1ahVXXXUV/f39LFiwgFtuueUl7dOnT2fhwoVdqk5VMig0buvXrwfgN7/5DUuXLn1Z+7x58wyKg0Rk9ubtM2q1Wg4PD3e7DOmgFhH3Zmat1M8XlJKKDApJRQaFpCKDQlJRzy5mRsQO4NctdJ1G/a2WelUv19fLtUFv19fLtUHr9c3KzL5Sp54NilZFxHArq7bd0sv19XJt0Nv19XJtUH19vvSQVGRQSCo6GIJisNsFFPRyfb1cG/R2fb1cG1Rc34Rfo5B04B0MMwpJB9iECYqIWBQRmyNiS0Rc0aT9sIj4WqP97oiY3WP1LY+IByNiY0SsjYhZvVLbiH7vi4iMiI6u5rdSX0Sc23j+NkXELc36dKO2iOiPiDsj4r7G13ZxB2v7ckQ8GREPjNEeEfHvjdo3RsRbxj1YZvb8B3AI9dtRvYH67ad+DgyM6vNh4EuN7fOBr/VYfe8Epja2P9Sp+lqprdHvcOpvUbUOqPXYczeH+t2GjmzsH91DtQ0CH2psDwCPdvC5OxN4C/DAGO2LgTuAoH7rtrvHO9ZEmVHMBbZk5iOZ+RywBlgyqs8SYHVj+xvA/IiIXqkvM+/MzF2N3XXAjF6preHTwOeA3R2qa59W6vs7YFVmPg2QmZ26YUgrtSXw2sb2EcC2DtVGZv4I+N1+uiwBvpJ164DXRcSx4xlrogTF8cBjI/a3No417ZOZe4DfA6/vSHWt1TfSpdSTvhOKtUXEacDMzPx2h2oaqZXn7o3AGyPiJxGxLiIW9VBtVwPvj4itwO3A33emtJb8ud+XY5oob1zTbGYw+s81rfQ5UFoeOyLeD9SAeQe0ohFDNjn2Ym0RMQn4AnBJh+oZrZXnbjL1lx/voD4T+3FEnJKZz/RAbRcAN2Xm5yPi7cBXG7XtPcC1taKyn4mJMqPYCswcsT+Dl0/xXuwTEZOpTwP3Ny2rUiv1ERELgCuBszLz2R6p7XDq9wC8KyIepf5adqiDC5qtfm2/lZnPZ+avgM3Ug6MXarsUuA0gM39K/R6LvXJDk5a+L1vSqYWXNhdtJgOPACfwp0Wlk0f1uYyXLmbe1mP1nUZ9YWxOrz13o/rfRWcXM1t57hYBqxvb06hPp1/fI7XdAVzS2D6p8YMYHXz+ZjP2Yubf8NLFzHvGPU6nPqEKnpDFwEONH7YrG8dWUP/tDPUk/zqwhfpdaN7QY/V9H3gC2ND4GOqV2kb17WhQtPjcBfBvwIPA/cD5PVTbAPCTRohsAN7VwdpuBbYDz1OfPVwKfBD44IjnbVWj9vvb+bp6ZaakoomyRiGpiwwKSUUGhaQig0JSkUEhqcigkFRkUEgqMigkFf0/epVocLYI18gAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "import numpy as np\n", "import matplotlib.pyplot as plt\n", "\n", "import pylbm\n", "\n", "L, W = 1., 1.\n", "dico_geom = {\n", " 'box': {'x': [0, L], \n", " 'y': [-.5*W, .5*W], \n", " 'label':list(range(4))\n", " }\n", "}\n", "geom = pylbm.Geometry(dico_geom)\n", "print(geom)\n", "geom.visualize(viewlabel=True);" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### The stencil\n", "\n", "The stencil of the $\\DdQq{2}{9}$ is composed by the nine following velocities in 2D:\n", "$$\n", "\\begin{gathered}v_0=(0,0),\\\\ v_1=(1,0), \\quad v_2=(0,1), \\quad v_3=(-1,0), \\quad v_4=(0,-1), \\\\ v_5=(1,1), \\quad v_6=(-1,1), \\quad v_7=(-1,-1), \\quad v_8=(1,-1).\\end{gathered}\n", "$$" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "+---------------------+\n", "| Stencil information |\n", "+---------------------+\n", " - spatial dimension: 2\n", " - minimal velocity in each direction: [-1 -1]\n", " - maximal velocity in each direction: [1 1]\n", " - information for each elementary stencil:\n", " stencil 0\n", " - number of velocities: 9\n", " - velocities\n", " (0: 0, 0)\n", " (1: 1, 0)\n", " (2: 0, 1)\n", " (3: -1, 0)\n", " (4: 0, -1)\n", " (5: 1, 1)\n", " (6: -1, 1)\n", " (7: -1, -1)\n", " (8: 1, -1)\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAATcAAAE/CAYAAAAjcYRfAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAFR5JREFUeJzt3HuQVvWdoPHnS0M3lwabm6CAgREVBJRWVrzVSLysxLJm2oxWIu4E1ySsM3FWsmyNsbITLzuVqTEmY3bRdUikiPfNxkiMMRoZgxbGC+EqiDFobAEFgthC2wjd8Ns/urXojGU39Guf7l8/n6qu9Hnf857zrVNvHs973j5ESglJyk2vogeQpE+DcZOUJeMmKUvGTVKWjJukLBk3SVkybup2IuKYiKiPiLKW5aUR8ZWi51LXYtx0SCLi7Ij4TUS8FxE7I+LZiPgPLc9dGRHLPu0ZUkpvppQqU0r727N+RHw9Ira2zLwwIio+7RlVPOOmdouIQcCjwP8GhgCjgJuAvUXO9Uki4kLgG8B5wFjgz2ieWZkzbjoUxwOklB5IKe1PKe1JKf0qpbQ2IiYCdwJntHxkrAOIiIqIuDUi3oyIbRFxZ0T0a3luRkRsjoh5EbE9It6OiP/84c4iol9EfDcialvOupa1PDY2IlJE9G7HzLOBu1JK61NK7wL/E7iyxMdFXZBx06F4FdgfET+KiM9FxOAPn0gpbQCuBp5r+chY1fLUP9McxanAeJrP9r510DZHAke0PP5l4PaDtnsrcCpwJs1nin8PHDjEmScBaw5aXgOMiIihh7gddTPGTe2WUtoFnA0k4AfAHyPikYgY8XHrR0QAXwW+nlLamVLaDXwb+OJBqzUCN6eUGlNKjwH1wAkR0Qu4Crg2pbSl5UzxNymlQ/0IXAm8d9Dyh78PPMTtqJtpz2m99JGWM7QrASJiAnAvcBtw+cesPhzoD6xo7hwAAZQdtM47KaWmg5YbaA7SMKAv8FoHR64HBh20/OHvuzu4XXVxnrnpsKWUXgEWAZM/fOhPVtkB7AEmpZSqWn6OSClVtmPzO4APgGM7OOZ64OSDlk8GtqWU3ungdtXFGTe1W0RMaLn4P7pleQzNZ2zPt6yyDRgdEeUAKaUDNH98/ZeIOLLlNaNavsH8RC2vXQh8LyKOjoiyiDjjMP6M427gyxFxYsu1vP9Bc5CVOeOmQ7EbmA68EBHv0xy1dcC8luefovlMaWtE7Gh57DpgI/B8ROwClgAntHN//x14CVgO7KT5y4lDes+mlB4HbgF+DdS2/NxwKNtQ9xT+Y5WScuSZm6QsdThuETEmIn4dERsiYn1EXFuKwSSpIzr8sTQijgKOSimtjIiBwAqgJqX0cikGlKTD0eEzt5TS2ymllS2/7wY20PzX5pJUmJJec4uIsUA18EIptytJh6pkdyhERCXwEDC35TadP31+DjAHoH///qcef/zxpdp1t7d//37KysraXrEHiKZ6Dhw4QJQPanvlHsL3R2urV6/ekVIa3tZ6JflTkIjoQ/M/hfNESul7ba1fXV2dVq1a1eH95qKuro6qqqq2V+wJti2lvr6eymMvLnqSLsP3R2sRsSKlNK2t9UrxbWkAdwEb2hM2SeoMpbjmdhbw18C5EbG65eeiEmxXkg5bh6+5pZSW0fwvPUhSl+EdCpKyZNwkZcm4ScqScZOUJeMmKUvGTVKWjJukLBk3SVkybpKyZNwkZcm4ScqScZOUJeMmKUvGTVKWjJukLBk3SVkybpKyZNwkZcm4ScqScZOUJeMmKUvGTVKWjJukLBk3SVkybpKyZNwkZcm4ScqScZOUJeMmKUvGTVKWjJukLBk3SVkybpKyZNwkZcm4ScqScZOUJeMmKUvGTVKWjJukLBk3SVkybpKyZNwkZcm4ScqScZOUJeMmKUvGTVKWjJukLBk3SVkybpKyZNwkZcm4ScqScZOUJeMmKUvGTVKWjJukLBk3SVkybpKyZNwkZal30QP0ZA8/vIFvf3sZ69Zto7y8N1OmHMnPf345gwf3K3q0TvfVrz7Cs89u4s3anZSXB6ef8R7f+c4FTJp0ZNGjFWbGjEU8/XRtq8cmTRrOunV/W9BE3UtJ4hYRC4GLge0ppcml2GbuHnjgJWbN+ikVFWVcdNGxDBkygOXL36KhobFHxu2HP1zF9OmjuPySI3ny6Z388pcbWbt2Gxs3/lf69u3Z/w2++uqpVFRUAHDUUZUFT9N9lOpdswiYD9xdou1lLaXEddctAeDxx/8TU6dWUVVVVfBUxXr22as488wxsG0p6195h8kz1rFly25efvmPnHLKUUWPV6h/+qdzevz743CU5JpbSukZYGcpttUT/P73O9m0aRf9+vXmllueZfToOxg//n9x++0vFj1aYc48c8xHv+9rTAD06hWeqQBjx97J4MH/zHnn3c3y5VuKHqfb8AuFAuzY0QDAnj1NvP76u9TUHMeWLbu55ppfsnjxKwVPV6z69/dz9d+/AcC8eWdw1FEDix2oQAMHVnDxxcdzySXHccwxR/DUU3/gwgvvZevW+qJH6xY67WJGRMwB5gCMGjWKurq6ztp1l1NR0fTR73fccT7HHTeAvn17c9dda/nJT15ixoyRBU5XnHfe2cMXPr+SFWsbmD17MtdfP61Hv0/uvnsmEcHu3bupqOjPtGk/YtOm3Tz66HouvfSEosfr8jotbimlBcACgOrq6tSTryFMmTKQQYMq2LVrLwMHDmTgwAGUl5cDMGRIZY+8vlJbW8dFFz3Eq6828N/+y0i+e+dfFT1SoRoaGqmr+4Cjj24+c+3ffyC9ejV/0Bo0qGe+Rw5Vz/4aqiDl5WXMnTudm29+hi99aTGnnnokDz30KmVlwRVXTCl6vEKceeZC3nprN8eMquCDvQeYO/dxAGbNmsJpp40qeLrOt337+5xwwnzOPXccI0f2ZeXKP1Jb+x4jRgzg3HPHFT1et1CqPwV5AJgBDIuIzcANKaW7SrHtXP3DP5zDvn37WbRoDQ8//CqTJx/JTTfNYPr00UWPVoi33toNwJtb9nLHou3AdgCmTh3ZI+M2dGg/vvSlk3jqqTdYunQXRxzRl5qaCfzjP36WYcP6Fz1etxAppU7faXV1dVq1alWn77erqqur82PGh7Ytpb6+nspjLy56ki7D90drEbEipTStrfX8tlRSloybpCwZN0lZMm6SsmTcJGXJuEnKknGTlCXjJilLxk1SloybpCwZN0lZMm6SsmTcJGXJuEnKknGTlCXjJilLxk1SloybpCwZN0lZMm6SsmTcJGXJuEnKknGTlCXjJilLxk1SloybpCwZN0lZMm6SsmTcJGXJuEnKknGTlCXjJilLxk1SloybpCwZN0lZMm6SsmTcJGXJuEnKknGTlCXjJilLxk1SloybpCwZN0lZMm6SsmTcJGXJuEnKknGTlCXjJilLxk1SloybpCwZN0lZMm6SsmTcJGXJuEnKknGTlCXjJilLxk1SloybpCwZtwLNnr2YUaO+x4gR8xk27BZmzryXVaveLnqsQnzwQRN/93ePceSk3zB80krOOmshL7ywueixCnPbbc9z0kn/h7Kymxk8+PvceOPSokfqdkoSt4iYGRG/i4iNEfGNUmyzJ6itreOccz7DFVecyNCh/Xniideoqfm/RY9ViLlzH2f+/OWMGF7OxRdU8dxzm7jggnvYsaOh6NEKsWLF2wwZ0o8xYwYVPUq31bujG4iIMuB24AJgM7A8Ih5JKb3c0W3nbunSKwGoq6vj9df3cOqpC9i8eReNjfvp06es2OE60fbt77Nw4Sp69Qr+7f+dRP9+++hb1ci9965l/vwXufHGGUWP2OnuuecSAGpqHqS29r2Cp+meSnHmdhqwMaX0ekppH/Ag8Jcl2G6PMH/+i8yb9xSXX/4QAPPmndGjwgawfv12GhsPcMwxR3Dk8HIApk07CoDVq7cWOZq6sQ6fuQGjgE0HLW8Gpn/SC6KpHrYtLcGuu7+f3L+ap59r/i/z6KMrOGvSez3u2Gx7dTsAlX0b4d3VlO3Zw4D9YwHYuumtHnc8Wtm7o/l/69/o2cfhMJQibvExj6V/t1LEHGAOwLjRw6ivry/Brru/R+8Zz8536/nNb5u44muvcelX1rN6yWQ+M7qi6NE6zaAB+wHYXd/Inj2NfPDBB7yzs/la27AhvXr0e6WpqQmAffv29ejjcDhKEbfNwJiDlkcDb/3pSimlBcACgOrq6lR57MUl2HX3tWdPI+XlZZSV9aKpro5LplRSed2t7Nq1l+0HTmHSseOKHrHTnFpZT58+/8KmtxrZ1ed0Bgzax9o/7AM2c+qZJ1N57GeLHrEwvQc8CLxH+ZDjqTx2RtHjdCuliNty4LiIGAdsAb4IzCrBdrP2wgtbmDXrIf78zz9D//69ePHFrezatZfhw/tzyilHFT1epxoxopIrr5zKD36wkvMuXcOE8RX89LF3qaws55prTit6vEL88IcrWbbsTVaubP7ToMWLX+GNN+qoqZlATc2EgqfrHjoct5RSU0RcAzwBlAELU0rrOzxZ5o4+eiDHHz+UJ598nd279zJ8+AAuu+xEvvWtczjiiL5Fj9fpvv/9mfTp04sfP7iGjX/Yw+mnj+a73/2PDB8+oOjRCrFs2Zv86EdrPlpes2Yba9ZsY+zYKuPWTpHSv7s89qmrrq5Oq1at6vT9dlV1dXVUVVUVPUbXsG0p9fX19PTLFgfz/dFaRKxIKU1raz3vUJCUJeMmKUvGTVKWjJukLBk3SVkybpKyZNwkZcm4ScqScZOUJeMmKUvGTVKWjJukLBk3SVkybpKyZNwkZcm4ScqScZOUJeMmKUvGTVKWjJukLBk3SVkybpKyZNwkZcm4ScqScZOUJeMmKUvGTVKWjJukLBk3SVkybpKyZNwkZcm4ScqScZOUJeMmKUvGTVKWjJukLBk3SVkybpKyZNwkZcm4ScqScZOUJeMmKUvGTVKWjJukLBk3SVkybpKyZNwkZcm4ScqScZOUJeMmKUvGTVKWjJukLBk3SVkybpKyZNwkZcm4ScqScZOUJeMmKUvGrQBLl75BxE0f/Qwe/P2Pfl+0aHXR4xXqgYe3M3D8CiJuYu7cx4sepzCrV2/lwgvvZejQWzj66Ns58cTbueOO5UWP1a307siLI+Iy4EZgInBaSum3pRgqd6NHD+Laa6d/tLxzZz333LMegPHjhxQ1VuE2b97F337j9/TuDU1NRU9TrJqaB6mtfY8pU45k3LhB/PznG/na1x5j4sRhfPaz44oer1voUNyAdcDngX8twSw9xvjxQ7jttpkfLd9yy1IAqqtHcvbZxxQ0VbFSSsyevZijR5Qz8biBPPSLd4seqTCNjfvZtGkXAPff/1eMHl3O+ef/mBUr3uaNN+oKnq776NDH0pTShpTS70o1TE+UUmLBgjUAfP3rpxc8TXFuu+15li17k/vumEhFRc++WtKnT9lHZ/ZXXPFTZs/+BStXvs3JJ4/gkksmFjxd99Gz30VdwKOPvsprr9UxcmQlX/jC5KLHKcS6ddu5/vp/4+abZzB1cmXR43QJNTUTGDu2irVrt/HIIxvp3bsXNTUTGDiwvOjRuo02P5ZGxBJg5Mc89c2U0s/au6OImAPMARg1ahR1dZ5eA9x66zIArrpqMg0Nu2loKHigAtx330r27dvPkiUb+fUTdbz08vsALF68gYj93HDDWQVP2Ll27tzD5z53Lw0NTTz22KWMHl3BlVcu4aabnqayMvjKV04uesRuoc24pZTOL8WOUkoLgAUA1dXVqaqqqhSb7dZeemkbzzyzmb59y5g792yqqgYUPVIhKir6khIsWVLb6vHa2l2sWvVHetp7ZePGBhoamujTpxfnnnsCe/bUM3nySFau3EZt7fs97ngcLj+WFui2254H4LLLJjB8eM8MG8CNN84gpRuaf7aew6zPDwXg2muns3TplcUOV4CJE4cxZEg/GhsPcN55d/M3f/MrHnjgJYAe+4XT4ehQ3CLikojYDJwB/CIinijNWPnbsaOB++9fB8DVV08teBp1JQMGlPPYY7M4//w/Y8OGHfzsZ79v+Yb9wh57XfZwREqp03daXV2dVq1a1en77arq6ur8qPGhbUupr6+n8tiLi56ky/D90VpErEgpTWtrPT+WSsqScZOUJeMmKUvGTVKWjJukLBk3SVkybpKyZNwkZcm4ScqScZOUJeMmKUvGTVKWjJukLBk3SVkybpKyZNwkZcm4ScqScZOUJeMmKUvGTVKWjJukLBk3SVkybpKyZNwkZcm4ScqScZOUJeMmKUvGTVKWjJukLBk3SVkybpKyZNwkZcm4ScqScZOUJeMmKUvGTVKWjJukLBk3SVkybpKyZNwkZcm4ScqScZOUJeMmKUvGTVKWjJukLBk3SVkybpKyZNwkZcm4ScqScZOUJeMmKUvGTVKWjJukLBk3SVkybpKyZNwkZcm4ScqScZOUJeMmKUsdiltEfCciXomItRHxcERUlWowSeqIjp65PQlMTimdBLwKXN/xkSSp4zoUt5TSr1JKTS2LzwOjOz6SJHVcKa+5XQX8soTbk6TD1rutFSJiCTDyY576ZkrpZy3rfBNoAu77hO3MAea0LO6NiHWHPm62hgE7ih6iC/F4tObxaO2E9qwUKaUO7SUiZgNXA+ellBra+ZrfppSmdWjHGfF4tObxaM3j0Vp7j0ebZ25t7GQmcB1wTnvDJkmdoaPX3OYDA4EnI2J1RNxZgpkkqcM6dOaWUhp/mC9d0JH9Zsjj0ZrHozWPR2vtOh4dvuYmSV2Rt19JylJhcfPWrdYi4rKIWB8RByKiR34zFhEzI+J3EbExIr5R9DxFi4iFEbHdP5uCiBgTEb+OiA0t/z+5tq3XFHnm5q1bra0DPg88U/QgRYiIMuB24HPAicDlEXFisVMVbhEws+ghuogmYF5KaSJwOvC1tt4fhcXNW7daSyltSCn9rug5CnQasDGl9HpKaR/wIPCXBc9UqJTSM8DOoufoClJKb6eUVrb8vhvYAIz6pNd0lWtu3rqlUcCmg5Y308abVz1TRIwFqoEXPmm9Dv0pSDuGKMmtW7loz/HoweJjHvOrfLUSEZXAQ8DclNKuT1r3U41bSun8T3q+5dati2m+dSv7N3Jbx6OH2wyMOWh5NPBWQbOoC4qIPjSH7b6U0k/bWr/Ib0s/vHXrL7x1S8By4LiIGBcR5cAXgUcKnkldREQEcBewIaX0vfa8pshrbt66dZCIuCQiNgNnAL+IiCeKnqkztXy5dA3wBM0Xi3+cUlpf7FTFiogHgOeAEyJic0R8ueiZCnQW8NfAuS29WB0RF33SC7xDQVKWusq3pZJUUsZNUpaMm6QsGTdJWTJukrJk3CRlybhJypJxk5Sl/w8zJYT2fOmD4AAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "dico_sten = {\n", " 'dim': 2,\n", " 'schemes': [\n", " {'velocities': list(range(9))}\n", " ],\n", "}\n", "sten = pylbm.Stencil(dico_sten)\n", "print(sten)\n", "sten.visualize();" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### The domain\n", "\n", "In order to build the domain of the simulation, the dictionary should contain the space step $\\dx$ and the stencils of the velocities (one for each scheme). " ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "+--------------------+\n", "| Domain information |\n", "+--------------------+\n", " - spatial dimension: 2\n", " - space step: 0.1\n", " - with halo:\n", " bounds of the box: [-0.05 -0.55] x [1.05 0.55]\n", " number of points: [12, 12]\n", " - without halo:\n", " bounds of the box: [ 0.05 -0.45] x [0.95 0.45]\n", " number of points: [10, 10]\n", " \n", " +----------------------+\n", " | Geometry information |\n", " +----------------------+\n", " - spatial dimension: 2\n", " - bounds of the box: [0. 1.] x [-0.5 0.5]\n", " - labels: [0, 1, 2, 3]\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAARIAAAEICAYAAACTenveAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJztnXl4VPd9r9/vzEhIQkISWpBYxCZ2zGIDxlu84S2L7RjbiQ3YTtPmpmnTtGn6NHmS9vFNbm7T5PamN83W2LGxwUmc2Bho6sQsBttAwGYTiwCzgyRAu4SE1pnf/eMcmdEwmjmj0WhG0vd9Hh7Ncn7z+WnQeTXn6JzzEWMMiqIo0eCK9wQURRn4qEgURYkaFYmiKFGjIlEUJWpUJIqiRI2KRFGUqFGRKHFDRIpEpElE3PGeixIdKpIhiIicEZEWEbksIvUiskNEvigi/frzYIw5Z4xJN8Z4+zNX6XtUJEOXTxljMoDxwPeAfwR+Gd8pKQMVFckQxxjTYIxZD3wGeFpEZotIpoi8LCJVInJWRL7V9WlFRJ4Rke0i8kP708wpEbnZfvy8iFSKyNNdry8inxCRfSLSaD//rN9zE0TEiIjHvr9VRL5jv/5lEdkgIrn9/JYovUBFogBgjHkfKANuA/4DyAQmAbcDTwGf81v8RuAAkAP8CvgNsBAoBpYDPxaRdHvZZnt8FvAJ4C9F5OEQU3nSzsoHkoGv9cG3p8QYFYniTwUwEuvTyTeMMZeNMWeAfwNW+C132hjzor1v41VgHPBtY0ybMWYD0I4lFYwxW40xB40xPmPMAeDXWHLqiReNMR8aY1qA3wLz+vh7VGKAikTxZwzgwfokcNbv8bP2c11c8rvdAmCMCXwsHUBEbhSRLfZmUgPwRSDU5spFv9tXul5HSWxUJAoAIrIQSxZrgQ6snbBdFAHlvXzpXwHrgXHGmEzg54BEMVUlAVGRDHFEZISIfBJrP8dqY0wJ1ibFd0UkQ0TGA18FVvcyIgOoNca0isgirH0gyiDDE+8JKHHjv0SkE/ABpcD/xfq0APBlrB2up4BW4DnghV7mfAn4NxH5MfAOlqSyopi3koCIXthIUZRo0U0bRVGiRkWiKErUqEgURYkaFYmiKFGTsH+1yc3NNRMmTIj3NBRlULNnz55qY0xetK+TsCKZMGECu3fvjvc0FGVQIyJnwy8VHt20URQlalQkiqJEjYpEUZSoUZEoihI1KhJFUaJGRaIoStSoSBRFiRoViaIoUaMiURQlalQkiqJEjYpEUZSoUZEoihI1KhJFUaJGRaIoStSoSBRFiRoViaIoUaMiURQlavpEJCJyv4gcE5ETIvL1HpZ5XERKReSwiPyqL3IVRUkMor7Uooi4gZ8A9wBlwAcist4YU+q3zBTgG8Atxpg6EcmPNldRlMShLz6RLAJOGGNOGWPasTpkHwpY5i+Anxhj6gCMMZV9kKsoSoLQFyIZA5z3u19mP+bPVGCqiGwXkZ0icn+wFxKRL4jIbhHZXVVV1QdTUxSlP+gLkUiQxwILhT3AFOAO4AngeRG5pkjaGPMLY8wCY8yCvLyor5CvKEo/0RciKQPG+d0fC1QEWWadMabDGHMaOIYlFkVRBgF9IZIPgCkiMlFEkoHPAusDllkL3AkgIrlYmzqn+iBbUZQEIGqRGGM6gb8G3gKOAL81xhwWkW+LyIP2Ym8BNSJSCmwB/sEYUxNttqIoiYEYE7g7IzFYsGCB6bemvX1fh9aLV++nFMD878V2rGZqZl+MjRIR2WOMWRDt6yRsZWe/0noRhk+4er/5TOzHaqZm9sXYBEEPkVcUJWpUJIqiRI1u2oC1Ter/cTKlIPZjNVMz+2JsgqA7WxVlCNNXO1t100ZRlKhRkSiKEjUqEgcYA2c3QWtd5GNPXqjhD9tORDzO6/OxdtMxLtVfjjyzopY/bj/Zq8w3NvYu83h5TVSZVQ1NEY/9sKyat3b0f+aGP0We6cPHBt6nkeaIxw4EVCQOqNgOFdugdGVwmVzgWS7w7DWPn2+rZs2qs+zfdJl3twaeftTzOB8+frf9IEe2XeHVladpqfc5z2ytZs3qM+zb2Mh7Wy9ElPnbbQc5uv0Kv115JqLMc61VrF19ln0bG9n2TmSZr753gKPbr/Cbladobbh2f11PY8+2VLL2lXPs3dDI9vciy/zNu1bmqy+djjxz9Xn2vNXIjm0XHY/z4WMN77KDQ/yaTZhrzmkd+KhIHDBqIWSMg7aGnmUSSAXVvDLsj7R8/Ax5kol7awHnt4Yf58PHerZzdOF+ZGwLU+oncWSli9Z6h5kpf6TlgbPkSSauraMcZ65jG8cW7cc1ppXi+omOM8up4lcpb32UKVtGUfaOs8y1vMeHN5bgGt3GlLpJlK4U2hrCjy2jkl+lbqD1/nPkSxayuYCy95xlvsG7HLczi2sncuSlSDPPki9ZsGkU5ducZa7hXQ5ximEk83EWI0FPmB/YqEgc4BkG05c7l0kF1bzMH2mlnWmzM3lw6Xhc4qJsKyFX7C6J7Oc4SSnC48snMXZsOm31dmaIFbtb5nUj+NQjzjPXsY0STpCc4uKxFRMZO8ZZZjlVrOItWmln+pxMPvVpK/P8FkLKpEsiBzhJcoqLx5+ayOjR6bTVweGVhFyxy6hkFRtoo50ZczP55MPjERHObyakTLokcpBTJKe6+MwKK7O1Fkpfcp45c14Wn3zIyjy3iZAyCZTIcu5hHKN6HjCAUZE4xKlM/FfoGYznUe4gf7abKUsBoccVu5tE8LCce5mUUsCM5ZA+lpArdrDMUde5mfJI+MyPJEISy7jHylwB6WNCZ/pLZCYTWMrtjJrjpvjTVmZPMukmEZJYzr1MTBnFzKdg+GhCysR/hZ7FRCtzrovih+3MHmTSTSIksYJ7mZA6ipkrrMxQMvHPnM0kHuFjjJrnYvJDVmZPMhlKEgEVSUSEk0mwFdqNG4Dc2fQok2ASGY91UJInhZAyCZl5HT3KJJhEumWGkEkwiXRl5s2hR5kEk0iRvXJ5Uggpk2AScdk/vnlz6VEmwSTStUJ7Ugkpk2AS6crMn0ePMhlqEgEVScQEk0l7XQoXaetxhe4iUCaXtk7Ah+lRIh9lBpFJe73DzACZdGX2JJFumQEyaa9P4QJtPUqki0CZVL4zHh+mR4n4ZwbKpL1hGBW09iiRjzIDZFL5XhE+TI8S+SgziEwCMwMl0kWgTKq2WZlDTSKgR7Za9OI07s42OLoaLp+H2uQX2fHYFjrG3MINacVBV2h/qg/B8dehpn0zJ27eS8WdI8hMSw0qkW6ZrXBkNTSVQcXwdWx78vckZS9wlnkQjq/pnpmVlhZUItdkroKm8u6ZC9KmBJWIP1UH4MQbVubxW/Zw4Y5MstLSgkokMLP0ZWiugPL0tWx/8vckZS1kYdrUoBLpllkCJ9bambd9wIWPZZOVlhZUIt0yW6B0lZVZkbGWbU9YmYvSpgWViD+V++HkOqjp2Mil+9/k2PXTyU4b7lwiehmBQYKT07ifbe921wNM74Rdh0+QkrGdW76fzNaHXdxx4yLcuT2vXGB9Mqmra2HrmnrqthZRe8nHEw8vZvzI0OdYdH0yefvndewua8P1wk20LU7mjvkOMq+zM9+wMyu9PPHQTc4yV8Dmn9Wxu9zKbL8xmduvXxg2M28O1Ne1sGVtPfVbxlNb6eXJB2+maGTolavrk8mmn9ayu6Id1ws3075oGLdfvxBXbugP0Xlzre9zy7p66jdPpOail2UP3sK4cJn2J5NNP6vlAzuzY9EwPnb9grCZ+fOgru4Kx05tIqe5nOaNeTSO+CS/dAXPfPb2gAf0MgJDG49UMKvg16QMa8blzmZ+4wLqa9ocjfXlN5O9ZD+eVB9zzUTSqjKcZaZA3t0teEbVkTuqhRtr5jjPHNVM9t12pm9SRJn5S1rw5NeTW9DCotoIMguaGXlXCUlpXub5JpNale44c9Q9rSTlNZBXcIVFNddRV9PqaKwpbCbbzpzvnUxK1XBnman295nXQF5hMwsjyGT0FeqLmsFtuDH/ApPliLNxgwT9ROKUZ5O736+vgJ2/IbnBR2fFXOoK7gM85OSkOXq5nJw0Mvc+zvRpnSS5nI8DGFeYRXHhIpI8bsRI/2WOXti7zKbHSJvay8wxC0hyuxEiy8xqeozhvcgsGp3NlNELrO8zwsxRe/8HknmYWaPe5cG5W8hIS8HNrY6zBzK6jwQi30atr4CdL0NHKxTMoHr8x6mpayMnJ43cXOc/tNXVV6ipuRLxuGjGambsM3Pzj5ORvQGDwcOS8DIZBPtIVCSREiARbngUXKH3FShDDy/78bLOuUziREJdRsBJibi93KMiYkQk6onHBZWI4hA383DzEILQySa8ODiefgATtUj8SsQfAGYCT4jIzCDLZQB/A+yKNjMetJdVYP6kElGcEyiTy2yK95RiRn+ViAN8B/g+4HA3eOLQfKyCSz95mcbjrRiViBIBXTJppJn3+Sn7WBnvKcWEfikRF5H5wDhjzO9DvVCiloh3VNZDRzuXymdQ1vgoRlQiinMqKGQtWXTipZmz+Lj2Eg0Dnb7482/IEnERcQE/BJ4J90LGmF8AvwBrZ2sfzK1PyLptJj7305RvHIt5zw0uGHsHyOA7G1zpY85zidVspI1CRvMEN/FoyKNkByr9USKeAcwGtorIGWAxsH6g7XAdefN4ipe6EZd1IlrZVuvKaYrSE1cl0s51TOIBHsc9SA/dinmJuDGmwRiTa4yZYIyZAOwEHjTGJODfdkOTOxuKH0FlooQlUCKfDnO+zkCnv0rEBw0qEyUcQ00i0EeHyBtj3gTeDHjsn3tY9o6+yOxTIjyyMHe29fXEGksmDfUteKY3kZs3fMAcfamZsck81HCWTQV/wjMc5xKJ45GtfcXg3GCLlF6c/ZsLUCkcOeCledMbpKSnsmHFZO69r9jRD1919RU2bjyJyyX4fIZ77pns+Ie2t2M1M7aZf/v2Dibc8CJ/PDqXMXID5vLHOBhEInr2r9KN3Dwvkye+QWbGQdLc+0lutn6TOaGm5goulzBqVDoulzgeF81YzYzd2OqaJibMeoP84fU8lnqc4vPTkSG0euknEqcEnv3r88GB9WQdO8Sl/a2czbuHjozhEZ0t6vMZLl1qwuczEZ2h2tuxmhm7sbk56Ux750nSZq1ivBcm3/1rstL/AiHTcfZARk/ag8i3UW2JcH4/uJOonbKUKl/ugNmO18wYZtbWkF+4jrSMaoSRJPF0eJno2b+xI2HP/g2QCDcuh5zx8Z6VkkAYWuhkFT4qnMskTiTU2b9DBpWI4gAhFQ8rcDEaQy0dvITBQQvXAEZF4hSViBIBQ00mKhIHmE4flS+vp/NDlYjiHH+Z+KjhIN+mgfJ4TysmqEgcUL1+B20H91NZkkTrdJWI4hwhFRfL2E8TFznBbr47KM/+VZE4IPvehUj+ZC7IckrXjXdUIq4o0NW69z6bmEIzBUznC4PycPnB9x3FAE/6MAq+spxhE8c7KhFXFLAk8jrvcJjTuElnId9mDPPiPa2YoCJxiCdFHJWIKwp0l8gwklnBvYwdxNWdKpIICFcirijQk0Ty4z2tmKIiiRCViRKKoSgR0HNtLCI8RLlLJl0l4vt/dIYRd3nJn1E4MA7j1syYjK2sbuRN7ypKMzxkpw13LhG9jMAgIYoS8SP7ynF1vAQ7ktm89G7ufuS6hD7VXTNjM7aquokfnPsuo0cdoPzUTC5X/T3Pu4JLRC8joHTDIxVMGPcSw4a14EoqhPSUhD7VXTNjN7a2ppW6yutwSxJLOluZLIcdZw4G9BOJU3ooEXc3dHCxPI/yguvxiSuhT3XXzNhm3r13LmmHPYwteJfJi7baJeK3Oc4eyOjZv6Al4prZp5m5+SfIyH7L7v29O7xM9DICsSNhLyOg/b+KA7yU4GWtc5nEiYS6jEC4EnER+aqIlIrIARHZLCID82QVlYjiEDdzcfOw3fu7GS/vxXtKMaW/SsT3AQuMMXOA17A6gAcUbee1RFyJjECZNLIh3lOKGf1SIm6M2WKM6dr9vROrjW/A0HysgkotEVd6QZdMGmjmA37GXl6M95RiQr+UiAfweeAPwZ5I6BLxTqtE/HyDlogrkVFOgV0i7uMKZYPyMgIxLxHvtqDIcmABEHhIjjUokUvEk56h/K0xmG1Wifi4O7VEXAnPOS6xmg20U8gYLREPSbgScQBEZAnwTaze37Y+yO1XRi4uovhRq0S8/F04v0WrOpXQXJVIh5aIOyBkiTiAiMwH/hNLIpV9kBkXcmdB8VJUJkpYAiUy2Pt/+6tE/AdAOvA7EdkvIut7eLmER2WihGOoSQT6qUTcGLOkL3JiRqQl4rOsrydet2TSUN9C0gwtEdfMKxxsOMOmgp0kaYn4ECTKEvErm9aQmp7KhhXFWiI+hDP/bst2JtzwAm8dna8l4opzcvO8FE9cQ2bGIVLcB7REfAhnVtc0MX7mWvLSGnlcS8SVHglWIl6yjswPD3Nxf5uWiA/xzNycdKZtWU7adS/ZJeK/Iiv9Cwlb1dnX6El70LsS8ZJ1UFYC7mRqpzyiJeKaaZeI15JfuNYuEc8miWe0RDyeJOzZvwES4cZlWpildMPQSicv2yXiDmUSJxLq7N8hg0pEcYCQgoen7N7fOjpYOah7f0FF4hyViBIBQ00mKhIHmE4flS+to/O4SkRxjr9MrBLx/6kl4kOZ6vU7aDtUwqWSZFqnqUQU5wgpCMvYRxMXOakl4kOZ7HsXwqhiLsoyDq8bT2ttvGekDBS8eFnDLjYxjWYKmcEXB+Xh8oPvO4oBnvRhFH5lOcMmjae9EQ6vRGWihMWLl9d5h1LO4GE4i/gOo5kT72nFBBWJQzzDYPoyyChCZaKExV8iw0jmKe5jDHnxnlbMUJFEgMpEccJQkwioSCJGZaKEYihKBPRcG4velIgvg6OvwOVzsP9Hp60S8ZmjB8Zh3JoZk7GXquv5g281pelJZKcNdy4RvYzAICGaEvH9ZbjaX4Y/JbF56RItER+imVXVTfyf8/+b0fkHuXB6BpcvfY3n3MElopcRULrhkXImjH3ZKhFPHqMl4kM4s7amlbqLc3FLMnd1tFHsPug4czCgn0icEngZgbpy2PWqX4n4fC0RH+KZd++dQ9ohD2ML32HywnftEvGghQmDDj37FyLfRq0rh12rrNa9wplUFz2gJeKa6VcifpKM7D/avb93hpeJXkYgdiTsZQQCJML1S7V1T7kGLwfw8oZzmcSJhLqMgIMS8WEi8qr9/C4RmdAXuf2OSkRxiJs5uPm03fu7BS/vxHtKMaW/SsQ/D9QZY4qBHwL/Gm1uf9N2rhyzUyWiOCdQJg38Md5Tihn9UiJu33/Jvv0acLfIwCm8bD5STuVPV9H4YSumQCWiOKdLJg00s5v/ZA/Px3tKMaG/SsQ/WsYu1GoAcgJfKGFLxGsbrRLxipmcr1+qJeJKRJSRzxt2iXgLFwflZQT6q0TcUdF4wpaI3zID47ZLxLfbJeJ3aYm4Ep6zXOQVNtJOIUUs4yYe0csI9ICTEvGPlhERD5AJDKgzVLIXFzHlMbtE/D04/7ZWdSqhuSqRDuZSzH08qiXiIQhbIm7ff9q+/SjwtknUvzuHIGcmTHkUlYkSlkCJPMStg/KTSBf9VSL+SyBHRE4AXwWu+RPxQEFlooRjqEkE9IA0i14cWVhTCsdfA+OD9OtaSJqpJeKaeYUDDafZXLCLpOE4l8ggOLJ1cG6wRUovzv7NAaiC0hIfVza9Tkp6mpaID/HMv9u6nYk3vMCGo/MYI4swl2+lREvElXDk5PoonvQ6mRmHSdUS8SGdWV3TRNGMteSmNvJY6nEml03REnElCMFKxPevJfN4KRdL2jiTd6+WiA/hzNycdKZvWU7adasY7/Ux+a6uEvEsx9kDGd1HAr0rEd+/FsoPWCXiU5dS5c0ZMNvxmhnDzNpa8gvXkZZRhZBl9/6Gkckg2EeiIomUAImweDmMLIr3rJQEwioRX4WPcucyiRMJdfbvkEElojjAqupcgYsxGOrt3t/6eE8rpqhInKISUSJgqMlEReIA0+mjcuVaOk+oRBTn+MvERy0HeJb6bue3Dh5UJA6oXr+DtsMHrBLx6SoRxTlWifiT7KWJS5xmD/8yKM/+VZE4YOT9i5CCqVx0Lefw2iItxFIcYxVm7WQz02hiNDP50qA8XH7wfUcxwJ2WTOFXniRlUpG26ymO8eLlNbZyhLMkMZzFfIdCZsd7WjFBReIQd7LVrjdivFZ1KuHxl0gKyTzF/YwmN97TihkqkghQmShOGGoSARVJxKhMlFAMRYmAnmtjEeEhyl0yOfoKNJ7VEnHNtLBKxF+mNH0YI9OGO5eIlogPEnpxGQE3MN0LpfvO2yXiyVoiPoQzrRLx7zI6/xAXT82gqfIf+IU7uET0MgJKN9yU2SXirbiSxmLSUxP6VHfNjN3Y2ppWai9cj4dk7uxop9hd4jhzMKCfSJxyTYl4Gex8FXdjJxcq8ikfNR8jktCnumtmbDOX7J1N2qFljB29lckL3yM9LQUPdzjOHsjo2b/QixLxMti5CjrbYPQsqsc9QE1d64DYjtfM2Gfm5p8iI/sPGAxu7ggvE72MQOxI2MsIBEiE+UvBpVuISne8HMTLGucyiRMJcRkBERkpIhtF5Lj9NTvIMvNE5E8iclhEDojIZ6LJjCsqEcUhbq7DzSMIgpetdLI13lOKKdGuBV8HNhtjpgCbCV4zcQV4yhgzC7gf+HcRScyrvISg9UyZVSKuElEcEiiTBt6M95RiRrRrgn85+EvAw4ELGGM+NMYct29XAJVAXpS5/UpTaRlVP1tFw4dtmEKViOKcLpk00MwHPMcenov3lGJCtGvDKGPMBQD7a36ohUVkEZAMnOzh+YQsEe+sbwJvB5cqZnGubilGVCKKc8rIYw1ZePHRStWgvIxA2D//isgmoCDIU9+MJEhECoFVwNPGmKDvZMKWiN88HeP5HOV/GIPZ7gKBoru1RFwJz1kuspoNdFDIeFawmIcH5WUEworEGLOkp+dE5JKIFBpjLtiiqOxhuRHAfwPfMsbs7PVs40j2onFMyYDjv4OKbdZjKhMlFFcl0sk8pnAftwxKiUD0mzb+5eBPA+sCF7CLxd8AXjbG/C7KvLiSMwOmPGb1/lZsg3ObtfdXCU6gRB4cxBKB6EXyPeAeETkO3GPfR0QWiMjz9jKPAx8DnhGR/fa/eVHmxg2ViRKOoSYR0APSLHpTIn7E2swxPkiffYWkWc1aIq6ZHGg4xaaCXSQPF+cSGQRHtuq5NhB1iXjL5tfxDU9jw4opWiI+hDP/bus2Jt3wAhuPzmcsizBNt7BfS8SVcOTk+pgy6TVGpJeS6j6oJeJDONMqEV9HTuplHks9waRyLRFXghG0RPwNRhw/woWSNs7m3acl4kM486MS8dmrGe/1aol4opDQ+0hsiVB+0C4Rf5Qq78gBsx2vmTHMrK0lv3A9aRmVWiKeCCTs2b8BEmHxChg5Lt6zUhIIq0R8NT7KtERcCYJKRHGAVdW5HBdjh0TvL6hInKMSUSJgqMlEReIA0+mj8sU36DiuElGc4y8TLRFXqFq3g7bSg1SWJNM6TSWiOMcqEX+CvTRrifhQJ+eBRVAwjYvuFRxeO46WmnjPSBkoXC0Rn6ol4kMdd1oyo7/yBCmTxtF+GUpXojJRwqIl4so1dKvqVJkoYRhq1Z0qkghQmShOGGoSARVJxKhMlFAMRYmAnmtjEWWJeMmPTjHiTi/5s8cMjMO4NTMmY60S8ZcoTU/REvEhSV+UiO9M4u1HlnDX0jkJfaq7ZsZmrJaIK72me4l4Eb6MtIQ+1V0zYzfWKhG/QUvElTCELRGfpyXiQzxzyd5ZpB1aztjRW7REPFFI6MsIdKvunE31uPu1RFwz/UrET5OR/aaWiCcCCXsZgQCJMP8Rbd1TrsHLIby8riXiDicRtkTcb9kRIlIuIj+OJjOuqEQUh7iZjZulWiLuECcl4l18B3gnyry40b1EXCWihCdQJloi3jNhS8QBROQGYBSwIcq8uPBRifixNkyhSkRxTpdM6mnSEvEQhC0RFxEX8G/AP4R7sYQtEW9otkrEL8zmXO0jWiKuRMR5cniDkVoiTnQl4l8C3jTGnJcwRbkJWyJ+0zSM588of3M0ZoddIr5Ee3+V8JzhAq+wkQ4KmMAKbuKRQXkZgf4oEb8JuE1EvgSkA8ki0mSMCbU/JeHIXjiWqRnw4W+hYrv1mMpECcVViXQynyncx60Ig/MHJuYl4saYZcaYImPMBOBrWGXiA0oiXYycDlMft3t/t8O5Tdr7qwQnUCIPDmKJQP+UiA8qVCZKOIaaREAPSLPoxZGFtUetzRzjg+GzmkiefYXcvPTIj4SsbiQnd0TkR19WNlJT20pOrmYmUmZJ/Qk2F35A8nBxLpFBcGSrnmsDvTr7dyQwtRoOl3RC1tO01/vYkPYs9943xfEZqnvOLKfYlLCj5P9x8913OT9DtbKRi6dvo4009u79OffcG0nmMorNQXbs/xE3L7kjgsx6Lp6+nVbpRebpZRTT28yP0Uo6e/f+LKLMvWeeZLI5xLb9/8GtS26PKPPS6dtokRHs3fvTiDJfPrmC2cM+YNeBzzEi+U5M063sCyIRPftX6cbIHB9Tr3+O7HEnyM47TWHWO9TUNDsa23DmOMXNB0hLbmPK2JXU1ZU5C/X5aC/5TzI91RSklVGYv9V55ukP7cxWpox7kdq6coeZXtpLniPTU01hahmF+e9GkHmMKVd6l9lR8nMyPTUUpp2nIP+9CDKPUtx8kLTkVqaNe5HauoqIMkd4ailMPUdB/jbnmadKuaFzD1nuVr7g/oDJ5cWDfnPGH/1E4pRgJeL71pBxopJTjal4bzKM6DhCXsFODJ8M/UNUe56Cs2s53TCCtnkptIzIZdLo1zF8ASGn53E+H+x7nazaes5czqFzoYsMjpJXsAvDJ8JknqPg7DpON2bSNq+dlhE5TCrsyhwZItML+9YEZJaSV/A+ho+Hzqw5S8HZ9ZxuzKQ10sw9r5FZe5kzzTl03uBiBIfJK/gAwwPhM8+s5/TlTFrnd9AyYqRfZo9ncHTLPNucQ8e7+5GoAAALj0lEQVQNLkZwiLyCSRjuD51ZfYaCs7+nqSmbpoUpFKcWce9dr9gl4iEyBxEqErC2Sf0/TqYEO2zGD1siVBwibUQ6nbcupDXtCtOSc0kbvgcvKbhZEvyHr/Y87FrN8GTDsBkTqRg7njmpU0kbfpEOVto9sUFkYkuEisOkjUinbe7ttKVeYdrIXNKG77Yz7+4h8xzsXM3wYQGZ6f6ZQVZsWyJXM++gLbXJzvzAzrwreGbNWdj1ipU5cyIVY4vszEvhM/e8BhePkJaZYWWmdGW+b2feGSJzNcNTILloEhfGjAuSGWTFDshsnXsHbSmX7cxdduYdwTOrz8D7rzA8BVKLiqkfNZZJKQ4y/Yn05y8BUZFAZDu2/CSCZxjcuBzPyF+STg6ZfJZOfosX60CTa2RiS4TONhhzHW3zDDkuIYs/o5NX8HE2uEz8JIJnGCxeQVL28ySRQyafoZPf4WWbnRkgE1sieNsdZPqt2H4SsTKfIin7OZIYSSaP08lreHnPzgyQiS0RK3NOkMxzPWfaKzRJKVZm1i9IYiQjeAwvr+PlXTvzziCZq8HbAePm0T7HP3M1Ps4HX7FDZj6KlzV47VPErpGJLRErcz7tcw054iAzkAF2WcVg6D6SSAgiEf/WPRfT8fA4ggsv2/GyCYP9V7EAiTDv0+CyfiiFZDwsw8V4DJftntgav8zuEiF7rF/mDDw8Zmduw8tmv8zuEuk5s9HOrLUzr5UI2WP8Mmfi4VE78z28vH01M0AizHvYL3OYnVkUPDNghSZr9EeZbmbZJ8C58PIuXrYEZF6VCHMeDMhcjotxGBrszDqHmbNx84id+Q5etl7NDJAIcx/86OjEkJmDFBWJU8JIpIugMqk9F0Qi3d/6oDLxVYWUyNXMIDKpPRtEIqEyG/0ye5bI1cwgMqk5E0QigZlBZOKrCrlCdxFUJjVngkgkWGbAiu2rdpgZRCbVp3uUSMjMQSwTFYkTHEqki24yaXkT74lnMZ2tPa7QXXRbsU0jHeV/j6n6IKRErmb6yaTlv/Ee/58Yb8/iCp5ZT0f51zBVu0NK5Gqmn0yu/BfeE9+xM4NL5Gqmn0xMPR1lf4+p2Rtyhe6im0yurMd74tsYb3uPEumeaa/Ypo6Osq9GkOknkyvr8J74X3ZmcIkEzRzkMlGROOHkdscS6cLFdDwdDyPlh/FmnsI3OTPkCt3FRyv2pTbMlXN0jDuKWbwspESuZs7A0/EgUnYIb9YpfJOzIsxsx1w5S+e4o5jFy0NK5GrmTDwdn7K+z6yT+CZnh5TI1UxbJhdbMS3n6Bx3DLN4RcgVugs3s3B3fMLOPIUpzg0pke6Zy3FdarEzP4wgczbujo8j5YfwZp/EFOeFlMg1mbZMOnn16ubRIEJF4oSJN0LBdMcS6cKVNAdP5t/iTr4J1/S/cXwNEyEZT+63cLmvwzPmn5Dsoggy5+HJ+oqd+eVeZbrH/BOSHcn3OR9P5pdxJ92Ma/pfR5A5zMr0zMI99p+RrPDi6sKddAPuzL/GnXwLMu0vI8vM+Sc7858izFyAe8RfWZnTv+T4jM2rMpmKm4cG5fEleoi8ogxhEuKarYqiKKAiURSlD1CRKIoSNXpkK0R3Gndvx2qmZvbF2ARBRQLRncbd27GaqZl9MTZB0E0bRVGiRkWiKErU6KYNRHcad2/HaqZm9sXYBEEPSFOUIUxCHJDmtERcRIpEZIOIHBGRUhGZEE2uoiiJRX+ViL8M/MAYMwNYRPAiLUVRBigxLxEXkZmAxxizEcAY02SMuRJlrqIoCUTMS8SBqUC9iKwRkX0i8gMRcQd7sUQtEVcUJTT9USLuAW4D5gPngFeBZ4BfBi6YqCXiiqKEpj9KxMuAfcaYU/aYtcBigohEUZSBScxLxIEPgGwRybPv3wWURpmrKEoCEfMScWOMF/gasFlEDgICPBdlrqIoCURUR7YaY2qAu4M8vhv4c7/7G4E50WQpipK46Lk2iqJEjYpEUZSoUZEoihI1KhJFUaJGRaIoStSoSBRFiRoViaIoUaMiURQlalQkiqJEjYpEUZSoUZEoihI1KhJFUaJGRaIoStSoSBRFiRoViaIoUaMiURQlahK2aU9ELgPH4j2PAHKB6nhPwg+dT2gSbT6QeHOaZozJiPZFErn791hfVAn2JSKyO5HmpPMJTaLNBxJvTiLSJ724ummjKErUqEgURYmaRBbJL+I9gSAk2px0PqFJtPlA4s2pT+aTsDtbFUUZOCTyJxJFUQYIKhJFUaImriIRkZEislFEjttfs4MsM09E/iQih0XkgIh8xu+5lSJyWkT22//m9XIe94vIMRE5ISJfD/L8MBF51X5+l4hM8HvuG/bjx0Tkvt7k92I+XxWRUvv92Cwi4/2e8/q9H+v7Yj4O5/SMiFT5Zf+533NP2//Hx0Xk6cCxMZrPD/3m8qGI1Ps91+fvkYi8ICKVInKoh+dFRH5kz/eAiFzv91ws3p9w81lmz+OAiOwQkbl+z50RkYP2++Psz8PGmLj9A74PfN2+/XXgX4MsMxWYYt8eDVwAsuz7K4FHo5yDGzgJTAKSgRJgZsAyXwJ+bt/+LPCqfXumvfwwYKL9Ou5+mM+dQJp9+y+75mPfb4rB/5OTOT0D/DjI2JHAKftrtn07O9bzCVj+y8ALMX6PPgZcDxzq4fmPA3/AqqxdDOyK1fvjcD43d+UAD3TNx75/BsiNJC/emzYPAS/Zt18CHg5cwBjzoTHmuH27AqgE8gKXi4JFwAljzCljTDvwG3tePc3zNeBuERH78d8YY9qMMaeBE/brxXQ+xpgtxpgr9t2dwNgoM6OeUwjuAzYaY2qNMXXARuD+fp7PE8Cvo8wMiTHmXaA2xCIPAS8bi51AlogUEpv3J+x8jDE77Dzog5+heItklDHmAoD9NT/UwiKyCOs30Em/h79rfzz7oYgM68UcxgDn/e6X2Y8FXcYY0wk0ADkOx8ZiPv58Hus3XRcpIrJbRHaKyDVijvGcltr/F6+JyLgIx8ZiPtibfROBt/0ejsV7FI6e5hyL9ydSAn+GDLBBRPaIyBecvEDMD5EXkU1AQZCnvhnh6xQCq4CnjTE+++FvABex5PIL4B+Bb0c6xSCPBf5NvKdlnIyNFMevKSLLgQXA7X4PFxljKkRkEvC2iBw0xpwMNr6P5/RfwK+NMW0i8kWsT3B3ORwbi/l08VngNWOM1++xWLxH4ejPnyHHiMidWCK51e/hW+z3Jx/YKCJH7U84PRLzTyTGmCXGmNlB/q0DLtmC6BJFZbDXEJERwH8D37I/Fna99gX7o2Ib8CK926woA8b53R8LVPS0jIh4gEysj41OxsZiPojIEiwZP2h//8BHm38YY04BW4H5Uc7H0ZyMMTV+83gOuMHp2FjMx4/PErBZE6P3KBw9zTkW748jRGQO8DzwkDGmputxv/enEngDJ+tVX+90inCH0A/ovrP1+0GWSQY2A38b5LlC+6sA/w58rxdz8GDt4JrI1R13swKW+Su672z9rX17Ft13tp4i+p2tTuYzH2vzbkrA49nAMPt2LnCcEDsh+3hOhX63Pw3stG+PBE7bc8u2b4+M9Xzs5aZh7TiUWL9H9utNoOedm5+g+87W92P1/jicTxHWPr2bAx4fDmT43d4B3B82qy8mHMU3mmNL4rj9daT9+ALgefv2cqAD2O/3b5793NvAQeAQsBpI7+U8Pg58aK+c37Qf+zbWb3uAFOB39hv/PjDJb+w37XHHgAf66H0JN59NwCW/92O9/fjN9vtRYn/9fB/+X4Wb078Ah+3sLcB0v7F/Zr93J4DP9cd87PvPEvDLJVbvEdanngv2z2oZ1ubCF4Ev2s8L8BN7vgeBBTF+f8LN53mgzu9naLf9+CT7vSmx/z+/6SRPD5FXFCVq4v1XG0VRBgEqEkVRokZFoihK1KhIFEWJGhWJoihRoyJRFCVqVCSKokTN/wc3tqIOcGFGQwAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "dico_dom = {\n", " 'space_step': .1,\n", " 'box': {'x': [0, L], \n", " 'y': [-.5*W, .5*W], \n", " 'label': list(range(4))\n", " },\n", " 'schemes': [\n", " {'velocities': list(range(9))}\n", " ],\n", "}\n", "dom = pylbm.Domain(dico_dom)\n", "print(dom)\n", "dom.visualize(view_distance=True);" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### The scheme\n", "\n", "In pylbm, a simulation can be performed by using several coupled schemes. In this example, a single scheme is used and defined through a list of one single dictionary. This dictionary should contain:\n", "\n", "* 'velocities': a list of the velocities\n", "* 'conserved_moments': a list of the conserved moments as sympy variables\n", "* 'polynomials': a list of the polynomials that define the moments\n", "* 'equilibrium': a list of the equilibrium value of all the moments\n", "* 'relaxation_parameters': a list of the relaxation parameters ($0$ for the conserved moments)\n", "* 'init': a dictionary to initialize the conserved moments\n", "\n", "(see the documentation for more details)\n", "\n", "In order to fix the bulk ($\\mu$) and the shear ($\\eta$) viscosities, we impose\n", "$$ s_\\eta = \\frac{2}{1+\\eta d}, \\qquad s_\\mu = \\frac{2}{1+\\mu d}, \\qquad d = \\frac{6}{\\lambda\\rho_0\\dx}.$$\n", "\n", "The scheme velocity could be taken to $1$ and the inital value of $\\rho$ to $\\rho_0=1$, $q_x$ and $q_y$ to $0$.\n", "\n", "In order to accelerate the simulation, we can use another generator. The default generator is Numpy (pure python). We can use for instance Cython that generates a more efficient code. This generator can be activated by using 'generator': pylbm.generator.CythonGenerator in the dictionary." ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "+--------------------+\n", "| Scheme information |\n", "+--------------------+\n", " - spatial dimension: 2\n", " - number of schemes: 1\n", " - number of velocities: 9\n", " - conserved moments: [ρ, qx, qy]\n", " \n", " +----------+\n", " | Scheme 0 |\n", " +----------+\n", " - velocities\n", " (0: 0, 0)\n", " (1: 1, 0)\n", " (2: 0, 1)\n", " (3: -1, 0)\n", " (4: 0, -1)\n", " (5: 1, 1)\n", " (6: -1, 1)\n", " (7: -1, -1)\n", " (8: 1, -1)\n", "\n", " - polynomials\n", " \n", " ⎡ 1 ⎤\n", " ⎢ ⎥\n", " ⎢ LA⋅X ⎥\n", " ⎢ ⎥\n", " ⎢ LA⋅Y ⎥\n", " ⎢ ⎥\n", " ⎢ 2 2 ⎥\n", " ⎢ 3⋅X + 3⋅Y - 4 ⎥\n", " ⎢ ⎥\n", " ⎢ 2 ⎥\n", " ⎢ 2 2 ⎛ 2 2⎞ ⎥\n", " ⎢ 21⋅X 21⋅Y 9⋅⎝X + Y ⎠ ⎥\n", " ⎢- ───── - ───── + ──────────── + 4⎥\n", " ⎢ 2 2 2 ⎥\n", " ⎢ ⎥\n", " ⎢ ⎛ 2 2⎞ ⎥\n", " ⎢ 3⋅X⋅⎝X + Y ⎠ - 5⋅X ⎥\n", " ⎢ ⎥\n", " ⎢ ⎛ 2 2⎞ ⎥\n", " ⎢ 3⋅Y⋅⎝X + Y ⎠ - 5⋅Y ⎥\n", " ⎢ ⎥\n", " ⎢ 2 2 ⎥\n", " ⎢ X - Y ⎥\n", " ⎢ ⎥\n", " ⎣ X⋅Y ⎦\n", "\n", " - equilibria\n", " \n", " ⎡ ρ ⎤\n", " ⎢ ⎥\n", " ⎢ qx ⎥\n", " ⎢ ⎥\n", " ⎢ qy ⎥\n", " ⎢ ⎥\n", " ⎢ 2 2⎥\n", " ⎢ 3.0⋅qx 3.0⋅qy ⎥\n", " ⎢-2⋅ρ + ─────── + ───────⎥\n", " ⎢ 2 2 ⎥\n", " ⎢ LA LA ⎥\n", " ⎢ ⎥\n", " ⎢ 2 2 ⎥\n", " ⎢ 3.0⋅qx 3.0⋅qy ⎥\n", " ⎢ ρ - ─────── - ─────── ⎥\n", " ⎢ 2 2 ⎥\n", " ⎢ LA LA ⎥\n", " ⎢ ⎥\n", " ⎢ -qx ⎥\n", " ⎢ ──── ⎥\n", " ⎢ LA ⎥\n", " ⎢ ⎥\n", " ⎢ -qy ⎥\n", " ⎢ ──── ⎥\n", " ⎢ LA ⎥\n", " ⎢ ⎥\n", " ⎢ 2 2 ⎥\n", " ⎢ 1.0⋅qx 1.0⋅qy ⎥\n", " ⎢ ─────── - ─────── ⎥\n", " ⎢ 2 2 ⎥\n", " ⎢ LA LA ⎥\n", " ⎢ ⎥\n", " ⎢ 1.0⋅qx⋅qy ⎥\n", " ⎢ ───────── ⎥\n", " ⎢ 2 ⎥\n", " ⎣ LA ⎦\n", "\n", " - relaxation parameters\n", " \n", " ⎡ 0.0 ⎤\n", " ⎢ ⎥\n", " ⎢ 0.0 ⎥\n", " ⎢ ⎥\n", " ⎢ 0.0 ⎥\n", " ⎢ ⎥\n", " ⎢1.13122171945701 ⎥\n", " ⎢ ⎥\n", " ⎢1.13122171945701 ⎥\n", " ⎢ ⎥\n", " ⎢0.025706940874036⎥\n", " ⎢ ⎥\n", " ⎢0.025706940874036⎥\n", " ⎢ ⎥\n", " ⎢0.025706940874036⎥\n", " ⎢ ⎥\n", " ⎣0.025706940874036⎦\n", "\n", " - moments matrices\n", " \n", " ⎡1 1 1 1 1 1 1 1 1 ⎤\n", " ⎢ ⎥\n", " ⎢0 LA 0 -LA 0 LA -LA -LA LA ⎥\n", " ⎢ ⎥\n", " ⎢0 0 LA 0 -LA LA LA -LA -LA⎥\n", " ⎢ ⎥\n", " ⎢-4 -1 -1 -1 -1 2 2 2 2 ⎥\n", " ⎢ ⎥\n", " ⎢4 -2 -2 -2 -2 1 1 1 1 ⎥\n", " ⎢ ⎥\n", " ⎢0 -2 0 2 0 1 -1 -1 1 ⎥\n", " ⎢ ⎥\n", " ⎢0 0 -2 0 2 1 1 -1 -1 ⎥\n", " ⎢ ⎥\n", " ⎢0 1 -1 1 -1 0 0 0 0 ⎥\n", " ⎢ ⎥\n", " ⎣0 0 0 0 0 1 -1 1 -1 ⎦\n", "\n", " - inverse of moments matrices\n", " \n", " ⎡1/9 0 0 -1/9 1/9 0 0 0 0 ⎤\n", " ⎢ ⎥\n", " ⎢ 1 ⎥\n", " ⎢1/9 ──── 0 -1/36 -1/18 -1/6 0 1/4 0 ⎥\n", " ⎢ 6⋅LA ⎥\n", " ⎢ ⎥\n", " ⎢ 1 ⎥\n", " ⎢1/9 0 ──── -1/36 -1/18 0 -1/6 -1/4 0 ⎥\n", " ⎢ 6⋅LA ⎥\n", " ⎢ ⎥\n", " ⎢ -1 ⎥\n", " ⎢1/9 ──── 0 -1/36 -1/18 1/6 0 1/4 0 ⎥\n", " ⎢ 6⋅LA ⎥\n", " ⎢ ⎥\n", " ⎢ -1 ⎥\n", " ⎢1/9 0 ──── -1/36 -1/18 0 1/6 -1/4 0 ⎥\n", " ⎢ 6⋅LA ⎥\n", " ⎢ ⎥\n", " ⎢ 1 1 ⎥\n", " ⎢1/9 ──── ──── 1/18 1/36 1/12 1/12 0 1/4 ⎥\n", " ⎢ 6⋅LA 6⋅LA ⎥\n", " ⎢ ⎥\n", " ⎢ -1 1 ⎥\n", " ⎢1/9 ──── ──── 1/18 1/36 -1/12 1/12 0 -1/4⎥\n", " ⎢ 6⋅LA 6⋅LA ⎥\n", " ⎢ ⎥\n", " ⎢ -1 -1 ⎥\n", " ⎢1/9 ──── ──── 1/18 1/36 -1/12 -1/12 0 1/4 ⎥\n", " ⎢ 6⋅LA 6⋅LA ⎥\n", " ⎢ ⎥\n", " ⎢ 1 -1 ⎥\n", " ⎢1/9 ──── ──── 1/18 1/36 1/12 -1/12 0 -1/4⎥\n", " ⎣ 6⋅LA 6⋅LA ⎦\n", "\n", " \n" ] } ], "source": [ "import sympy as sp\n", "X, Y, rho, qx, qy, LA = sp.symbols('X, Y, rho, qx, qy, LA')\n", "\n", "# parameters\n", "dx = 1./128 # spatial step\n", "la = 1. # velocity of the scheme\n", "L = 1 # length of the domain\n", "W = 1 # width of the domain\n", "rhoo = 1. # mean value of the density\n", "mu = 1.e-3 # shear viscosity\n", "eta = 1.e-1 # bulk viscosity\n", "# initialization\n", "xmin, xmax, ymin, ymax = 0.0, L, -0.5*W, 0.5*W\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", "dico_sch = {\n", " 'box': {'x': [xmin, xmax], \n", " 'y': [ymin, ymax], \n", " 'label':0\n", " },\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", "}\n", "sch = pylbm.Scheme(dico_sch)\n", "print(sch)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Run the simulation\n", "\n", "For the simulation, we take\n", "\n", "* The domain $x\\in(0, L)$ and $y\\in(-W/2,W/2)$, $L=2$, $W=1$,\n", "* the viscosities $\\mu=10^{-2}=\\eta=10^{-2}$, \n", "* the space step $\\dx=1/128$, the scheme velocity $\\lambda=1$, \n", "* the mean density $\\rho_0=1$.\n", "\n", "Concerning the boundary conditions, we impose the velocity on all the edges by a bounce-back condition with a source term that reads\n", "$$q_x(x, y) = \\rho_0 v_{\\text{max}} \\Bigl( 1 - \\frac{4y^2}{W^2} \\Bigr), \\qquad q_y(x, y) = 0,$$\n", "with $v_{\\text{max}}=0.1$.\n", "\n", "We compute the solution for $t\\in(0,50)$ and we plot several slices of the solution during the simulation. \n", "\n", "This problem has an exact solution given by\n", "$$ q_x = \\rho_0 v_{\\text{max}} \\Bigl( 1 - \\frac{4y^2}{W^2} \\Bigr), \\qquad q_y = 0, \\qquad p = p_0 + K x, $$\n", "where the pressure gradient $K$ reads\n", "$$K = -\\frac{8 v_{\\text{max}} \\eta}{W^2}.$$\n", "\n", "We compute the exact and the numerical gradients of the pressure." ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAdAAAAEICAYAAADvMKVCAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzsnXecFEX2wL+PXYKAIBgRUUD0FM+MomfACOiJCVQ8f4oKYj6V8xSzhxjPcOepp4gBMeGBAfV0RRE5EwIiUYEFlCCKBJGgwLLv98frgWGY3Qk7Mz07+76fT32mu7qq+vX0TL+uqlfviariOI7jOE5q1ApbAMdxHMepjrgCdRzHcZw0cAXqOI7jOGngCtRxHMdx0sAVqOM4juOkgStQx3Ecx0kDV6BOVhGRG0VkYJbP8Y6I9Ai2zxeRj6OOqYi0SbPd/iKyWER+EJGWQVvFmZLbcZzqjStQJylE5FsR+VVEVorIjyLyjIg0TFRPVe9S1V7ZlE1VT1DVQZlsU0RaAH8B2qrqDplsO015jhKR+QnKPCsi/bN0fhWRVcH9Xxn9UiTGvSKyJEj3iYhU0tafROS7oL3XRaRp1LGmIvJacOw7EflTPtR1nHi4AnVSoYuqNgQOAA4Cbg5ZnmyyC7BEVReFLUgesa+qNgxS9EtRb+BUYF9gH+Ak4OJ4DYjIXsATwLnA9sBq4LGoIo8Ca4Nj5wD/DuqEVtdxKkRVPXlKmIBvgeOi9v8OvBVs7wgMB5YCpcBFUeVuB54PtusBzwNLgJ+BscD2wbHGwFPAQmAB0B8oim0j2G8JKFAc7I8CegXb5wMfR5VVoE2wXRe4H5gL/Ag8DmwR51qPA34FyoGVwLNxzhn3moNr/BXYJti/GSgDGgX7/YF/VPAdXwB8DawAZgMXB/kNYuRZCewYU7c3sA5TAiuBNzN8/zd8j3GOfQr0jtrvCXxeQdm7gBej9ncNZN4yuM61wO5RxwcD94RZ15OnipL3QJ2UCYY3TwQmBFkvAfMxpdINuEtEjo1TtQemKFsAWwOXYIoBYBCmaNoA+wMdgUwP/d4L7A7sF5ynOXBrbCFVfR84Afherbd1fpy24l6zqv6GvRh0CModCXwHHBa1/1EF8i3Cem+NMGX6kIgcoKqrYuRpqKrfx8g8AHgBuC843iXeCURkkoj8XEF6LF6dKEYH88GvikjLqPy9gIlR+xODvHhsUlZVZxEoryCtV9UZFbQVVl3HiYsrUCcVXheRn4GPMSVwV6BMDweuV9XfVPUrYCA2VBbLOkxxtlHV9ao6XlV/EZHtMQVxtaquUhs2fQjoninBgzm5i4BrVHWpqq7AeiUpnyOJa/4I6BAYHO0DPBzs18OGvv8Xr11VfVtVZ6nxEfAecESq8lWGqu6jqltVkC6rpGoHrBe+B/A98FaUQVVDYHlU2eVAwwrmQWPLRspvmeBYmHUdJy5uUeikwqlB72wDIrIjEFFIEb4D2sWpPxjrfb4sIlthw7k3YfONtYGFUc/cWsC8DMq+LVAfGB91DgGK0mgr0TV/BDyIzRVPBkZgw9OHAKWqujheoyJyAnAb1iOqFcg7OQ35Mo6qjg4214rIVcAvwJ6YfCuxXnOERsBKVY0XqSK2bKT8CmyIuqJjYdZ1nLh4D9SpKt8DTUUk+m19Z2wecxNUdZ2q/k1V2wJ/wIYrz8MU5Rps3jDSG2qkqpEhtFWYMomQjlXsYmy4eK+oczRWM4pKlUTX/CnwO+A04CNVnRYc/yMVDN+KSF1gGDZHu72qbgX8F1PyYHOQiUhYRkSmRlnSxqbHkzhH9Lkisk3FDIgi7BvkxWOTsiLSGpubnhGkYhHZrYK2wqrrOHFxBepUCVWdhymMu0WknojsgxmRvBBbVkSOFpG9RaQI68Gsw+aeFmLDlQ+ISCMRqSUiu4pIZB7xK+BIEdlZRBoDN6QhZznwJDavuF0gT3MR6ZTpa1bV1cB44HI2KsxPMcvUiuY/62AP9J+AsqA32jHq+I/A1sH1V8SPQOsEsu8VNY8amy6JV0dE9hKR/USkSGzp0gPYy8LXQZHngD7B97kjtvzn2QpEeAHoIiJHiEgDoB/wqqquCOZ6XwX6iUgDETkMOAUbuQizruPExRWokwnOxubHvgdeA25T1RFxyu0ADMWU59eYMnk+OHYepkSmAcuCcs0AgraGAJMwxfRWmnJej1nMfi4ivwDvYz3FdEh0zR9hw9JfRO1vCYwmDsFw8J+BV7Dr/xNm5Rs5/g1muDQ7MPjZMU4zTwFtg+Ovp3ld8dge+/5/wayDWwInqeq64PgTwJvYcO4U4O0gD4Cgd3tEcB1TMeOxFzCjqS2B6LnXy4AtgmMvAZcGdUKr6zgVIfGnKRzHcRzHqQzvgTqO4zhOGoSqQEXoLMJ0EUpF6Bvn+JEifClCmQjdYo6tF+GrIA2Pym8lwhgRZoowRIQ6ubgWx3Ecp2YRmgIVoQhzn3UC0BY4W4S2McXmYp5lXozTxK+q7Bekk6Py7wUeUmU3bC6pZ8aFdxzHcWo8YfZADwZKVZmtylrgZczybQOqfKvKJGydVkJEEOAYzAAFzLvNqZkT2XEcx3GMMB0pNGfThfLzgfYp1K8nwjjM/ds9qryOebn5WZWyqDabx6ssIr0x/6EAB9avXz9eMcdxHKcCVq9erapaY21pwlSg8dx8pWISvLMq34vQGhgpwmTMzD6pNgPfoQMAGjRooKtWrUrh1I7jOI6I/Jq4VOES5pvDfMytW4SdsDV1SaFqZVWZjUXj2B/zNrOVyIYXg5TadBzHcZxkCVOBjgV2C6xm62BOvYcnqAOACE1EqBtsb4NFupimigIfwgaL3R7AGxmX3HEcx6nxhKZAg3nKK4ASzCvNK6pMFaGfiFnVinCQCPOBM4AnRDb4ptwTGCfCRExh3qPKtODY9UAfEUqxOdGncndVjuM4ziaIdEZkOiKliGy2XBGRuogMCY6PIRIqT2RrRD5EZCUij8TUGRW0+VWQtsvBlWyGeyLC50Adx3HSQURWq2qDSgoUYc76j8em7cYCZ2MBFiJlLgP2QfUSRLoDp6F6FuazeH/g98DvUb0iqs4o4FpUx2X8olKgxlpPOY7jOFnnYKAU1dmoxl2uGOwPCraHAsciIqiuQvVj4LfciZsarkAdx3GcdCkWkXFRqXfM8XjLFWOXFm4so1qGBTPfOolzPxMM395C/ODtWccDajuO4zjpUqaq7So5nsxyxXSWNJ6D6gIsJu8w4FwsrF5O8R6o4ziOky2SWa64sYxIMdAYWFppq6oLgs8VmKvXgzMhbKq4AnUcx3GyxVhgN0RaIVLRcsXh2JJDsCWII6nMulWkGJFtgu3awElYHNqc40O4juM4TnZQLUMkslyxCHga1amI9APGoTocW2o4GJFSrOfZfUN9kW+BRkAdRE4FOgLfASWB8iwC3geezN1FbcSXseDLWBzHcdIh4TKWAseHcB3HcRwnDVyBOo7jOE4auAJ1HMdxnDRwBeo4juM4aeAK1HEcx3HSwBWo4ziO46SBK1DHcRzHSQNXoI7jOI6TBq5AHcdxHCcNXIE6juM4Thq4AnUcx3GcNAhVgYrQWYTpIpSK0DfO8SNF+FKEMhG6ReXvJ8JnIkwVYZIIZ0Ude1aEOSJ8FaT9cnU9juM4Ts0htGgsIhQBjwLHY/HgxoowXJVpUcXmAucD18ZUXw2cp8pMEXYExotQosrPwfG/qjI0u1fgOI7j1GTCDGd2MFCqymwAEV4GToGNClSVb4Nj5dEVVZkRtf29CIuAbWGDAnUcx3GcrBLmEG5zYF7U/vwgLyVEOBioA8yKyr4zGNp9SIS6VRPTcRzHcTYnTAUqcfJSCk4qQjNgMHCB6oZe6g3AHsBBQFPg+vh1pbeIjBORcWVlZamc1nEcx3FCVaDzgRZR+zsB3ydbWYRGwNvAzap8HslXZaEqqsoa4BlsqHgzVHWAqrZT1XbFxWGOZDuO4zjVkTAV6FhgNxFaiVAH6A4MT6ZiUP414DlV/hNzrFnwKcCpwJSMSu04juM4hKhAVSkDrgBKgK+BV1SZKkI/EU4GEOEgEeYDZwBPiDA1qH4mcCRwfpzlKi+IMBmYDGwD9M/hZTmO4zg1BFFNadqxIGnQoIGuWrUqbDEcx3GqFSKyWlUbhC1HWLgnIsdxHMdJA1egjuM4jpMGrkAdx3EcJw1cgTqO4zhOGrgCdRzHcZw0cAXqOI7jZA+RzohMR6QUkc2ibiFSF5EhwfExiLQM8rdG5ENEViLySEydAxGZHNR5GJF4nu2yjitQx3EcJzuIRKJunQC0Bc5GpG1MqZ7AMlTbAA8B9wb5vwG3sHk0LoB/A72B3YLUOfPCJ8YVqOM4jpMtDgZKUZ2N6lrYEHUrmlOAQcH2UOBYRATVVah+jCnSjYg0Axqh+hnmyOA5zOtcznEF6jiO46RLcSQoR5B6xxxPJurWxjKqZcByYOtKztk8aKeyNnOCe1F3HMdx0qVMVdtVcjyZqFupRuaqciSvTOE9UMdxHCdbJBN1a2MZkWKgMbA0QZs7JWgzJ7gCdRzHcbLFWGA3RFohUlHUreFAj2C7GzCSypy0qy4EViBySGB9ex7wRsYlTwIfwnUcx3Gyg2oZIpGoW0XA06hORaQfMA7V4cBTwGBESrGeZ/cN9UW+BRoBdRA5FeiI6jTgUuBZYAvgnSDlHI/GgkdjcRzHSQePxuI4juM4Tsq4AnUcx3GcNHAF6jiO4zhp4ArUcRzHcdLAFajjOI7jpEGoClSEziJMF6FUhM289ItwpAhfilAmQreYYz1EmBmkHlH5B4owOWjzYZG4Xiscx3Ecp0qEpkBF2MxLvwixXvrnAucDL8bUbQrcBrTHnBXfJkKT4HBeeOl3HMdxCpswHSkcDJSqMhtAZIOX/mmRAqp8Gxwrj6nbCRihau6eRBgBdBZhFNBIlc+C/IiX/lAW2eYrM2fC44/DihWbH5szBxYvhtatQQTWrNmYfvsNyspg/XpL0dvr10NkSbHqxhRNJGJf9Gd0qlWr8v3ovHj1o1Ps8YrOH/0ZK2e8Y7lk0SL44Qcoj/31B2yxBeyyi31CfLkry0vms7LvuqL8eKmiexvvs7LtylJR0aafsXmxx6PLRe/Hk18VvvoKli2z+xGb6tWDfv2gYcP07rVTPQlTgcbz0t++CnWbk4KX/iBqQG+AOnXqJHna6s2338Idd8CgQfaw2DqId1Cs6yiT2gCsXGmKdeJEO77jjlC/PtStC1ttBbVrb3wYFRdv3C4qqlx5RSvXyGdsKi9PnBdRJvHqRyvtePuVfUaI3g/Lx8jixfYis3q1PZjr1o1fLqJgd9gBdt55Y7l415bs9xD7/UVvV/YdJ5Pi3cvYz4ry4u1X9GIRFm+8Yf+vs86y/4NT+ISpQKviUb+iukm3qaoDgAFgnoiSPG+15Pvv4c474cknTaFdcQX07WsPXgA6d4FttrFCu+zC9Onwt7/Byy/D2rXQpw9ccw00ahTqZRQ0qvDf/8Ktt8LUqbD77nD77XDmmRU/jBcuhLvvhieegCVL4KKL4Kab7KWnpqBqox8RhRpvOzJCEpsXrYhj98vL4X//s+92xgxo0QJ69oS99rIXx9ie7IwZNqpzzjlw1112704/3Y47BYyqhpJADwUtidq/AfSGCso+C9otav9s0Cei9p8I8pqBflNRuYpS/fr1tRBZtEi1Tx/VevVUi4tVL75Yde7cmELr1qneeKMVqltX9frrVX/+WVVVJ01SPe006y80bap6zz2qq1bl/joKndGjVQ85xL7nVq1Un33WbkuyzJ2r2ru33eN69VSvuUZ12bLsyVvIlJerlpSotm+f+v1Yv1715ZdV99jD6u67r+obb1ibhQqwSkPSIfmQwlSgxaCzQVuB1gGdCLpXBWVjFWhT0DmgTYI0B7RpcGws6CGgAvoO6ImJZClEBfrtt6o776xaq5bq+eerzpqVoMLcuarnnWc/iW22Uf3Xv1TXrlVV1XHjVE880Q61b6+6fHn25a8pvPaaalGR6k47qT7xxIavPC1mzbJ7XauW6oEHuhJNlXXrVP/0J/udt2ihOmBAevejrEz1uedUd93V2urSRfW33zIvbz7gCjTMk6Mngs4AnQV6U5DXD/TkYPsg0Pmgq0CXgE6NqnshaGmQLojKbwc6JWjzEVBJJEehKdB581Rbt1bdaivVL75IsfL48apHH20/jd13tyd88Ao9dKj1cg4/XHXlyszLXdN46y3V2rUz/1ISafeQQ1R/+SVz7RYyZWWq55xjP/vbbsuMwlu7VvX++63Nk09WXbOm6m3mG65A80CIsFMhKdCFC03vbbml6pgxaTZSXq765psbx6IOP1z1k09UVXXIEOvhHHOM6urVmZO7plFSYiPm2eopvvqq9Wz9ZScx69dbzx1U77or8+0/+qi13bVr1UYY8hFXoHkgRNipUBTookWqbduqNmig+vHHGWhw3TrVxx9X3WEH+6mceqrqtGk6eLCqiGrnzoU7NJVNRo60ucp991VdsiR754m87Bx9tM9dV8T69aq9etnP+/bbs3eehx6yc3Tvntr8dr7jCjQPhAg7FYICXbLEHsj16tkDOqOsXKnav791a2vVUu3VS1/8+/wNQ1OF9ladTUaPVq1fX3WvveyFJ9s8/7y97HTsqPrrr9k/X3WivFz1ssvsKXjTTdk39rnvPjvXuefakHEh4Ao0D4QIO1V3Bfrzz6rt2qnWqWNDg1lj0SLVq66yCbYtttDxHftqY5bpGWcU1lt1tvjsM9WGDW2IfeHC3J33qafsn/7HPxbmPFw6lJfbTxlU//rX3FnK9u9v57zwQuv9VndcgeaBEGGn6qxAV6xQPfRQM+55880cnXTWrA3mir9usZVexz16wVmrCuKBkC0mTFBt3NgsM+fPz/35H3/c/u2nnVYYD+6qUF5uShNMieZ6mcmtt9q5L764+i9xcQWaB0KEnaqzAr34YhtVHTo0hJNPmLBhfcsCmunIMx7z8dw4/PabzU03b6763XfhyfHgg/aPf+SR8GTIB5591r6Hyy4LR4GVl6v27WsyPPRQ7s+fSVyB5oEQYafqqkBHjbI7eO214cpR/tFondzkcFXQtTu3Vn3hBe/mRBHpcbz9drhylJfbXGjDhrZOuCby/fe2vOuII8L9iZaX25B6/fpJrNHOY2q6AhX7Dmo2DRo00FWrVoUtRkr8+ivss4+5Mps0yfzVhsnc75Rr9niHvxffQOuVk2DvvaF/f+jSJVyP7CEzaRIceCB07w6DB4ctjflD/v3v4YgjzHVgTbs1XbvadU+caO4Sw2TePHMNeNBB8P771fNeiMhqVW0Qthxh4Z4aqym33w6lpebfNmzlCbDzLsIx959Im5UTGH3Jixa65ZRT4JBD7OlQA1/UysrMf2qTJvDQQ2FLY7Rsaf5z330Xnn8+bGlyy9Ch8Oqr5uc5bOUJ5l/373+HkSPhqafClsZJi7C7wPmQqtsQ7rhxtki+V6+wJdmU9etVDzvM/Ob+MH+d6sCB5hMNVI86aoMzhprC3/9ul/7yy2FLsimb3KcfwpYmNyxerLrddua4Ip8sxtevt79Go0bhGJdVFWr4EG7oAuRDqk4KdO1aW+/ZrFl++jr9+mtbTnPmmUHGb7+p/vOf9vQCMzr68stQZcwFM2famtxTTslPS8vIfTrjjLAlyQ3nnmuW6l99FbYkmzNzpuoWW5jP3Hz8rVSGK9A8ECLsVJ0U6J132l177bWwJamYyFq311+Pyly5UvXuu1WbNNENfs2mTAlNxmyyfr1qhw62bGXBgrClqZjIb+nVV8OWJLu8/bZd5y23hC1JxUR85r74YtiSpIYr0DwQIuxUXRTo11+b/9R87zWsXau6zz6qO+4Yp5f888/mrXvLLc1FzjnnqM6YEYaYWeOJJ+yf9eSTYUtSOWvXqu63n3lqXLo0bGmyw/LlNovQtm1+u50sK1M9+GALhJQLD1WZIikFCp0VpiuUKvSNc7yuwpDg+BiFllHHbgjypyt0isr/VmGywlcK4xLK4Aq0ZivQ9evNMXiTJtVj3mrsWFufetFFFRRYvNhij9avbxO6F16oOmdOLkXMCvPm2XzW0UdXj+G48ePt6+/ZM2xJssOll9rv8PPPw5YkMVOmmJOvs88OW5LkSahAoUhhlkJrhToKExXaxpS5TOHxYLu7wpBgu21Qvq5Cq6CdIt2oQLep9NyuQF2BRohEdHj22bAlSZ6It5dKffP+8IPq1Vdb17p2bdVLLjEtVA0pL1c96SSbzyotDVua5Iks6h8xImxJMktknXSfPmFLkjx/+5vJPHx42JIkRxIK9FCFkqj9GxRuiClTonBosF2ssFhBNiu7abm8UKC+DpT8Xwe6ZAm0amUrQkpKqs96sdWrba1q7dowZQoUFVVSeP58uOsuGDgQatWCiy+Gvn2hWbOcyVtVRoyAjh1tacK114YtTfL8+ivsuy/UqWPrVmsVwOK2sjJb77puHUyenB9LvZJh7Vpo1w6WLYOZM6FevbAlqhwRWQtMjsoaoKoDogp0Azqj2ivYPxdoj+oVUWWmBGXmB/uzgPbA7cDnqD4f5D8FvIPqUETmAMsABZ4g+pw5pAD+KoXPww/DihXw4IPVR3mCPbTuvhu++QaGDUtQeKed4LHH7Klx7rnw6KOw666miX76KSfyVpU77oDmzeHKK8OWJDW22MLWFU+dCq+/HrY0mWHIEJg+3V5mqovyBHuJeeghe598+umwpUmKMlVtF5ViFVm8J1Zsr62iMpXVPQzVA4ATgMsROTIlqTOEK9A8Z/lyU6Cnn25v1NWNrl1hzz3NKVF5eRIVdtnFvENMnw5nnGFPk1at4IYbrCuep4weDf/7H1x3HdStG7Y0qXPWWbDbbnafqvugVHk53Hmn/V9OPTVsaVLnmGPg0EPhnnusR1rNmQ+0iNrfCfi+wjIixUBjYGmldVUjn4uA14CDMy14UoQ9hpwPKZ/nQO+6y+ZExo8PW5L0GTzYruGNN9Ko/M03ZlUhYk5cb745L01GjzvOlrpW58DVTz9t9+mtt8KWpGq88opdx0svhS1J+rzzjlYLS24Sz4EWK8wOjIAiRkR7xZS5PMaI6JVge68YI6LZgVFSA4UtgzINFD5V6FypHFlKoSou0M6g00FLQTczbwatCzokOD4GtGWQfw7oV1GpHHS/4NiooM3Ise0SyZGvCnTlSjNrP/HEsCWpGuvWqbZurXrQQVWwTJ0yxdbvgJm53n67LYnJAz77zMS6776wJakaa9eq7rKLavv21cOCOB7r19sSqt/9rnoHrS4vtxi/rVvnl+ekWBIqUFNyJyrMCKxobwry+imcHGzXU/hPsFzlC4XWUXVvCupNVzghyGsdKNaJClM3tFmTFChoEegs0NagdUAngraNKXMZ6OPBdnfQIXHa2Rt0dtT+KNB2qciSrwr0gQfsDn36adiSVJ0nn7RrqXLA74kTLaglWFiNO+6wxX4h8sc/mlu8FStCFSMj/PvfWq0tcl9/3eQfNChsSarOG2/k/7UkpUALOIVmhSvCocDtqnQK9m8AUOXuqDIlQZnPRCgGfgC2Vd04CS3CXdh7wE3B/ijgWlXGJStLPlrh/vabTf21bQsffBC2NFVn7VqzCWrVyuYLq8yECXDbbfDmm9C0KfzlL2a9s+WWGWg8NTEOOMAMiG6+Oaenzgpr1kDr1jYfOmpU2NKkhqpFNlm2zKbQi4vDlqhqqML++5uV9LRpCazYQ8KjsYRHc2Be1P78IC9uGVXKgOXA1jFlzgJeisl7RoSvRLhFJK4lFyLSW0TGici4srKydK8hazz9NPzwQ2E8lMGsC6+/3gxtPvooAw3uvz8MHw5ffGHre266ybTzvffCypUZOEFy9O8PjRtXP8vbiqhb1wyhPvrI7lV14t13Yfx4szer7soTzOL+5pthxgz4z3/ClsaJS1hdX9AzQAdG7Z8L+q+YMlNBd4ranwW6ddR+e9DJMXWaB59bgr4Hel4iWfJtCHfNGtWdd1b9wx+q71xUPFavVt1+ezO4yTiff67aubONeW2zjU1IrlyZhRNtZMoUO93NN2f1NDln1SoziOrYMWxJkqe8XPXQQ+1/s2ZN2NJkjvXrzQ3hXnvlZ4x6avgQbpg90JTMm4Mh3Ih5c4TuxPQ+VVkQfK4AXiQs8+YqMHgwzJ1rb5/Vad1nIrbYwpZ1vv8+fP55hhtv3x7eeQc+/dTGVK+7zsYiH3jAPDpkgTvvhIYN4eqrs9J8aNSvbyPi771nHfzqwMiR8Nln5nujTp2wpckctWrZ4MrUqfDaa2FL42xGWJobtBh0NmirKCOivWLKXB5jRPRK1LFaoPNBW8e0uU2wXRt0KOgliWTJpx7ounWqu+6qesABhdX7jLBihRncnHRSlk/08cfW1QXr9j74oHWBM8T06eZj9brrMtZkXvHLL3afunQJW5Lk6NDBghf8+mvYkmSesjLV3XYzx//59kyghvdAwz05eiLojGBo9qYgrx/oycF2PdD/BMtYvohRlkeBfh7TXgPQ8aCTguHff4IWJZIjnxTo88/bXSnkEFN33GHXmJOwoKNHqx5zjJ1whx1UH3ooI4r0/PPN5211cOyfLv362dc2YULYklTO6NEm5z/+EbYk2eOZZzQvfeTWdAXqvnDJHyvc8nLznlJUBBMnFoZP0ngsX24Oh44/PofGEaNHm7+6Dz80/7p9+8JFF9m4cop8+y20aQNXXAH/+EfGJc0bfv7Z7lPHjvltxNKxo/1f5sypXm77UmHdOth9d9h2WxgzJn+mdtwK18kbXnsNvv7a5jwKVXnCRqvVYcPMPD8nHHmkTZR9+KE9ia66ytbV/OtftmYoBe69115y/vrXLMmaJ2y11cb79PXXYUsTnzFjzIn/tdcWrvIEC8hw440wdqzNTTv5gfdAyZ8e6JFHwoIFZraej2u+MsnixdCiBVxwgfmQzzmjRtk60tGjYccdbe1Dr14Jw1+sWGEd2DPPrDbOvqvETz+Zn/9LL83P3vaf/mS2Y/PmmUFXIbN2ra3U2ndf+O9/w5bG8B6okxdMn27r7i66qPCVJ8A220C3bvDCC1kzkq2co44yJTpypPVEr7zSxmUffbTSHun8D97KAAAgAElEQVSQIbBqld2nmsC225pD9sGDzclCPrFkifWO/+//Cl95glkXX3ihrXedNy9xeScNRFJ6GXAFmic8/bQpzh49wpYkd/TqBb/8kkSos2whAkcfbV4DPvjAXu+vuKJSRfrUU+Yd6pBDQpA3JHr1gqVL8y/U2QsvWK+sV6+wJckdF15oHoqeeSZsSQoMkT8gMg34OtjfF5GEY2MJh3BFuB94RpWpmZAzHwl7CHfdOhvObN8e3ngjNDFyjqpNRzZvnidu41RtjvS22+Djj02wG26Anj2hXj2mTjUjrwcegD59whY2d5SXb3TvN2JE2NIYqhasfYstqs9a1UzRsaONWM2eHf5oVcEM4YqMAboBw1HdP8ibgmqlQSST6YF+AwwQYYwIl4jQuMrCOpvw9tvw44/2nK5JiNgb9UcfWRzt0BGxYIyjR8ftkQ564jdq17Z43zWJWrXsPr3/vlm65gNjx8KUKTWr9xmhVy9ztFIIPrLzCtXYgfH1iaokVKCqDFTlMOA8oCUwSYQXRTg6LSGdzXjqKTNMOfHEsCXJPT162Ft0XhnkVKBIr360DY/t9Sjbbpma1W4hcP759rXky9DhwIFmddu9e9iS5J5TToGtt7bvwMkY8xD5A6CI1EHkWiLDuZWQ1ByoCEXAHkFaDEwE+ojwchUEdoDvvzeLuh49CsMBdqrsuKO9ODz7LOSdT/8oRfrRLe8zq7wVvb66woyOHnkk5eUv1Zmdd4ZOnUyBrk/4Xp5dVq6El16Cs86CRo3ClSUM6taF886zOemffgpbmoLhEuByLIDJfGC/YL9SEipQER7EhnFPBO5S5UBV7lWlC7B/lUR2GDTI5pguvDBsScKjZ0+LPJMvpvmbIcJdY47l/1qMZv17H2y02q1hirRXL5g/P/x1iK+8Ykq0Jg7fRujZ02wnBg8OW5ICQKQIOBfVc1DdHtXtUP0/VJckrJqEEdGFwMuqbLbYQITGqixPW/A8ISwjovJyM6LZaac8MaIJiXXrrIdz8MH5aUT13Xc2invrrebMCFW7YbffvnEdacSzUYJ1pNWZtWvtt3rEESFaTgN/+IN5SZo6NX888oRBPnwPBWRENArVo1KtlswQ7jmxylOEDwAKQXmGyejRMGtWzTMeiqV2bRvCfvttWLgwbGk2JzLvd8EFQUb08peRI8089c9/NlPVhx+2CMgFSJ06NnQ4fDgsWhSODFOnWtSVXr1qtvIE+w6+/tq+D6fKfILII4gcgcgBG1ICKlSgItQToSmwjQhNRGgapJbAjhkUvMby1FM2h9O1a9iShM+FF9rc2qBBYUuyKevXm4HT8cebX9jNOPpo641++CH87nfmIrB1a3PbU4CKtGdPm6t+7rlwzj9wIDXSEjoeZ55pDiTcmCgj/AHYC+gHPBCk+xNVqnAIV4SrgKsxZRkdp/MX4ElVHqmiwHlDGEO4P/9slrehubLLQzp0MKOqGTPyp3dRUgKdO9u82xlnJFHho4/gb38zhbrDDhaX9OKLC8pR62GHmWOFadNye5/WrLHR8mOPtfvhQO/e5lBi4cJwDKoKZgg3TSrsgaryT1VaAdeq0ioq7VtIyjMsXnzRbE9q+vBtND17QmmpDW3nCwMH2pKBk09OskKHDjasO3q0uSzq02djYO888LecCXr2hG++sdjlueT1101x12TjoVh69TJXmC/7eoiqIXJr3JSoWiU90GNUGSnC6fGOq/JqFUXOG8LogR54oA0PTpiQP72tsFm92nrlp5wS3hBhND/9ZM6IrrgCHnwwzUY+/th6pO+/D9ttZ2FDLrsMGlTfl/aVK+0+nXFGbtfvHn+8OdyYPbuwoxWlgqo5l69XLxyPTAXTAxX5S9RePeAk4GtUK10fUdnPsEPw2SVOOil9SZ0JE+DLL90QIpb69S26xtChNsQdNoMHm4VwlUYJDj/c/N99/LE96a67Dlq2tJhoK1dmStSc0rChOTAYMsR8GeeCOXPsHeTCC115RiNiz5GxYy0mqpMmqg9EpTuBo7A1oYnqhR/VO+xUv359zSWXX65at67q0qU5PW21YOxYVVB97LFw5SgvV91zT9VDDslww59+qtqpk13k1lur3nWX6i+/ZPgk2efzz+0SnngiN+e7+WZVEdW5c3NzvurEkiX2PLnyytyfG1ilefAMz3iCJgozE5VLxpHCXSJsFbXfRIT+VdD1NZpff7VJ/9NPhyZNwpYm/zjwQHMS/tRT4crx+ee2RCDj822HHmrxqD77zBa+3nij9UjvvDN33bkMcPDB5lg/F/dp/XpbStS5swVdcDalaVN7ngweXJCG37lBZDIik4I0FZgO/DNRtWQGQ05QZcOAmirLMK9EVUaEziJMF6FUhL5xjtcVYUhwfEywhAYRWorwqwhfBenxqDoHijA5qPOwCHk1SPrf/9rw5IY1hc4mRBzMjx9vhiphMWSIuUw788wsneCQQ+zHMGaMKdWbbzZFescdsDz/l1eL2NrdL77IvoP5jz6yQPP+n6mYCy6w58q774YtSRxEOiMyHZFSRDZ7ziNSF5EhwfExiLSMOnZDkD8dkU5Jt5k6J7FxirIjsCOqiY1lE3VRQSeB1o3a3wJ0alW7yKBFoLNAW4PWAZ0I2jamzGWgjwfb3UGHBNstQadU0O4XoIeCCug7oCckkiWXQ7hnn20jd+vW5eyU1Y758214sH//cM5fXq7aooXqySfn8KRjx6p26WIXvtVWqrffrrpsWQ4FSJ3Zs03c++/P7nkuu0x1iy1UV67M7nmqM+vW2XPlnHNye14SDeFCkcIshdYKdRQmKrSNKXOZwuPBdneFIcF226B8XYVWQTtFSbWZ+pDtrgp1g+2jFP6ssFWiesn0QJ8HPhChZ+DWbwSQieXuBwOlqsxWZS3wMnBKTJlTos41FDi2sh6lCM2ARqp8Zt8DzwGnZkDWjLBmDbz1Fpx6as10HJ8szZtbBy0sd3Fjx8K8eTl2cNGunbn4GT/elsLcfrt5brj1Vlu7kYe0agX775/d+1ReDq++CiecUK0Nl7NOcbFZr7/5pj1n8oiDgVJUZ6Oa9HMeEQnyX0Z1DapzgNKgvWTaTJVhwHpE2gBPAa2AFxNVSiac2X1Af2BPoC1wR5BXVZoD0fHX5rO51dOGMqqUAcuBrYNjrUSYIMJHIhwRVX5+gjYBEJHeIjJORMaV5SgMyIgRsGKFex5Khq5dzVp59uzcn3vYMHsgdemS+3NzwAG24HHCBDjuOBvSbdnShniXJPRtnXO6drXp3AULstP+p59aoIFu3bLTfiHRrZtNo7//fk5PWxx5jgapd8zxlJ7zqEY/5yuqm0ybqVIenPt04B+oXgM0S1QpWYPwCcBHwKhgOxPE60nGLkqtqMxCYGdV9gf6AC+K0CjJNi1TdYCqtlPVdsU56g4OGwaNG5snFadyIi8Zr+Z4tbGq3adjjw3ZyGu//UyQSZPMeuauu0yR3nBDXsWwityn117LTvvDhpkP3j/+MTvtFxLHHmvPlxyP3JRFnqNBGhBzvCrP+VTzq8I6RM7G4l6/FeTVTlQpGSvcM4EvgG7AmcAYETLxPjgfiLap24lNXQZuUkaEYqAxsFSVNaosAVBlPDAL2D0ov1OCNkNh3TqLNNKliz0QnMrJxfBgPCZNMgf/eTNKsPfe5rdu8mTTIvfea1/OddeF59E9ij32MIdL2bhPkZeZjh1rZtzPVKlTx54vb7xhz5s8IaXnPCIbnvOV1E2mzVS5ADgUuBPVOYi0wqYvKyfRJGlg3LNd1P62oBOrNGFr7RSDzgZtFWVEtFdMmctjjIheiZKhKNhuDboAtGmwPxb0kCgjohMTyZILI6IRI8zg4rXXsn6qgqF/f/vO5s/P3TlvuUW1Vi3VH3/M3TlTYto0sxSpVcssa/r0UV24MFSRIt/ZokWZbXfMGLv/zz6b2XYLmddes+/svfdycz4SGxEVK8wOjIAiBj97xZS5PMaI6JVge68YI6LZgQFR4jarZlDURGGfZMomo+gmx+zXis1LN4GeCDojsMa9KcjrB3pysF0P9D+gpYF1besgvyvo1EDpfgnaJarNdqBTgjYfAZVEcuRCgV5yiWr9+qqrV2f9VAXD11/bL/Rf/8rdOdu2Ve3QIXfnS5tvvlE97zzTXPXqqV51leqCBaGIMmGC3acBAzLb7nXXqRYXu8ORVFi9WrVBA9WLL87N+RIqUFNIJyrMCCxnbwry+imcHGzXU/iPQqnCFwqto+reFNSbrnBCpW1WTWmOUmik0FRhrsJ4hQcT1UsmoPbfgX2Al4Kss4BJqlyfSv84n8m2L9z1682y9Igj4D//ydppCpK2bWH77S24Sbb55hvYc08L6Xnlldk/X0YoLbX50eeeM8uniy6C66+3yNc5QhXatLGwqJlah6hq7e26q0XEcZLnrLMswt7330NRUXbPVUC+cCeguj8ivYAWqN6GyCRU96msWjJWuH8FBmBKdF9gQCEpz1zw6afw4495NK9Wjeja1QKb5MJuJjKPd3rc8Al5Sps25tF9xgwLkvn446Z1LrsM5s7NiQgidp8++ACWLctMmxMn2ly0W9+mTteuNj3+8cdhS1KtKEakGWbn81aiwhGSssJVZZgqfVS5RpUs2dsVLsOGmVcbtyRMna5dbS3g669n/1zDhtn60+ZVNYgPg9at4cknLVzJBRdYHLY2bSwW6bffZv30XbtaoO0338xMe8OGmdP4U/NmFXf14cQTLTpLWOuoqyn9gBJgFqpjEWkNzExUqbJwZiuIbxos2NRpwdjFZXMIV9XWw++3n62Td1IjG8OD8Zg92zpuf/+7RRyr9sydC/fcY85qy8vN796NN5qizQLl5fY7P+AAswKtKm3bWjzykSOr3lZN5LTTzCHI3LnZjV5TMEO4aVJZQO0tVWkUJ21ZSMoz24Ti1aaAyMbwYDwi600L5j7tvDM89piNg156KTz/POy+u/VOZyZ8sU6ZWrVs6LukxJyFVIVp08yRf8HcixDo2tWcW4wZE7Yk1QSR3RH5AJEpwf4+iNycqFpS7yYiHC7CBcH2NiK0qpKwNYiIV5uTTw5bkupLpocH4/Hqq7butFWh/bJ32smsoubMMcuol1+2xZvnnQfTp2f0VF27mhu5d96pWjuRocfTTqu6TDWVLl2gdm0fxk2BJ4EbAFtBqzoJ6J6oUjKOFG4Drg8aB6hDMgtMHVTtwXzMMR66rCocdJDpgWx5JVqwwNzRFXSPp1kzeOghU6R9+tiTdc89LYL5tGkZOcVhh8F221X9oT1smLW1444ZEatG0rgxHH+8fZcJFlo4Rn1Uv4jJS+jjNZke6GnAycAqAFW+B7ZMWbwayOTJtsqgoB/MOSB6eHDlysy3H3FDVyPu0w472ETvnDnw17/axPzvf29rHyZPrlLTRUXWa3z77fTjUpaWmgVujbgXWaZrV7Mf+/LLsCWpFixGZFcidj8i3TCXsZWSjAJda+tMrWERauyEcaoMG2ZzeKdUNU6Aw+mnw2+/WQjNTDNsmBmt7LFH5tvOW7bbztwCfvut+dd95x2LZN6tm2mwNOnaFVatgvfeS69+tVxKlKeccoq91PgwblJcDjwB7IHIAuBq4JJElZJRoK+I8ASwlQgXAe9j48VOAoYNM+cJ228ftiTVn8MPz8zwYCw//WTrTGtsj2ebbeDOO02R3nqrhfLYbz9bPzJ+fMrNHXWUTVeke5+GDbMh+112Sa++s5Gtt4ajj4ahQ30Yt1JEagHtUD0O2BbYA9XDUf0uUdVkFGg58D8sXtruwK2q/Ksq8tYEpk+HqVNr8IM5wxQV2TO9KsOD8Xj9dVuCUePvU9Om8Le/mSK9/Xb46COLUdqlC3wROzVUMbVrW89n+HBYuzY1Eb77zqzWa/y9yCBdu5rR9ZQpYUuSx6iWA1cE26tQTdqOPBkFuiVmQHQI8C0wKXUJax4+FJV5qjo8GI9hw2z95z6VOuyqQWy1Fdx2m2mz/v3NjVb79hbR+rPPkmqia1dYvjz1NZwFt5QoDzj1VJtG8mHchIxA5FpEWiDSdENKQEJfuBsKCvtgfnC7AvNVOa5q8uYP2XCkcNBB1mv6/POMNlujWbfOhsO7dIFBgxKXT8SKFTbMdfXVcF8mQsQXIitW2HrS+++HxYstyPett9rcRAWsWWMjw+ecY54Fk6VDB/j55ypNwTpxOPJIe6HJxvdaMI4UROYQz3GQaqWeR1LxUbEI+AFYAmyXimw1jUWLYNw4e9A7maN2bYstXVJiw65VZeRIU8onnlj1tgqWLbc05/Rz5pj17qRJ9kQ+5hjzWB6HunUtuHNJSfJzb8uXW2fX3V1mnj/+0W7b93kRGTlvaQs8CkwEvgL+BeyVqFIy60AvFWEU8AGwDXCRKj7gVQmRIcZOncKVoxDp1Mkc82fibbqkBBo2hD/8oeptFTwNG5qPwzlz4MEHzVXQ0Udbt3HkyM00ZadONp2arNOjkSPNWYb/ZzJP5DvN5NRHATII2BN4GFOeewZ5lZJMD3QX4GpV9lLlNlUys+q6gCkpsSGsAw4IW5LCo2NH+8xEiKuSEtMBdepUva0aQ/36cM015jz4n/+0hZvHHmtDuu+9t0GRRh7ayd6nyMvMoYdmSe4azD772NSHh4WrlN+h2gvVD4PUG/hdokrJhDPrq8pXGRGxBlBebs+Rjh2z68S5ptKsGey7b9UfBqWlpgO8x5MmW2wBf/6z+dp99FEzOurUybrz77xD61ZKmzbJ3SdVK3fssf4ykw1q1bJbM2KExSZ24jIBkUM27Im0Bz5JVMkf8Rlm4kSbA/UHc/bo1Ak++aRqTssjD/bOnTMjU42lXj2LPVpaahZD339vk8rt2/OX373FhyOVNWsqb2LmTBvu9f9M9ujUCZYsSWtpb02hPfApIt8i8i3wGdABkcmIVLjyxBVohomE3IoMNTqZp1MnM/758MP023j3XVu+suuumZOrRlO3rsUenTnT4pIuXswlb3dh9K/t+ObeNyq1Jor8Z1yBZo/jj7flLD6MWyGdgVZAhyC1Ak4ETgIqNAcNVYGK0FmE6SKUitA3zvG6IgwJjo8RoWWQf7wI40WYHHweE1VnVNDmV0HKqcVwSYkNMe6wQy7PWrM47DCbikv3YbB2rSlff2BngTp1oFcvmD6d3x57msYsZ9/bTrVQN8OGxTWfLimxeK9ZClXqANtuCwce6Aq0QlS/qzRVQGgKVIQizGz4BMyE+GwR2sYU6wksU6UN8BBwb5C/GOiiyt5AD2BwTL1zVNkvSIuydhExrFhhQ4s+LJhd6ta1VRTpPgw++cQcMrgCzSK1a1Pv0gu4pMM33NRikLmP6tbN3AT+5z8bFOlvv9lqGL8X2adTJ1uX/vPPYUtSOITZAz0YKFVltiprgZeBWLfrp7DRlHgocKwIosqEICoMwFSgngh1cyJ1JXz4oZvi54pOncx+pbQ09bolJRaj9eijMy+XsynHn1DMXfPO4/v3p1lQ73Xr4MwzYe+94eWX+fij9axe7f+ZXNC5sxkRffBB2JIUDmEq0ObAvKj9+UFe3DKqlAHLga1jynQFJqgSbarwTDB8e4sIEu/kItJbRMaJyLiysoRh35Li3XehQQMbYnSyS6rLJKJ5911zTr+lB+XLOpHRmPc+KDLXRFOmwEsvWebZZ7P3n37PeUUvcNThmfkPOhXTvj00arRxztmpOmEq0HiKLdbSoNIyIuyFDeteHHX8nGBo94ggnRvv5Ko6QFXbqWq74uLilASvCF9XmDvatIFWrVJXoD/8YJbS3uPJDfvsY/YAG+5TURF0726xR195hV9WFTNo/f/R8OC28NxzNoTjZIXatVP3EOVUTpgKdD7QImp/JyDW2dSGMiIUA42BpcH+TsBrwHmqzIpUUGVB8LkCeBEbKs46kXWFPv+ZG0Tsux45MrWoH+4lKreImEX6ZmsQa9ViwR/O4HdrJvLaucPMKqxHDwvK+swzNtTrZJzOnWHePPjmm7AlKQzCVKBjgd1EaCVCHaA7MDymzHDMSAigGzBSFRVhK+Bt4AbVjYtdRSgWYZtguzZmgpyTQD5uip97OnUyY6BPEi533khJicUV3Xff7MnlbEpkDeKXX26a/957oNRi12tPt4Ovv25jjBdeCL/7HQwc6Io0w0SeTz6MmxlCU6DBnOYVQAnwNfCKKlNF6CfCyUGxp4CtRSgF+sCGpS5XAG2AW2KWq9QFSkSYhDkEXkCOgn+XlJgZfps2uTibAzZcXlyc/DCue4kKh4rWIJaUmGepvffGbsgpp9hK/zfftDA5F11k61sGDEg9uKgTl112sU5+XixnsZBhIxCZGXw2qaBcj6DMTER6ROUfGDg6KEXkYUQkyL8dkQWIfBWk7IWLUNUan+rXr69VYc0a1QYNVC+7rErNOGnQoYPqfvslV3bsWFVQff75rIrkxKFdO9XDDtu4X1am2qSJao8eFVQoL1f9739V27e3m9aihepjj6n+9lsuxC1orrpKtV491dWrq94WsErTffbCfQp9g+2+CvfGKdNUYXbw2STYbhIc+0LhUAVReEfhhCD/doVr05YrheTv4RnA1xWGR6dO8NVXZhyUiMhb9/HHZ1cmZ3MiaxCXL7f9ceNg2bJKbAZENgbxfvdd2Gkncxm4667wyCO2gNRJi06d7OsbPTpsSTZZpjgIODVOmU7ACFSXoroMGAF0RqQZ0AjVzzBF/lwF9bOKK9AM8O67vq4wLFIJ1VRSYhFytvNotjmnU6dN1yCWlJiOTPgyI7LR+fGIEWZ6feWVpkgfftgcNDgp0aGDOSPJ0DBucWQ5YJB6p1B3e1QXAgSf8f6ZFS13bB5sx+ZHuAKRSYg8XeHQcAZwBZoBSkps7aevK8w9++1nbsoSPQyWL7fOjI8ShMMhh9j/I3Kf3n0X2rWzqc6kEIHjjrNu08iRNjd61VVmePDQQ7B6ddZkLzTq17eY6BlSoGUaLAcM0oBNjoq8j8iUOCnWaU5FVLSUsbIljv8GdgX2AxYCDyR5rpRxBVpFIusKfflKOERCNb33Xlw3qxvwgM3hElmD+O67sHQpjBmT5n9GxIZ6Ro2y1LYt9OljPdP777e5FCchnTrBtGm2pCWrqB6H6u/jpDeAH4OhWILPeG5XK1ruOD/Yjs0H1R9RXY9qOWZEmrWljK5Aq4ivKwyfTp1g8WKYMKHiMiUl1gPygM3h0bkzzJ0Lgwfby06V/zMdOtiY8P/+Zx4b/vpXU6T33gsrV2ZE5kIl8vISsjVu9DLFHsAbccqUAB0RaRIMxXYESoIh3xWIHBJY3563oX5EKRunkcWljK5Aq8i77/q6wrCJhI6raG2bBgGbjznGvUSFSURhvvwyNG5sruUywuGH2/zoJ5/YJHffvtCyJdx9N/zyS4ZOUli0bQvNm4euQO8BjkdkJnB8sA8i7RAZCIDqUuAOzG/AWKBfkAdwKTAQKAVmAe8E+fdFxfE8GrgmWxcg6j6daNCgga5KY+invBy2397e5gbHxoNxcsoBB0DDhvEtC2fMsHX5jz0Gl16ae9mcjey+uw0b/vGPMHRolk4yZgz06wf//S80aWJDvFdeaVrb2UDPnhZhbvFiM4JMBxFZraoNMitZ9SEzTmBrKO+8Yz++Vq3g44839S9pi9fM8rC83D6jt2M/I9vJJNXNtyPni96Ol6Jli7edzGf0NcbbzjYSx3ygVi0byeva1ebbouWZMcM+33zT5kJjvyvYuA+bfoeRY+l+j/G+u8q+x1TKVvSdJ2q/srYqkzfR7yEeRUV2byJp6VJbQrFoEVx7bfw68e5vZfmbH2sPv3+b5luP47jP+tH2llv4tf8DfNLuKj5pdxW/1muSsL1USedeJPPbSfZ/Xdn/GOxaa9e2VFxsn0uXmnHdoEGmTJ3U8R4o6fdAO3cOfQgkJUQ2PjQq207mM7rNeNvZoqKf67p1sGZN/GNO+mTinsa7Z3Xrxu/1JKuIUqmzb/kEblzfj1PKX2c5jXis6M/8s9Y1LJOmlQueIskq/nj/mYr+gyL24hHZTibFtgGmgNetM0O66M9168y3fyRATurXXLN7oK5ASV+Bjh9vvdBDDtmYF/3niLx1R97Ci4o23Y53TGTzt/bIHyhyPDov+s8V749WU1zW/fQT/Pvftqylfn3Li74Xa9dafmXffUXH4t2Liu5H9H2o6MFXkykrMwcK22wTwncxcSL0729jxw0b2rBunz4mTA3l559tZDvde+EK1BVo2grUcZxqyJQppkhfecXeqi6/HP7yF/ewkQY1XYHWkP6J4zhOwO9/b6bAU6fCqafa+tFWrUyJJuMT0nECXIE6jlMz2XNPeP55+PprOOMM+Oc/TZFedRUsWBC2dE41wBWo4zg1m913h2efhenT4ZxzbL1T69Y2tDt3btjSOXmMK1DHcRwwB/UDB8LMmXDBBfDkkxbgt3dvmDMnbOmcPMQVqOM4TjQtW8Ljj8OsWaY8n3vOnNeff/7GRcWOgytQx3Gc+LRoYbFHZ8+2JS+vvGLzpuecY57YnRqPL2PBl7E4jpMEP/4IDzxgc6SrV8Of/mQ+PGvw4l5fxhIiInQWYboIpSL0jXO8rghDguNjRGgZdeyGIH+6CJ2SbdNxHCcttt8e7rsPvv0WbrwRdtyxRitPJ8QeqAhFwAzMC/98zNP+2apMiypzGbCPKpeI0B04TZWzRGgLvITFedsReB/YPahWaZvx8B6o4zhO6ngPNDwOBkpVma3KWuBlIDZK+SnAoGB7KHCsCBLkv6zKGlXmYOFsDk6yTcdxHMepMmEq0OZAdDz0+UFe3DKqlAHLga0rqZtMmwCISG8RGSci48rKyqpwGY7jOE5NJEwFGm/yIHY8uaIyqeZvnqk6QFXbqWq74nSD4TmO4zg1ljAV6HygRdT+TsD3FZURoRhoDCytpG4ybTqO4zhOlQlTgY4FdhOhlQh1gO7A8Jgyw4EewXY3YKQqGnGMIhIAAAsrSURBVOR3D6x0WwG7AV8k2abjOI7jVJnQxi5VKRPhCqAEKAKeVmWqCP2AcaoMB54CBotQivU8uwd1p4rwCjANKAMuV2U9QLw2c31tjuM4TuHjjhTwZSyO4zjp4MtYHMdxHMdJGVegjuM4jpMGrkAdx3EcJw1cgTqO4zhOGrgCdRzHcZw0cAXqOI7j5B6RpoiMQGRm8NmkgnI9gjIzEekRlX8nIvMQWRlTvi4iQxApRWQMIi2zdQmuQB3HcZww6At8gOpuwAfB/qaINAVuA9pjwUJui1K0bwZ5sfQElqHaBngIuDfzohuuQB3HcZwwiI62NQg4NU6ZTsAIVJeiugwYAXQGQPVzVBcmaHcocCySncCtrkAdx3GcdCmORLUKUu8U6m6/QQHa53ZxyiQdYStuHdXoKF4Zx8OQOI7jOOlSpqrtKjwq8j6wQ5wjNyXZftIRtqpYJy1cgTqO4zjZQfW4Co+J/IhIM1QXItIMWBSn1HzgqKj9nYBRCc4aico1H5HoKF4Zx4dwHcdxnDCIjrbVA3gjTpkSoCMiTQLjoY5BXrLtdgNGkiWn765AHcdxnDC4BzgekZnA8cE+iLRDZCAAqkuBO7BQlWOBfkEeiNyHyHygPiLzEbk9aPcpYGtESoE+xLPuzRAejQWPxuI4jpMOHo3FcRzHcZyUcQXqOI7jOGngCtRxHMdx0sAVqOM4juOkQSgKVISmIowQYWbwGdeJsAg9gjIzRcwsWYT6IrwtwjciTBUJLLfs2Pki/CTCV0HqlatrchzHcWoWYfVA+wIfqFKhE2ERNnMiHKVo71dlD2B/4DARToiqOkSV/YI0MKtX4TiO49RYwlKgSTsRVmWpKhucCKuyWpUPAVRZC3yJeadwHMdxnJwRlgLdXpWFAMFnWk6ERdgK6IL1YiN0FWGSCENFaJFZsR3HcRzHyJovXBGy6kRYhGLgJeBhVWYH2W8CL6myRoRLsN7tMfHlk95Ab4A6deokKZLjOI7jGFlToKpU6ERYhB9FaKbKQhHSdSI8AJipyj+izrkk6viTVBJIVVUHBG3QoEEDd8fkOI7jpERYQ7hJOxEWoUlgPLTBibAI/TEP+1dHVwiUcYSTga8zLLfjOI7jACH5whVha+AVYGdgLnCGKktFaAdcomrLT0S4ELgxqHanKs+IsBM2N/oNsCY49ogqA0W4G1OcZVj4mktV+SaRPO4L13EcJ3Vqui9cdyaPK1DHcZx0qOkK1D0ROY7jOE4auAJ1HMdxnDRwBeo4juM4aeAK1HEcx3HSwBWo4ziO46SBK1DHcRzHSQNXoI7jOI6TBq5AHcdxHCcNXIE6juM4Thq4AnUcx3GcNHAF6jiO4zhp4ArUcRzHcdLAFajjOI6Te0SaIjICkZnBZ5MKyvUIysxEpEdU/p2IzENkZUz58xH5CZGvgtQrW5fgCtRxHMcJg77AB6juBnwQ7G+KSFPgNqA9cDBwW5SifTPIi8cQVPcL0sCMSx7gCtRxHMcJg1OAQcH2IODUOGU6ASNQXYrqMmAE0BkA1c9RXZgLQSvCFajjOI6TLsUiMi4q9U6h7vYbFKB9bhenTHNgXtT+/CAvEV0RmYTIUERapCBTShRnq2HHcRyn4ClT1XYVHhV5H9ghzpGbkmxf4uRpgjpvAi+hugaRS7De7TFJni8lXIE6juM42UH1uAqPifyISDNUFyLSDFgUp9R84Kio/Z2AUQnOuSRq70ng3iSlTRkfwnUcx3HCYDgQsartAbwRp0wJ0BGRJoHxUMcgr2JMGUc4Gfi66qLGJxQFKkJTEUaIMDP4jGu+LEKPoMxMkQ1fNCKMEmG6CF8Fabsgv64IQ0QoFWGMCC1zdEmO4zhOatwDHI/ITOD4YB9E2iFilrOqS4E7gLFB6hfkgch9iMwH6iMyH5Hbg3b/jMhURCYCfwbOz9YFiGqi4eQsnFS4D1iqyj0i9AWaqHJ9TJmmwDigHTbmPR44UJVlIowCrlVlXEydy4B9VLlEhO7AaaqclUieBg0a6KpVqzJybY7jODUFEVmtqg3CliMswhrCTdp8WZWlqmxqvpxcu0OBY0XiTkI7juM4TpUIy4hoe1UWAqiyMDIEG0Mi8+VnRFgPDAP6q6LRdVQpE2E5sDWwOLbxwNw6YnKtIvJrApmLgbKEV1a98GvKfwrtesCvqbqQzDVtkQtB8pWsKVARsmm+fI4qC0TYElOg5wLPJaizaabqAGBAkrIgIuMqNdeuhvg15T+Fdj3g11RdKMRryjRZU6CqVGi+LMKPIjQLep8pmy+rsiD4XCHCi5g7p+eCOi2A+SIUA42BpVW/GsdxHMfZlLDmQJM2XxahSWCl2xEoEaFYhG0ARKgNnARMidNuN2BkMLTrOI7jOBklrDnQe4BXROgJzAXOABChHXCJKr1UWSqywXwZoF+Q1wBTpLX5//buH0SOMozj+PfHcYmFgsqBCho1kMI/TU4IkYCksJAUWmhxlUSwUUQtg4WClVqkEIuIJIUiMSCipySFomCVQwmJSQzqJZWcqERIPJTI6WMx74Vz3bl9b7I7s2/u94Fh527m4Hl4lnt3Z2efByaAz6i+LAuwH3hHYp7qnefMEGPOvtxbEOc0/q62fMA5leJqzGmoOvkai5mZWencicjMzKwBL6BmZmYNeAGtIelGSZ9K+iE91rQb1N+Sjqdttu04B5H0kKTvJM1L+t/AWkkbJR1Kx+ck3dF+lGuTkdNuSb+uqMvIJtIPi6QDkn6RdKrmuCS9nnL+RtJ02zGuVUZOOyVdWFGnF9uOcS0k3SbpC0lnJJ2W9Fyfc4qqU2ZORdWpVRHhrc8GvAbsSft7gFdrzlvsOtZVcpgAzgKbgQ3ACeDunnOeBval/RngUNdxDyGn3cAbXce6xrweAKaBUzXHdwFHqL7rvB2Y6zrmIeS0E/ik6zjXkM8twHTavw74vs9zr6g6ZeZUVJ3a3PwOtF5Ou8Fxtw2Yj4hzEfEX8B5VXiv1aX+ocW5/mJNTcSLiS1b/zvIjwNtROQpcr/9OnRg7GTkVJSJ+iohjaf93qikfvcOdi6pTZk5WwwtovZsiTUuP+mnpANekSexHJY3bIpszzX1F+8NYgsvtD8dV7oT6R9MltPc1won0LcrNuzT3Szoh6Yike7oOJlf6qGMrMNdzqNg6rZITFFqnUVvXA7V15dPSATZFxIKkzcDnkk5GxNnhRHjFclobNpn43qWceD8GDkbEJY14In2LSqtTjmPA7RGxKGkX8CGwpeOYBpJ0LVUL0ecj4mLv4T5/MvZ1GpBTkXVqw7p+BxoRD0bEvX22j4Cfly+9qH5aOhGxkB7PUbUa3NpS+DmWWxsuuxVYqDtHUgntDwfmFBHnI+JS+vEt4L6WYhulnFoWJSIuRsRi2j8MTEqa6jisVUmapFpo3o2ID/qcUlydBuVUYp3asq4X0AEGthuUdIOkjWl/CtgBfNtahIN9BWyRdKekDVQ3CfXeKdyn/WGM8yvmgTmpxYn0LZoFHk93eW4HLix/xFAqSTcvf94uaRvV/6Pz3UZVL8W6HzgTEXtrTiuqTjk5lVanNq3rS7gDpHaD6mk3qNRuMJ4E7gLelPQP1ZPqlYgYmwU0IpYkPUPVV3gCOBARpyW9DHwdEbNcbn+oUbQ/HLrMnJ6V9DDVKKbfGOFE+mGRdJDqbscpST8CLwGTABGxDzhMdYfnPPAH8EQ3kebLyOkx4ClJS8CfwMyYv3jbQTX56aSk4+l3LwCboNg65eRUWp1a41Z+ZmZmDfgSrpmZWQNeQM3MzBrwAmpmZtaAF1AzM7MGvICamZk14AXUzMysAS+gZmZmDfwLXKP8adM7SfUAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ "Exact pressure gradient : -8.000e-03\n", "Numerical pressure gradient: -7.074e-03\n" ] } ], "source": [ "X, Y, LA = sp.symbols('X, Y, LA')\n", "rho, qx, qy = sp.symbols('rho, qx, qy')\n", "\n", "def bc(f, m, x, y):\n", " m[qx] = rhoo * vmax * (1.-4.*y**2/W**2)\n", " m[qy] = 0.\n", "\n", "def plot_coupe(sol):\n", " fig, ax1 = plt.subplots()\n", " ax2 = ax1.twinx()\n", " ax1.cla()\n", " ax2.cla()\n", " mx = int(sol.domain.shape_in[0]/2)\n", " my = int(sol.domain.shape_in[1]/2)\n", " x = sol.domain.x\n", " y = sol.domain.y\n", " u = sol.m[qx] / rhoo\n", " for i in [0,mx,-1]:\n", " ax1.plot(y+x[i], u[i, :], 'b')\n", " for j in [0,my,-1]:\n", " ax1.plot(x+y[j], u[:,j], 'b')\n", " ax1.set_ylabel('velocity', color='b')\n", " for tl in ax1.get_yticklabels():\n", " tl.set_color('b')\n", " ax1.set_ylim(-.5*rhoo*vmax, 1.5*rhoo*vmax)\n", " p = sol.m[rho][:,my] * la**2 / 3.0\n", " p -= np.average(p)\n", " ax2.plot(x, p, 'r')\n", " ax2.set_ylabel('pressure', color='r')\n", " for tl in ax2.get_yticklabels():\n", " tl.set_color('r')\n", " ax2.set_ylim(pressure_gradient*L, -pressure_gradient*L)\n", " plt.title('Poiseuille flow at t = {0:f}'.format(sol.t))\n", " plt.draw()\n", " plt.pause(1.e-3)\n", "\n", "# parameters\n", "dx = 1./16 # spatial step\n", "la = 1. # velocity of the scheme\n", "Tf = 50 # final time of the simulation\n", "L = 2 # length of the domain\n", "W = 1 # width of the domain\n", "vmax = 0.1 # maximal velocity obtained in the middle of the channel\n", "rhoo = 1. # mean value of the density\n", "mu = 1.e-2 # bulk viscosity\n", "eta = 1.e-2 # shear viscosity\n", "pressure_gradient = - vmax * 8.0 / W**2 * eta\n", "# initialization\n", "xmin, xmax, ymin, ymax = 0.0, L, -0.5*W, 0.5*W\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", "dico = {\n", " 'box': {'x': [xmin, xmax],\n", " 'y': [ymin, ymax],\n", " 'label': 0\n", " },\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}\n", " },\n", " 'generator': 'cython',\n", "}\n", "\n", "sol = pylbm.Simulation(dico)\n", "while (sol.t