{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Transport in 1D\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 test the most simple lattice Boltzmann scheme $\\DdQq{1}{2}$ on two classical hyperbolic scalar equations: the advection equation and the Burger's equation.\n", "\n", "## The advection equation\n", "\n", "The problem reads\n", "$$\\drondt u + c\\drondx u = 0, \\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." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The numerical simulation of this equation by a lattice Boltzmann scheme consists in the approximation of the solution on discret points of $(0,1)$ at discret instants.\n", "\n", "The spatial mesh is defined by using a numpy array. To simplify, the mesh is supposed to be uniform.\n", "\n", "First, we import the package numpy and we create the spatial mesh. One phantom cell has to be added at each edge of the domain for the treatment of the boundary conditions." ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "%matplotlib inline " ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYYAAAD8CAYAAABzTgP2AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAADvJJREFUeJzt3H+s3XV9x/HnSzrYjI5frS2j1MtGzVY1meYENfvFBkIxkZJJFliMdWFr4saSybasxmQwNJlsUxYzNleFrDOZwkw2b8JcgyAxMcI4VeesDlsRpQhSKSMhRBn63h/ny3Y/d6fcc+85nNN7fT6Sm57v+X7a7/uTW/LsOd9zSVUhSdKzXjDrASRJxxfDIElqGAZJUsMwSJIahkGS1DAMkqSGYZAkNQyDJKlhGCRJjXWzHmAl1q9fX3Nzc7MeQ5JWlf3793+nqjYstW5VhmFubo5+vz/rMSRpVUnyjVHW+VaSJKlhGCRJDcMgSWoYBklSwzBIkhqGQZLUMAySpIZhkCQ1DIMkqWEYJEkNwyBJahgGSVLDMEiSGoZBktQwDJKkhmGQJDUMgySpYRgkSQ3DIElqGAZJUsMwSJIahkGS1DAMkqSGYZAkNSYShiTbk9yX5FCS3UPOn5Tklu78PUnmFp3fkuTJJH8wiXkkSSs3dhiSnADcCFwMbAOuSLJt0bIrgcer6hzgBuD6ReffB3xi3FkkSeObxCuGc4FDVXV/VT0NfBTYsWjNDmBv9/hjwPlJApDkUuDrwIEJzCJJGtMkwnAm8OCC48Pdc0PXVNUzwBPA6UleBPwR8CcTmEOSNAGzvvl8LXBDVT251MIku5L0k/SPHDny/E8mST+k1k3gz3gIOGvB8ebuuWFrDidZB5wMPAa8BrgsyZ8BpwA/SPLdqvqrxRepqj3AHoBer1cTmFuSNMQkwnAvsDXJ2QwCcDnw64vWzAM7gc8ClwF3VlUBv/DsgiTXAk8Oi4IkaXrGDkNVPZPkKmAfcAJwc1UdSHId0K+qeeAm4MNJDgFHGcRDknQcyuAf7qtLr9erfr8/6zEkaVVJsr+qekutm/XNZ0nSccYwSJIahkGS1DAMkqSGYZAkNQyDJKlhGCRJDcMgSWoYBklSwzBIkhqGQZLUMAySpIZhkCQ1DIMkqWEYJEkNwyBJahgGSVLDMEiSGoZBktQwDJKkhmGQJDUMgySpYRgkSQ3DIElqGAZJUsMwSJIahkGS1DAMkqSGYZAkNQyDJKlhGCRJjYmEIcn2JPclOZRk95DzJyW5pTt/T5K57vnXJ9mf5D+6X39lEvNIklZu7DAkOQG4EbgY2AZckWTbomVXAo9X1TnADcD13fPfAd5YVa8EdgIfHnceSdJ4JvGK4VzgUFXdX1VPAx8FdixaswPY2z3+GHB+klTV56vqW93zB4AfS3LSBGaSJK3QJMJwJvDgguPD3XND11TVM8ATwOmL1rwJ+FxVfW8CM0mSVmjdrAcASPJyBm8vXfgca3YBuwC2bNkypckk6YfPJF4xPAScteB4c/fc0DVJ1gEnA491x5uBfwLeUlVfO9ZFqmpPVfWqqrdhw4YJjC1JGmYSYbgX2Jrk7CQnApcD84vWzDO4uQxwGXBnVVWSU4DbgN1V9ZkJzCJJGtPYYejuGVwF7AO+AtxaVQeSXJfkkm7ZTcDpSQ4BVwPPfqT1KuAc4I+TfKH7esm4M0mSVi5VNesZlq3X61W/35/1GJK0qiTZX1W9pdb5k8+SpIZhkCQ1DIMkqWEYJEkNwyBJahgGSVLDMEiSGoZBktQwDJKkhmGQJDUMgySpYRgkSQ3DIElqGAZJUsMwSJIahkGS1DAMkqSGYZAkNQyDJKlhGCRJDcMgSWoYBklSwzBIkhqGQZLUMAySpIZhkCQ1DIMkqWEYJEkNwyBJahgGSVJjImFIsj3JfUkOJdk95PxJSW7pzt+TZG7BuXd0z9+X5KJJzCNJWrmxw5DkBOBG4GJgG3BFkm2Lll0JPF5V5wA3ANd3v3cbcDnwcmA78NfdnzdRmzZtIsn/+9q0adOkLzX1663lvU37emt5b9O+3lre27SvN+29AaSqxvsDktcB11bVRd3xOwCq6k8XrNnXrflsknXAI8AGYPfCtQvXPdc1e71e9fv95cx4zHPj7n/W11vLe5v29dby3qZ9vbW8t2lfb5LXSrK/qnpLrZvEW0lnAg8uOD7cPTd0TVU9AzwBnD7i75UkTdGqufmcZFeSfpL+kSNHZj2OJK1ZkwjDQ8BZC443d88NXdO9lXQy8NiIvxeAqtpTVb2q6m3YsGECY0uShplEGO4FtiY5O8mJDG4mzy9aMw/s7B5fBtxZgzfH5oHLu08tnQ1sBf5tAjNJklZo7DB09wyuAvYBXwFuraoDSa5Lckm37Cbg9CSHgKv5v5vOB4BbgS8D/wr8TlV9f9yZFtu4ceOynl9N11vLe5v29dby3qZ9vbW8t2lfb9p7gwl8KmkWlvupJEnSdD+VJElaQwyDJKlhGCRJDcMgSWoYBklSwzBIkhqGQZLUMAySpIZhkCQ1DIMkqWEYJEkNwyBJahgGSVLDMEiSGoZBktQwDJKkhmGQJDUMgySpYRgkSQ3DIElqGAZJUsMwSJIahkGS1DAMkqSGYZAkNQyDJKlhGCRJDcMgSWoYBklSwzBIkhqGQZLUGCsMSU5LcnuSg92vpx5j3c5uzcEkO7vnXpjktiT/meRAkveMM4skaTLGfcWwG7ijqrYCd3THjSSnAdcArwHOBa5ZEJC/qKqfBl4F/FySi8ecR5I0pnHDsAPY2z3eC1w6ZM1FwO1VdbSqHgduB7ZX1VNV9SmAqnoa+Bywecx5JEljGjcMG6vq4e7xI8DGIWvOBB5ccHy4e+5/JTkFeCODVx2SpBlat9SCJJ8ENg059c6FB1VVSWq5AyRZB3wEeH9V3f8c63YBuwC2bNmy3MtIkka0ZBiq6oJjnUvy7SRnVNXDSc4AHh2y7CHgvAXHm4G7FhzvAQ5W1V8uMceebi29Xm/ZAZIkjWbct5LmgZ3d453Ax4es2QdcmOTU7qbzhd1zJHk3cDLwe2POIUmakHHD8B7g9UkOAhd0xyTpJfkQQFUdBd4F3Nt9XVdVR5NsZvB21Dbgc0m+kOQ3x5xHkjSmVK2+d2V6vV71+/1ZjyFJq0qS/VXVW2qdP/ksSWoYBklSwzBIkhqGQZLUMAySpIZhkCQ1DIMkqWEYJEkNwyBJahgGSVLDMEiSGoZBktQwDJKkhmGQJDUMgySpYRgkSQ3DIElqGAZJUsMwSJIahkGS1DAMkqSGYZAkNQyDJKlhGCRJDcMgSWoYBklSwzBIkhqGQZLUMAySpIZhkCQ1xgpDktOS3J7kYPfrqcdYt7NbczDJziHn55N8aZxZJEmTMe4rht3AHVW1FbijO24kOQ24BngNcC5wzcKAJPlV4Mkx55AkTci4YdgB7O0e7wUuHbLmIuD2qjpaVY8DtwPbAZK8CLgaePeYc0iSJmTcMGysqoe7x48AG4esORN4cMHx4e45gHcB7wWeGnMOSdKErFtqQZJPApuGnHrnwoOqqiQ16oWT/CzwU1X19iRzI6zfBewC2LJly6iXkSQt05JhqKoLjnUuybeTnFFVDyc5A3h0yLKHgPMWHG8G7gJeB/SSPNDN8ZIkd1XVeQxRVXuAPQC9Xm/kAEmSlmfct5LmgWc/ZbQT+PiQNfuAC5Oc2t10vhDYV1V/U1U/UVVzwM8DXz1WFCRJ0zNuGN4DvD7JQeCC7pgkvSQfAqiqowzuJdzbfV3XPSdJOg6lavW9K9Pr9arf7896DElaVZLsr6reUuv8yWdJUsMwSJIahkGS1DAMkqSGYZAkNQyDJKlhGCRJDcMgSWoYBklSwzBIkhqGQZLUMAySpIZhkCQ1DIMkqWEYJEkNwyBJahgGSVLDMEiSGoZBktQwDJKkhmGQJDUMgySpYRgkSQ3DIElqGAZJUiNVNesZli3JEeAbU77seuA7U77mtLi31Wst78+9Td5Lq2rDUotWZRhmIUm/qnqznuP54N5Wr7W8P/c2O76VJElqGAZJUsMwjG7PrAd4Hrm31Wst78+9zYj3GCRJDV8xSJIahmGRJNuT3JfkUJLdQ86flOSW7vw9SeamP+XKjLC3q5N8OckXk9yR5KWzmHMlltrbgnVvSlJJjttPhCw2yt6S/Fr3vTuQ5B+mPeM4Rvh7uSXJp5J8vvu7+YZZzLkSSW5O8miSLx3jfJK8v9v7F5O8etozDlVVfnVfwAnA14CfBE4E/h3YtmjNbwMf6B5fDtwy67knuLdfBl7YPX7bWtpbt+7FwKeBu4HerOee4PdtK/B54NTu+CWznnvC+9sDvK17vA14YNZzL2N/vwi8GvjSMc6/AfgEEOC1wD2znrmqfMWwyLnAoaq6v6qeBj4K7Fi0Zgewt3v8MeD8JJnijCu15N6q6lNV9VR3eDewecozrtQo3zeAdwHXA9+d5nBjGmVvvwXcWFWPA1TVo1OecRyj7K+AH+8enwx8a4rzjaWqPg0cfY4lO4C/r4G7gVOSnDGd6Y7NMLTOBB5ccHy4e27omqp6BngCOH0q041nlL0tdCWDf8msBkvurXuJflZV3TbNwSZglO/by4CXJflMkruTbJ/adOMbZX/XAm9Ochj4F+B3pzPaVCz3v8upWDfrAXT8SfJmoAf80qxnmYQkLwDeB7x1xqM8X9YxeDvpPAav8j6d5JVV9V8znWpyrgD+rqrem+R1wIeTvKKqfjDrwdYqXzG0HgLOWnC8uXtu6Jok6xi8tH1sKtONZ5S9keQC4J3AJVX1vSnNNq6l9vZi4BXAXUkeYPBe7vwquQE9yvftMDBfVf9dVV8HvsogFKvBKPu7ErgVoKo+C/wog//X0Fow0n+X02YYWvcCW5OcneREBjeX5xetmQd2do8vA+6s7i7ScW7JvSV5FfC3DKKwmt6nfs69VdUTVbW+quaqao7B/ZNLqqo/m3GXZZS/k//M4NUCSdYzeGvp/mkOOYZR9vdN4HyAJD/DIAxHpjrl82ceeEv36aTXAk9U1cOzHsq3khaoqmeSXAXsY/BpiZur6kCS64B+Vc0DNzF4KXuIwU2ly2c38ehG3NufAy8C/rG7n/7NqrpkZkOPaMS9rUoj7m0fcGGSLwPfB/6wqlbDq9hR9/f7wAeTvJ3Bjei3rpJ/jJHkIwyivb67R3IN8CMAVfUBBvdM3gAcAp4CfmM2k7b8yWdJUsO3kiRJDcMgSWoYBklSwzBIkhqGQZLUMAySpIZhkCQ1DIMkqfE/3LkV4hS13xoAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "import numpy as np\n", "import pylab as plt\n", "\n", "def mesh(N):\n", " xmin, xmax = 0., 1.\n", " dx = 1./N\n", " x = np.linspace(xmin-.5*dx, xmax+.5*dx, N+2)\n", " return x\n", "\n", "x = mesh(10)\n", "plt.plot(x, 0.*x, 'sk')\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "To simulate this equation, we use the $\\DdQq{1}{2}$ scheme given by\n", "\n", "* two velocities $v_0=-1$, $v_1=1$, with associated distribution functions $\\fk{0}$ and $\\fk{1}$,\n", "* a space step $\\dx$ and a time step $\\dt$, the ration $\\lambda=\\dx/\\dt$ is called the scheme velocity,\n", "* two moments $\\mk{0}=\\sum_{i=0}^1\\fk{i}$ and $\\mk{1}=\\lambda \\sum_{i=0}^1 v_i\\fk{i}$ and their equilibrium values $\\mke{0} = \\mk{0}$, $\\mke{1} = c\\mk{0}$,\n", "* a relaxation parameter $s$ lying in $[0,2]$.\n", "\n", "In order to prepare the formalism of the package pylbm, we introduce the two polynomials that define the moments: $P_0 = 1$ and $P_1=\\lambda X$, such that\n", "$$ \n", "\\mk{k} = \\sum_{i=0}^1 P_k(v_i) \\fk{i}.\n", "$$\n", "\n", "The transformation $(\\fk{0}, \\fk{1})\\mapsto(\\mk{0},\\mk{1})$ is invertible if, and only if, the polynomials $(P_0,P_1)$ is a free set over the stencil of velocities.\n", "\n", "The lattice Boltzmann method consists to compute the distribution functions $\\fk{0}$ and $\\fk{1}$ in each point of the lattice $x$ and at each time $t^n=n\\dt$.\n", "A step of the scheme can be read as a splitting between the relaxation phase and the transport phase:\n", "\n", "* relaxation: $$\\mks{1}(t,x)=(1-s)\\mk{1}(t,x)+s\\mke{1}(t,x).$$\n", "* m2f: \n", "$$\\begin{aligned}\\fks{0}(t,x)&\\;=(\\mk{0}(t,x)-\\mks{1}(t,x)/\\lambda)/2, \\\\ \\fks{1}(t,x)&\\;=(\\mk{0}(t,x)+\\mks{1}(t,x)/\\lambda)/2.\\end{aligned}\n", "$$\n", "* transport: $$\\fk{0}(t+\\dt, x)=\\fks{0}(t,x+\\dx), \\qquad \\fk{1}(t+\\dt, x)=\\fks{1}(t,x-\\dx).$$\n", "* f2m: \n", "$$\\begin{aligned}\\mk{0}(t+\\dt,x)&\\;=\\fk{0}(t+\\dt,x)+\\fk{1}(t+\\dt,x), \\\\ \\mk{1}(t+\\dt,x)&\\;=-\\lambda\\fk{0}(t+\\dt,x)+\\lambda\\fk{1}(t+\\dt,x).\\end{aligned}\n", "$$\n", "\n", "The moment of order $0$, $\\mk{0}$, being the only one conserved during the relaxation phase, the equivalent equation of this scheme reads at first order\n", "\n", "$$\\drondt\\mk{0} + \\drondx\\mke{1} = \\grandO(\\dt).$$\n", "\n", "We implement a function equilibrium that computes the equilibrium value $\\mke{1}$, the moment of order $0$, $\\mk{0}$, and the velocity $c$ being given in argument." ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "def equilibrium(m0, c):\n", " return c*m0\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Then, we create two vectors $\\mk{0}$ and $\\mk{1}$ with shape the shape of the mesh and initialize them. The moment of order $0$ should contain the initial value of the unknown $u$ and the moment of order $1$ the corresponding equilibrium value.\n", "\n", "We create also two vectors $\\fk{0}$ and $\\fk{1}$." ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "def initialize(mesh, c, la):\n", " m0 = np.zeros(mesh.shape)\n", " m0[np.logical_and(mesh<0.5, mesh>0.25)] = 1.\n", " m1 = equilibrium(m0, c)\n", " f0, f1 = np.empty(m0.shape), np.empty(m0.shape)\n", " m2f(f0, f1, m0, m1, la)\n", " return f0, f1, m0, m1\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "And finally, we implement the four elementary functions f2m, relaxation, m2f, and transport. In the transport function, the boundary conditions should be implemented: we will use periodic conditions by copying the informations in the phantom cells." ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [], "source": [ "def f2m(f0, f1, m0, m1, la):\n", " m0[:] = f0 + f1\n", " m1[:] = la*(f1 - f0)\n", " \n", "def m2f(f0, f1, m0, m1, la):\n", " f0[:] = 0.5*(m0-m1/la)\n", " f1[:] = 0.5*(m0+m1/la)\n", "\n", "def relaxation(m0, m1, c, s):\n", " m1[:] = (1-s)*m1 + s*equilibrium(m0, c)\n", "\n", "def transport(f0, f1):\n", " #periodical boundary conditions\n", " f0[-1] = f0[1]\n", " f1[0] = f1[-2]\n", " #transport\n", " f0[1:-1] = f0[2:]\n", " f1[1:-1] = f1[:-2]\n", " " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We compute and we plot the numerical solution at time $T_f=2$." ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAEICAYAAACktLTqAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3Xt4VfWd7/H3NzuBhDsSQO63YpGiRURh1LH0aKvSKWjVipeq1ZbqHKe1zpk5ts5Yjz7OIzNTtY5SRx1rneOl1g6Wo45ap1qrU9Ao2haUigEhgCThfk+y8zt//PYOm7CTvfZ975XP63ny7NvK2t+VwCe//V1r/ZY55xARkXCpKHYBIiKSewp3EZEQUriLiISQwl1EJIQU7iIiIaRwFxEJIYW7hIqZXWlmrxf4Pf/czFYX8j1FUlG4S9kws1fNbLuZ9S5yHc7MPhV/7Jz7rXPu08WsSaQzhbuUBTMbD/w54IB5RS1GpAwo3KVcXA4sAx4Brog/aWZDzGypme0yszeBSQmv/djM/jlxJWb2SzO7IXZ/pJn9wsyazGytmX07YbmImX3fzD4ys91m9raZjTGz12KLvGdme8zsIjObY2YNCd97bOxTxg4zW2lm8xJee8TM7jOz52LrXW5mHTWL5IrCXcrF5cBjsa+zzGx47Pn7gAPACOCq2FfcE8BFZmYAZjYY+CLwpJlVAP8PeA8YBZwBXG9mZ8W+9wbgYmAuMCC23n3OudNjr3/WOdfPOfezxCLNrCq23peAYcBfAY+ZWWLbZgHwf4DBwBrg9kx/KCJdUbhLyTOz04BxwFPOubeBj4BLzCwCnA/c7Jzb65z7I/DThG/9Lb6N8+exxxcAv3PObQJOAoY65251zrU45+qBB/HBC/AN4O+cc6ud955zbmuAcmcD/YA7Yuv9NfAs/g9F3BLn3JvOuTb8H6vp6f5MRFJRuEs5uAJ4yTnXHHv8eOy5oUAlsCFh2Y/jd5yfFe9JDgXrJfgwBf/HYmSsdbLDzHYA3wfinwjG4P+IpGsksME5196pplEJjz9JuL8P/8dAJKcqi12ASHfMrAb4KhAxs3go9gYG4YO4DR/EH8ReG9tpFU8AL5nZHcAs4LzY8xuAtc65yV289QZ8//6PaZa8CRhjZhUJAT8W+FOa6xHJikbuUurOBaLAVHz7YjpwLL7lcjnwH8AtZtbHzKaSsLMVwDm3AmgGHgJedM7tiL30JrDbzP63mdXEdqBOM7OTYq8/BNxmZpPNO97MhsRe2wJM7KLe5fjR+N+aWZWZzQG+jP8EIVIwCncpdVcAP3HOrXfOfRL/Au4FLgWuw7c1PsEfSfOTJOt4HDgzdguAcy4K/AX+j8VaDv0BGBhb5E7gKfyO0V3AvwE1sdduAX4aa+d8NfGNnHMt+DA/J7bOxcDlzrkPECkg08U6RETCRyN3EZEQUriLiISQwl1EJIQU7iIiIVS049xra2vd+PHji/X2IiJl6e233252zg1NtVzRwn38+PHU1dUV6+1FRMqSmX2ceim1ZUREQknhLiISQgp3EZEQ0sRhIlJ2WltbaWho4MCBA8UuJW+qq6sZPXo0VVVVGX2/wl1Eyk5DQwP9+/dn/PjxxK7FEirOObZu3UpDQwMTJkzIaB1qy4hI2Tlw4ABDhgwJZbADmBlDhgzJ6pOJwl1EylJYgz0u2+1TW0ZKm3OwZAm8+65/PGYMXH01VGhcItId/Q+R0rVvH3z963D++XDbbf5r4UK45ZZiVybCKaecknKZb3zjG6xatQqAf/iHf8h3SYdRuEtpWbIEvvIV/zV9Ojz6KPzgB9DWBu3tPuxvuw0efzz1ukTy6L//+79TLvPQQw8xdepUoATD3cweNrNGM0t6LcnYJcjuMbM1ZvZ7M5uR+zKlR3j+ebjwQnjzTVizBgYPhhde8CP1SATM4P774fTT4aqr4KKLYMECWLy42JVLD9Svn7+u+auvvsqcOXO44IILmDJlCpdeeinxiyDNmTOHuro6brzxRvbv38/06dO59NJLC1JfkJ77I/hLmj3axevnAJNjX7OAH8dupUS0t7ezc+fOYpfRrcjbb9P/wguJTpvG7qVLoX//Qy9u337YsvZv/0bfhQupWLEC27oVe+45dnz1q/4PQJb69OlD7969s16PFM7111/Pu/F9Mjkyffp07r777sDLr1ixgpUrVzJy5EhOPfVU3njjDU477bSO1++44w7uvffenNfZnZTh7px7zczGd7PIfOBR5/9ULTOzQWY2wjm3OUc1SpYuu+wynnjiiWKXcZiBwD8BVwLxUzTqgVPee48t48YFXs+lwP8FTh86lKQfLdM0YsQINmzYQCQHfyik5zj55JMZPXo04P8wrFu37rBwL4ZcHC0zCtiQ8Lgh9twR4W5mC4GFAGPHjs3BW0sQ9fX1fPrTn+baa68tWg2DGhsZu3o1kbY2Im1tfPa3v6XPnj2snDWLfQMG0F5RwfsnncT3jjoqvfU2NcHtt3PXRRex8s/+LKsaX375ZZ599lkOHjxInz59slqXFE46I+x8Sfy0F4lEaGtrK2I1XkEPhXTOPQA8ADBz5kxdmbtAotEokyZN4jvf+U5h3nDbNr/Tc+NG//j99+GPncbV06fDQw9x3IkndjyVUTQ7Bz/+MWf278+ZWW5fa2srzz77LNFoNKv1iCRTVVVFa2trxtMJpCsX4b4RGJPweHTsOSkR0Wi0cG2GF1/0R7Q0NcHkyf654cPhRz+CefMgPjLv39/vIM2WGcyaBcuWZb2q+M9I4S75sHDhQo4//nhmzJjBY489lvf3y0W4LwWuM7Mn8TtSd6rfXlryGu7OwXnnwS9/eei5qVPhuefghBPy856dzZrlj6rZvfvwHbFpUrhLOvbs2QP4I2LmzJnT8fy9997bcf/VV1/tuL9o0SIWLVpUqPJSh7uZPQHMAWrNrAH4AbF9YM65+4HngbnAGmAf8PV8FSuZaWtry1+4P/GED/YrroBx46C2Fr75Taiuzs/7JTNrlv8jU1cHn/98xquJ/4xKoV8qkq0gR8tcnOJ1B/zPnFUkOReNRqmszMPulb174W//Fk48ER5+uHhTApx8sr9dtiyrcI//jDRylzDQ3DI9QN7aMnfc4Xea/uxnxZ3r5aij4JhjYPnyrFajtoyEicK9B8hpuD/8MMRPo/74Y7jkEjj11NysOxuzZsGvfuXbMxnuqFW4S5hobpkeIKc997vu8vO8zJ7tJ/G6667crDdbs2fDJ5/ANdfAt7+d0dEz6rlLmGjk3gPkrOdeX++PV7/zTvjud7NfXy6ddRaMHg0//zns2gXr18Mzz6S1CvXcJUw0cu8BctaWiR/uOH9+9uvKtUmTYMMGfwLV2Wf7cE+T2jKSjnvuuYdjjz2WwYMHc8cdd2S8nvgEZLmmkXsPkNNwnzYNJk7Mfl35NHZsVm0ZhbsEsXjxYl5++eWOOWVKjUbuPUBOeu5bt8JvfwvnnpubovJp7Fhf7969aX2beu4S1DXXXEN9fT3nnHMOd911F9dddx0AV155Jd/+9rc55ZRTmDhxIk8//TTgT3g644wzmDFjBscddxy/TDzpL080cu8BctJzf+45f7GMUmzJdBaflG7DBpgyJfC3qedepq6//tBlGHNl+nToZkKy+++/nxdeeIFXXnmFZ5999rDXNm/ezOuvv84HH3zAvHnzuOCCC6iurmbJkiUMGDCA5uZmZs+ezbx58/J6HViFew+QVVtm82ZobfU7KkeN8icslbp4uK9fn1a4qy0juXDuuedSUVHB1KlT2bJlCwDOOb7//e/z2muvUVFRwcaNG9myZQtHH3103upQuPcAGYf700/7KyPF/eVf5mayr3xLDPc0KNzLVAlM+Zsocfrf+BWZHnvsMZqamnj77bepqqpi/PjxHDhwIK91KNx7gIx77r/5DfTtC/fc489A/Yu/yH1x+TBypK83w3BXz11ybefOnQwbNoyqqipeeeUVPv7447y/p8K9B8i4575ihe89XnVV7ovKp8pK30JKM9zVc5d8ufTSS/nyl7/Mcccdx8yZM5mSRrswUwr3kGtvbwdIf+Te3u53UpVbsMeNHau2jOTVunXrAH+EzJVXXgnAI488ctgy8WmBa2tr+d3vfpd0PfFlck2HQoZcPKjSDvcPP/SHEhZqTvZcU7hLD6dwD7l4/zjtcF+xwt+Wc7hv2OA/gQSknruEicI95OKj0LR77itWQK9e/qpK5WjsWGhpgcbGwN+innt5iR+JElbZbp/CPeQybsu8846faqBXrzxUVQAZHA6ptkz5qK6uZuvWraENeOccW7dupTqLK5pph2rIZRTuzvmRezlMNdCVxHCPX6kpBYV7+Rg9ejQNDQ00NTUVu5S8qa6uzmreGoV7yGXUc29o8HOzzJiRp6oKIIuRu3rupa+qqooJEyYUu4ySprZMyGXUc3/nHX9brjtTAQYOhP790wp39dwlTBTuIZdRW2bFCj/NwPHH56mqAjBL+3BItWUkTNSWCbnA4b5vH1x7LezY4UfuU6b4qQfKmcJdejCN3EMucM/9zTfh0Ufh/fehttZfi7TcjR0L69b5gN+0KeXi6rlLmGjkHnKBe+5r1/rbF14o/SstBTVxot8xPG6cf/zoo/C1r3W5uHruEiYK95AL3Japr/czKY4ZU4CqCuRb34Kjj4a2Nli4EFav7nZxtWUkTNSWCbm0wn3sWKiqKkBVBTJwIFx+uZ/8bOjQlGerKtwlTBTuIRe45752bXjaMckMGxY43NVzlzBQuIdc4J57fT2E+aSQAOGunruESaBwN7OzzWy1ma0xsxuTvD7WzF4xsxVm9nszm5v7UiUTgdoye/fCli3hH7nHrmfZFbVlJExShruZRYD7gHOAqcDFZtZ5qsC/A55yzp0ALAAW57pQyUygcI9ddCDU4T58uHru0qMEGbmfDKxxztU751qAJ4H5nZZxwIDY/YFA6oOKpSAC9dzr6/1t2Nsye/b4k7W6oJ67hEmQcB8FbEh43BB7LtEtwGVm1gA8D/xVshWZ2UIzqzOzujDP5lZKAvXc4+Ee5pH7sGH+tpt/dxUV/r+DRu4SBrnaoXox8IhzbjQwF/h3Mzti3c65B5xzM51zM4cOHZqjt5buBGrLrF0L/fr5M1PDKh7u3bRmzIxIJKJwl1AIEu4bgcQzW0bHnkt0NfAUgHPud0A1EOKkKB+Bwj1+pIxZgaoqggDhDijcJTSChPtbwGQzm2BmvfA7TJd2WmY9cAaAmR2LD3f1XUpA4J57mFsykFa4q+cuYZAy3J1zbcB1wIvA+/ijYlaa2a1mNi+22F8D3zSz94AngCtdWK9/VWZS9tydC/8JTODPUIVAx7pr5C5hEGhuGefc8/gdpYnP3ZxwfxVwam5Lk1xI2ZZpbPRHkIT5SBnw0xf37au2jPQYOkM15FKGe084UiYu4BQECncJA4V7yKXsucen+g37yB0Ch7t67hIGmvI35Lrsuf/iF/7CHK+/7h+PH1/Ywoph2LCUV2ZSz13CQuEecknbMm1tsGCBvwU46STo06cI1RXY8OFQV9ftImrLSFioLRNyScO9sdEH+733QmsrLF9epOoKbNgwf4Zqe3uXiyjcJSwU7iGXtOe+ebO/HTUKKivDffJSomHD/B+1HTu6XEQ9dwkLhXvIJe25x8N95MgiVFREAU5kUs9dwkLhHnJJ2zLxcB8xoggVFVGAcFdbRsJC4R5y3Yb78OFFqKiIFO7SgyjcQy5pz33TJj8DZK9eRaqqSAKGu3ruEgYK95Drsufe01oyAEOG+J3H6rlLD6BwD7ku2zI9MdwrK33Aqy0jPYDCPeS6DPeedqRMXIoLZSvcJSwU7iF3RM+9vR0++aRnjtzBh/uqVfDoo7BkyREnNKnnLmGhcA+5+Cg0fn1Qtm71J/L01HCfMgU++ACuuAK+8pUjzs5Vz13CQuEectFolEgkgsXPQt20yd/21HD/l3+Bjz6C557zj+OHhcaoLSNhoYnDQi4e7h166glMcZWVfu763r394+bmw15WuEtYaOQecm1tbQr3ZIYM8bdJwl09dwkDhXvIRaPR5PPK9PRwr66Gfv2OCHf13CUsFO4hl7QtM2gQ1NQUr6hSUVurtoyElsI95I4I902bNGqPGzJE4S6hpXAPuaQ9d4W718XIXT13CQOFe8gl7bkr3L0k4a6eu4SFwj3kDmvLONezpx7orLbWn9SVQG0ZCQuFe8gdFu47dsDBgxq5x9XWwq5d0NLS8ZTCXcJC4R5yh/XcdRjk4Wpr/W3C6F09dwkLhXvIRaNRekcifibEVav8kwp3Lx7uCX139dwlLDT9QMhFo1Hu2rwZjj760JOjRxevoFKSJNzVlpGwULiHXDQaZXJLC5x4Ilx9NQwdCpMmFbus0qBwlxAL1JYxs7PNbLWZrTGzG7tY5qtmtsrMVprZ47ktUzLV1tbG4LY2OPlkuPZauOCCYpdUOroId/XcJQxSjtzNLALcB3wBaADeMrOlzrlVCctMBr4HnOqc225mw/JVsKTHtbYyKBr1I3Y5XJLJw9Rzl7AIMnI/GVjjnKt3zrUATwLzOy3zTeA+59x2AOdc1xeplILqe+CAv6NwP1JVFQwcqLaMhFKQcB8FbEh43BB7LtExwDFm9oaZLTOzs5OtyMwWmlmdmdU1NTVlVrGkpb/CvXudzlJVuEtY5OpQyEpgMjAHuBh40MwGdV7IOfeAc26mc27mUIVNQfTbv9/f0c87uSTh3t7ejnOuiEWJZC9IuG8ExiQ8Hh17LlEDsNQ51+qcWwv8CR/2UmQDDh70dxTuyXUK9/g8PBq9S7kLEu5vAZPNbIKZ9QIWAEs7LfMMftSOmdXi2zT1OaxTMqRwT6HTtL/xs3kV7lLuUoa7c64NuA54EXgfeMo5t9LMbjWzebHFXgS2mtkq4BXgb5xzW5OvUQppQHzelPiRIXK4JG0ZULhL+Qt0EpNz7nng+U7P3Zxw3wE3xL6khAxsaWF3VRX9q6qKXUppqq2Ffftg/36oqekIdx3rLuVOc8uE3MDWVnb16lXsMkpXp8nD1HOXsFC4h9zg1lZ29+5d7DJKV6ezVNWWkbBQuIfcoLY2dlVXF7uM0qVwl5BSuIfc4LY2divcu9ZFuKvnLuVO4R5m7e0MjkbZW1NT7EpKV6dwV89dwkLhHmY7dlAJ7FG4d23wYDBTW0ZCR/O5h1ls/p49ffoUuZASVlnpA37JEmhu5jPbtwMKdyl/Cvcwi4X7PoV79846C156CX7yE07at4/+qOcu5U9tmTCLh3vfvkUupMQ9/rhvyyxeDMBQNHKX8qdwDzOFe3pi8+8o3CUMFO5hFgv3/f36FbmQMqFwlxBRuIdZUxO7AaczVINJCHf13KXcKdzDrKmJZrOOY7clBY3cJUQU7mHW1EQTh47dlhT69iXau7fCXUJB4R5mTU00OadwT0PLwIEKdwkFhXuIuaYmGtHIPR1tgwap5y6hoHAPK+c62jLquQfXOniwRu4SCgr3sNq9G2tpUc89TVGFu4SEwj1snnkGPvtZmDULQOGeprajjlK4Sygo3MPm2Wfhww/h2GNpveACfoXCPR3Ro46iL+D27i12KSJZUTM2bJqbYfJk+I//YO+OHWx6+mn13NPghgwBILJtW5ErEcmORu5h09zccQGKeGtBI/fg2mPhXrVjR5ErEcmOwj1sFO5ZcbGfncJdyp3CPWyamyE2+owfq61wDy4e7pUKdylzCvcwiUZh27YjRu7quadh2DAAeincpcwp3MNkxw5/8pLaMhmrGDiQFqDX7t3FLkUkKwr3MIld5FnhnrlIZSXNQO9du4pdikhWFO5h0inc1XNPXyQSoQmFu5Q/hXuYdDFyV889uMrKSpqAarVlpMwFCnczO9vMVpvZGjO7sZvlzjczZ2Yzc1eiBKa2TNbiI/cahbuUuZThbmYR4D7gHGAqcLGZTU2yXH/gO8DyXBcpASncsxYP9+o9e4pdikhWgozcTwbWOOfqnXMtwJPA/CTL3QYsAg7ksD5JR3Mz1NRAnz6Aeu6Z6Aj3AwegpaXY5YhkLEi4jwI2JDxuiD3XwcxmAGOcc891tyIzW2hmdWZW19TUlHaxkkLCCUygnnsm4j134NAnIZEylPUOVTOrAO4E/jrVss65B5xzM51zM4fGLkYsOZQw9QCoLZOJ+MgdAA1ApIwFCfeNwJiEx6Njz8X1B6YBr5rZOmA2sFQ7VYtg61aFe5YU7hIWQT6vvwVMNrMJ+FBfAFwSf9E5txPoSBQzexX4X865utyWKik1N8O4cR0P1XNPX0VFxaFwf+EF2L0bjj8eJk0qZlkiaUs5cnfOtQHXAS8C7wNPOedWmtmtZjYv3wVKGrpoy6jnnp4tkQhtkQj88Ifwla/AhRcWuySRtAXquTvnnnfOHeOcm+Scuz323M3OuaVJlp2jUXsRtLXB9u1qy+TAvspKFi1cCO++CwsWwLp1xS5JJG06QzUs4lcOUrhnLRKJsL1PH38t2s98xv/RPKAjfKW8KNzDotMJTKCee6YikUjHz44RI/ztJ58UryCRDCjcwyJJuKvnnpnKysqOn11HuG/eXLyCRDKgcA+LeLgnOYlJI/f0RCIRhbuUPYV7WHQzcle4p0fhLmGgcA+LrVv9bcLIXT33zBzWcx86FCoqFO5SdhTuYdHcDH37+onDYtRzz8xhPfdIBI4+WuEuZUfhHhadTmACtWUydVhbBnxrRuEuZUbhHhYK95xRuEsYKNzDIkm4q+eemcN67qBwl7KkcC937e0QjXY7clfPPT2H9dzBh3tjo5/iQaRMKNzL2euvQ3U1VFZCfb0/siOB2jKZSdqWcQ62bCleUSJp0pCunK1YAa2t8L3v+SNlLrnksJcV7plJGu7gWzOjRiX/JpESo3AvZ42N/hjs227zh+x1op57ZpL23EF9dykrasuUs8ZG32fvIrzVc89M0p47KNylrCjcy1ljIwwb1uXLastk5oi2zPDh/lbhLmVE4V7OFO55cUS49+rlPyEp3KWMKNzLWYpwj/eNKyr0a07HET130LHuUnb0v76cBRi5q9+eviN67qBwl7KjcC9XBw7Arl2H+sFJRKNRtWQycERbBmDkSIW7lBWFe7lqbPS3KUbuCvf0JQ33ESP8pfba24tTlEiaFO7lKkC4t7W1Kdwz0GXPva3t0Lz5IiVO4V6uAo7c1XNPX5c9d4Dp02H8eHjooYLXJZIOhXu5Ulsmb5K2Zc44A669Fr7wBdi7F375y+IUJxKQhnXlSuGeN0nDffBgWLzY3z/vPPjTnwpfmEgaNHIvV42N/pJ6fft2uYh67plJ2nNPNGECrF3rZ4oUKVEK93LV2OgPgzTrchH13DOTtOeeaOJE2L9fUwBLSVO4l6sUJzCB2jKZStqWSTRxor+try9MQSIZULiXK4V73qQM9wkT/O3atYUpSCQDgcLdzM42s9VmtsbMbkzy+g1mtsrMfm9m/2Vm43JfqhwmQLir556ZlD338eP9rUbuUsJShruZRYD7gHOAqcDFZja102IrgJnOueOBp4F/zHWhksC5wCN39dzTl7LnXlPjpyNQuEsJCzJyPxlY45yrd861AE8C8xMXcM694pzbF3u4DBid2zLlMDt2+MvrqS2TFynbMnDoiBmREhUk3EcBGxIeN8Se68rVwH8me8HMFppZnZnVNTU1Ba9SDhfgGHdQuGcqULhPnKiRu5S0nO5QNbPLgJnAPyV73Tn3gHNupnNu5tChQ3P51j1LPNy7mRES1HPPVMqeO/hwb2iAlpbCFCWSpiDhvhEYk/B4dOy5w5jZmcBNwDzn3MHclCdJpTFyV889fSl77uDbMs7Bxx8XpiiRNAUJ97eAyWY2wcx6AQuApYkLmNkJwL/ig70x92XKYdSWyavAbRlQa0ZKVspwd861AdcBLwLvA08551aa2a1mNi+22D8B/YCfm9m7Zra0i9VJLsTDvba228UU7pmJRCI452jvbu72eLhrp6qUqECf2Z1zzwPPd3ru5oT7Z+a4LulOYyMMGQIpWi5tbW1UVVUVqKjwiP9BjEajXV9/dsQI6N1bI3cpWWrIlpMVK2DbNli5MmVLBnw41dTUFKCwcInvp4hGo13/cayo8CczKdylRCncy8XatTBjxqHH55yT8lvUlslM4si9WzrWXUqYwr1cfPihv128GKZNg6mdTxI+ksI9M4HDfeJE+M1v4IoroKoKbrrp0LwzIkWmcC8X69f727lzYVywqXt0nHtm4j+zlMe6f+lL8OKL8NprsG4djB0LN9/c/feIFIhmhSwX69f7Pu/IkYG/Rce5Zyax596tuXNhzRrfmpkyBd55pwDViQSjcC8X69f7YE/j6Be1ZTITuC2T6IQT/A5vkRKhcC8X69f7j/1pULhnJqNwnzHD/46am/NUlUh6FO7lIoNwV889M4F77oniRzJp9C4lQuFeDtrbYcOGjEbu6rmnL3DPPdEJJ/hb9d2lRCjcy0Fjo599UG2ZgsioLTN4sD8MUuEuJULhXg7ih0Eq3Asio3AHP3pXuEuJULiXgwzDXT33zGTUcwffd1+zBnbuzENVIulRuJeDLEbu6rmnL6OeOxzaqfreezmuSCR9CvdysH499OsHgwal9W1qy2Qm47ZMPNzVmpESoGFdOYgfBmmW1rcp3DOTcbgPH+5PNHvuOT9j5IAB8PnPp/17E8kFjdzLQQbHuIN67pnKuOcOcMop8PLLcN55cMYZft4ZkSJQuJeDDMNdPffMZNxzB3jkEX8i05tvQnU1/OIXuS1OJCCFe6nbvx+amjIOd43c05dxWwagb1+YPh1OOgm++EV45hl/IW2RAlO4l7oNG/xtmuEevwaowj19WYV7ovPO87+/t9/OQVUi6VG4l7osDoMEFO4ZyKrnnujLX4ZIBJYsyUFVIulRuJe6LMNdPff0ZdVzTzRkCJx+usJdikLhXqq+/nX47Gf9pdvMYNSotL5dI/fM5awtA7418/77sHp19usSSYPCvRQ1NfmjLioqYPZs+Pu/h1690lqFwj1zOQ33c8/1t+efD1/4Anz3u36WT5E802f2UrR8ub/90Y/8x/oMxPvFCvf05aznDjBmDFx/vT80cts2uPtufzTNFVdkv242AXmBAAAJCklEQVSRbijcS9Hy5X5H3IknZrwK9dwzl7Oee9xdd/nb9nb/Seymm+DCC6FPn9ysXyQJtWVK0bJlcNxx/pjpDKktk7mctmUSVVTAnXfCxo3wwx/mdt0inWhYV2ra2/1H+Isvzmo1CvfM5S3cAU47zfffFy2CmppD+1VOOSX37yU9msK91KxeDbt2waxZWa1GPffM5bTnnsyiRX7Omb/5m/gbwsMPw+WX5+f9pEcK1JYxs7PNbLWZrTGzG5O83tvMfhZ7fbmZjc91oT3GsmX+dvbsrFajnnvmct5z72zSJGho8Bf1+OQT+Nzn/A7WRYv8/pbly2HPnvy8t/QYKf/nm1kEuA/4AtAAvGVmS51zqxIWuxrY7pz7lJktABYBF+Wj4JLgnL+u6b59/vGgQf4amrmwfDkMHAif/nRWq1FbJnN5bcvE9erlvwYM8FMEX3wx3Jgwbho4EL71LT+aj7dvRo8G/bGWgIL8SzkZWOOcqwcwsyeB+UBiuM8Hbondfxq418zMudzPmLT2V79i869/feiJxLmyE+4f9sZJlolWVtLavz+tffrgKirAOSIHD1K1bx+VBw4AUNHWxqD6empXraJm69aO52q2bSPS0tKxyvZIhIZTT+WjL32JPSNGANAa/w+Zpv/x0kscHDeON555Ju3vTbRp0yZA4Z6J+M+srq6O2tragrynXXIJQ044gUhLCxWtrYx5/XVG/fM/Y//4jx3LtPXqxc4JE9hXWwtmuIoKWvr35+CAAUR798bFnqOiouO+M/OPE59L9e8yxfzzLtP56fM5r32e1p2vKd9GzZ3L2M99Lk9r94KE+yhgQ8LjBqBzQ7hjGedcm5ntBIYAzYkLmdlCYCHA2AxmOQT4+O67mfP88xl9bybagT8Cf8L/oqP4Df0Y2B1bZno0ytWvvcbnE+bufgf4IrA1jffqA+wE/gW4+fzzs64dKFg4hUm/fv3o3bs3Dz74IA8++GDR6hgPnAYY0Av4TEsLM1evZnjsbNcqYBgwoFgFSsZe27ChJMI9Z5xzDwAPAMycOTOjP4rTFi3iw699Lb7CxJWnvJ/4t90OHCCyaxeR3bs7lnE1NUT796e9psaPTsw4OGkSDBjAMQnfOzVJXet372bgyy8T2buXiv37mf7AA2yYOJGPHnyQ9gHd//c76qmnGPDqq0T27aNyxQouu+ce5ufgF19TU8OnPvWprNfT0/Tr14/6+nqam5tTL1wEBxJudwN28CDW0oK1t/ujrdrbO+6bcxCN+tv48921m1J92O7m9W7Hzvmc9jhf685jzdOmT8/buuOChPtGYEzC49Gx55It02BmlcBA0hu0BlY7bRq106blY9XZO/XUQ/fnzqVm/nymXXedf97s0NV5Ej3+ONx+OxxzjO/dn3UWE668Evr3L2jpcriRI0cycuTIYpchkrEgTeG3gMlmNsHMegELgKWdllkKxM+nvgD4dT767WXlnHPgqadg+3Z/+9Ofwpln+sPf4v365cvhqqv8FAN/+IN//MILCnYRyZoFyWAzmwvcDUSAh51zt5vZrUCdc26pmVUD/w6cAGwDFsR3wHZl5syZrq6uLusNKBv798MNN8D99/uLKA8a5C/kUFvrT1pSb1xEAjCzt51zM1MuV6wBdo8L97hnnoEnnvC90epqP8/IlCnFrkpEykTQcNdBs4V27rmHpoEVEckTTRwmIhJCCncRkRBSuIuIhJDCXUQkhBTuIiIhpHAXEQkhhbuISAgp3EVEQqhoZ6iaWRN+5tx01dJpKuEeoCduM/TM7e6J2ww9c7sz3eZxzrmhqRYqWrhnyszqgpx6GyY9cZuhZ253T9xm6Jnbne9tVltGRCSEFO4iIiFUjuH+QLELKIKeuM3QM7e7J24z9Mztzus2l13PXUREUivHkbuIiKSgcBcRCaGSDXczO9vMVpvZGjO7Mcnrvc3sZ7HXl5vZ+MJXmVsBtvkGM1tlZr83s/8ys3HFqDPXUm13wnLnm5kzs7I/ZC7INpvZV2O/75Vm9niha8y1AP++x5rZK2a2IvZvfG4x6swlM3vYzBrN7I9dvG5mdk/sZ/J7M5uRszd3zpXcF/5arR8BE4FewHvA1E7L/CVwf+z+AuBnxa67ANv8eaBP7P615b7NQbc7tlx/4DVgGTCz2HUX4Hc9GVgBDI49HlbsuguwzQ8A18buTwXWFbvuHGz36cAM4I9dvD4X+E/AgNnA8ly9d6mO3E8G1jjn6p1zLcCTwPxOy8wHfhq7/zRwhplZAWvMtZTb7Jx7xTm3L/ZwGTC6wDXmQ5DfNcBtwCLgQCGLy5Mg2/xN4D7n3HYA51xjgWvMtSDb7IABsfsDgU0FrC8vnHOvAdu6WWQ+8KjzlgGDzGxELt67VMN9FLAh4XFD7Lmkyzjn2oCdwJCCVJcfQbY50dX4v/jlLuV2xz6qjnHOPVfIwvIoyO/6GOAYM3vDzJaZ2dkFqy4/gmzzLcBlZtYAPA/8VWFKK6p0/98HpgtklyEzuwyYCXyu2LXkm5lVAHcCVxa5lEKrxLdm5uA/ob1mZsc553YUtar8uhh4xDn3QzP7M+DfzWyac6692IWVo1IduW8ExiQ8Hh17LukyZlaJ/xi3tSDV5UeQbcbMzgRuAuY55w4WqLZ8SrXd/YFpwKtmtg7fl1xa5jtVg/yuG4ClzrlW59xa4E/4sC9XQbb5auApAOfc74Bq/ORaYRbo/30mSjXc3wImm9kEM+uF32G6tNMyS4ErYvcvAH7tYnsoylTKbTazE4B/xQd7ufdg47rdbufcTudcrXNuvHNuPH5fwzznXF1xys2JIP++n8GP2jGzWnybpr6QReZYkG1eD5wBYGbH4sO9qaBVFt5S4PLYUTOzgZ3Ouc05WXOx9yZ3s5d5Ln608hFwU+y5W/H/scH/4n8OrAHeBCYWu+YCbPPLwBbg3djX0mLXXIjt7rTsq5T50TIBf9eGb0etAv4ALCh2zQXY5qnAG/gjad4FvljsmnOwzU8Am4FW/Kexq4FrgGsSfs/3xX4mf8jlv21NPyAiEkKl2pYREZEsKNxFREJI4S4iEkIKdxGREFK4i4iEkMJdRCSEFO4iIiH0/wEiI11ElJrfUQAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# parameters\n", "c = .5 # velocity for the transport equation\n", "Tf = 2. # final time\n", "N = 128 # number of points in space\n", "la = 1. # scheme velocity\n", "s = 1.8 # relaxation parameter\n", "# initialization\n", "x = mesh(N)\n", "f0, f1, m0, m1 = initialize(x, c, la)\n", "t = 0\n", "dt = (x[1]-x[0])/la\n", "plt.figure(1)\n", "plt.clf()\n", "plt.plot(x[1:-1], m0[1:-1], 'k', label='init')\n", "while t0, \\quad x\\in(0, 1).$$\n", "\n", "The previous $\\DdQq{1}{2}$ scheme can simulate the Burger's equation by modifying the equilibrium value of the moment of order $1$ $\\mke{1}$. \n", "It now reads $\\mke{1} = {\\mk{0}}^2/2$.\n", "\n", "More generaly, the simulated equation is into the conservative form\n", "$$\\drondt u + \\drondx \\varphi(u) = 0, \\quad t>0, \\quad x\\in(0, 1),$$\n", "\n", "the equilibrium has to be taken to $\\mke{1}=\\varphi(\\mk{0})$.\n", "\n", "We just have to modify the equilibrium and the initialization of the previous example to simulate the Burger's equation. The initial condition can be a discontinuous function in order to simulate Riemann problems. Note that the function f2m, m2f, relaxation, and transport are unchanged." ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX8AAAEICAYAAAC3Y/QeAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3XmcVNWZ//HP001Dqywim5oWgYijrSCaFo3JiDOgUaPAREQUNSQqamKMySQ/icYlxjE6/kYT4wYxRjHKIi81JEFxVBSNMQhqjLgFFbQRWRoEiSzd9Jk/TlUo2m66qu6turdufd+vFy/qVt3luSxPnX7OueeYcw4RESkvFVEHICIixafkLyJShpT8RUTKkJK/iEgZUvIXESlDSv4iImVIyV+kxJnZYjM7Juo4pLSYxvlLsZjZUqAPsA1oBJ4HLnDOfRBlXKXEzO4B6p1zP446FiltavlLsZ3snOsM7AWsBH6Zz0nMrEOYQYV9PpG4U/KXSDjnNgOzgNr0e2b2tJmdm7E9wcyey9h2ZvZtM/s78PfUe8eZ2Vtmtt7MbjezZ1qc45tm9oaZrTOzuWa2b1vnM+9mM1tlZhvM7G9mdnBr8ZtZNzP7tZmtMLPlZnatmVWmPqs0s/9vZmvM7N3UNVz6C8bMlprZiIxzXW1mv83YftDMPkrd03wzOyj1/kRgPPD/zGyjmf2+5fnMrJOZ/dzMPkz9+rmZdUp9doyZ1ZvZf6bucYWZfSP3vz1JAiV/iYSZ7QqcBryQ46GjgSOAWjPrif8C+RHQA3gLOCrjGqOAy4CvAb2AZ4FpbZ0POA44Gtgf6AaMBRraiOMeoAnYDzg0dWz6S+c84KTU+3XAmBzv8VFgINAbeAm4H8A5NyX1+r+dc52dcye3cuzlwJHAEOAQYCiQWSLaM3VvnwPOAW4zs+45xicJoOQvxfaImX0MrAeOBW7M8fifOefWOuc2AScCi51zDznnmoBbgI8y9r0gtf8bqc+vA4Zktv5bnK8R6AIcgO8Pe8M5t6JlAGbWJ3XtS5xz/3DOrQJuBsaldhkL/Nw594Fzbi3ws1xu0Dl3t3PuE+fcFuBq4BAz65bl4eOBa5xzq5xzq4GfAGdlfN6Y+rzROTcH2Aj8Sy7xSTIo+UuxjXbO7Q5UAxcBz5jZnjkcn9k5vHfmtvOjF+ozPt8X+IWZfZz6wlkLGL7V+5nzOeeeAm4FbgNWmdkUM+vaSgz7AlXAioxzT8a31D8TF7As25tLlYyuN7N3zGwDsDT1Uc8sT7F3i+stS72X1pD6Ikz7FOicbXySHEr+Egnn3Dbn3EP4kT9fTr39D2DXjN1a+1LIHJ62AqhJb5iZZW7jE/D5zrndM37t4px7vo3z4Zy7xTn3BXwZaH/gh63E8AGwBeiZcd6uzrmDMuLaJ2P/vi2O39l9ngGMAkbgyzP90rfXWryt+BD/5ZR57Q/bOUbKkJK/RCLVuToK6A68kXr7FeBrZrarme2Hr0nvzB+BQWY2OtWZ+m12TKR3Aj/K6DDtZman7iSmw83sCDOrwifozUBzy/1SpaDHgf8xs65mVmFmnzezYaldZgIXm1lNqp4+qcUpXgHGmVmVmbXsE+iC/2JpwH9BXNfi2JXAgDb/RHyfxo/NrFeqT+RK4Lc72V/KlJK/FNvvzWwjsAH4L+DrzrnFqc9uBrbiE9y9pDo62+KcWwOcCvw3PlnWAgvxyRPn3MPADcD0VAnlNeCEnZyyK/ArYB2+XNJA230SZwMdgddT+8/CD18ldY65wF/xHbYPtTj2CuDzqeN+AjyQ8dnU1LWXp87dskP81/jO7o/N7JFW4roW/2fwKvC31PWvbeuGpXzpIS9JDDOrwNf8xzvn5kUdT5qZ9QPeA6pa1NtFIqOWv5Q0M/uKme2eGst+Gb42nuvwUZGyo+Qvpe6LwDvAGuBk/GiiTdGGJBJ/KvuIiJQhtfxFRMpQbCez6tmzp+vXr1/UYYiIlJRFixatcc71am+/2Cb/fv36sXDhwqjDEBEpKWaW1RPlKvuIiJQhJX8RkTKk5C8iUoZiW/MXEWlPY2Mj9fX1bN68OepQiq66upqamhqqqqryOl7JX0RKVn19PV26dKFfv374SV3Lg3OOhoYG6uvr6d+/f17nUNlHRErW5s2b6dGjR1klfgAzo0ePHoF+4lHyF5GSVm6JPy3ofavsI1IozsFtt8GqVX57+HAYNmznx4gUiVr+IoXy6qvwne/AT3/qf/3oR1FHJAVw1FFHtbvPueeey+uvvw7AddftuD5PNsd37hz+SptK/iKFsmCB//2dd+CrX4WtW6ONRwri+eefb3efu+66i9raWuCzyT+b4wtByV+kUBYsgB49oH9/6NABmrSOSxKlW+VPP/00xxxzDGPGjOGAAw5g/PjxpGdNPuaYY1i4cCGTJk1i06ZNDBkyhPHjx+9w/MaNGxk+fDiHHXYYgwYN4ne/+11B41bNX6RQFiyAww8HM6ishG3boo4o0S65BF55JdxzDhkCP/959vu//PLLLF68mL333psvfelL/OlPf+LLX/7yPz+//vrrufXWW3mllUCrq6t5+OGH6dq1K2vWrOHII49k5MiRBevQVstfpBD+8Q947TUYOtRvq+VfFoYOHUpNTQ0VFRUMGTKEpUuXZn2sc47LLruMwYMHM2LECJYvX87KlSsLFqta/iKF8NJL0NzsW/6g5F8EubTQC6VTp07/fF1ZWUlTDn/n999/P6tXr2bRokVUVVXRr1+/gj65rJa/SCGkO3vTyV9lH0mpqqqisbHxM++vX7+e3r17U1VVxbx581i2LKuZmfOm5C9SCC++CPvuC336+G21/CVl4sSJDB48+J8dvmnjx49n4cKFDBo0iKlTp3LAAQcUNI7YruFbV1fntJiLlKwBA+ALX4AHH/Tb550Hc+bA8uXRxpUwb7zxBgceeGDUYUSmtfs3s0XOubr2jlXLXyRsq1fDe+9t7+wFtfwldhLX4fvJJ/D970cdRXnr2PQplc2N4BxdN62kz4a/s/s/lmM4cI5OTf+gunEDHZs2RR1qzioqjbqxA+h+zCFQU+OHcaaH4qV/f+YZ/7uSv8RY4pL/1q3+p2spjoO3vkTv5o8A2K/pTU7Y/BB1W5+ngp2XE5sxtlg1jhKalMtBJU10eiWLJ3UrK+Gww3bcVoevxEjikn+PHiqrFkV9PXzvezBr1o7vDxkCJ10Oe+zht3v0gIEDoW9fnwABdtuNit12Y5eK0qo6rloFffo47r3mfc4+5K/+jXSfWcvfBwyALl22H6yWv8RM4pK/FMjWrXDRRfDcc3572TI/jv2nP4XjjvPv9enjR7gklP/uMtZ13RdG5niflZVK/hIrSv7Svi1b4NRT4fe/h5NPhupqOPpouPRSP29NmeiQ+t+SV/WmQweVfSRWlPylde+840s6zsETT8CTT8Ltt8OFF0YdWWTSVau8k79a/ol1yy23cMcdd/DRRx9x6aWXMmnSpLzO07lzZzZu3BhydK1T8pfPWr/eLzySfsKwY0eYPBkmTow2rogFSv7pg5ubocT6OqR9t99+O0888QQ1NTVRh5I1/SuUz/r2t32H7vz5sGkTbNxY9okftpd98mrABzpY4uyCCy7g3Xff5YQTTuDmm2/moosuAmDChAlcfPHFHHXUUQwYMIBZqcERxZ66uS1q+YvX0OBr+48+CvffDz/5Cfzrv0YdVayE0vJvavI/SUn4IprT+c477+Sxxx5j3rx5/OEPf9jhsxUrVvDcc8/x5ptvMnLkSMaMGVP0qZvbouQvcPfdcO6524cpHnUUXHZZtDHFULpak3fNP++DpVSNHj2aiooKamtr/zk9c3rq5vnz51NRUfHPqZv33HPPosam5F/u1q/3o3aGDoVvftMnqdGjtycr2UHeIzZV9im8OMzp3ELmFM/pedSKPXVzW/Q/vNz97Ge+5DN37o5PpEqr8h6xGahmJElS7Kmb26LkX86WLvWtpbPOUuLPUt6zNKjlLynjx4/n5JNPZtCgQdTV1RV86ua2KPmXm9WrYdw42LABVq70k5Fde23UUZWMvJN/ZoevJE56ucYJEyYwYcIEAO65554d9kmP3+/Zsyd//vOfWz1Pscb4g4Z6lp+bb4Z586BXLxg0CO65B/bZJ+qoSkbez2qpw1diRi3/crJ+Pdx2G4wZAzNnRh1NSVLZR5JCLf9ycscdvtyT56PnEkLZRy3/0MV1NcJCC3rfSv7lYtMmX/I57jh17gagoZ7xUl1dTUNDQ9l9ATjnaGhooLq6Ou9zqOyTZM7Bd78Lb78Na9f6+ed/9KOooyppgYd6KvmHqqamhvr6elavXh11KEVXXV0daC4hJf8kmz8ffvlLqK31C4ucfz4MGxZ1VCUtcM1fZZ9QVVVV0b+MphUPk5J/kt16K3TvDi++CLvuGnU0iaAOX0mKUGr+Zna8mb1lZkvM7DO9iWb2fTN73cxeNbMnzSy5yz3FxQcfwMMP+zl7lPhDk/dQT3X4SswETv5mVgncBpwA1AKnm1lti91eBuqcc4OBWcB/B72utGPyZD93/Le+FXUkiaKWvyRFGC3/ocAS59y7zrmtwHRgVOYOzrl5zrlPU5svAKWz4kEp2bbNd+yuXAlTpvglF/v1izqqRNETvpIUYST/zwEfZGzXp95ryznAo619YGYTzWyhmS0sx977wE48EXr0gD339NM4fOc7UUeUOIGHeqrsIzFR1A5fMzsTqANaHXLinJsCTAGoq6srr4G7Qb3xBjz+OJxxBhxxBPTs6ZdilFDlPdRTZR+JmTCS/3Igc3KYmtR7OzCzEcDlwDDn3JYQriuZ7rrLJ5ibboI+faKOJrH0hK8kRRhlnxeBgWbW38w6AuOA2Zk7mNmhwGRgpHNuVQjXlExbtsC998KoUUr8BaYOX0mKwMnfOdcEXATMBd4AZjrnFpvZNWY2MrXbjUBn4EEze8XMZrdxOsnHI4/4BVnOOy/qSBIv8KyeSv4SE6HU/J1zc4A5Ld67MuP1iDCuIy2k5zOZMgX23ReOPTbaeMpAZaX/QSuvA0FlH4kNTexWqm64wa8oXlEBTz0F55yzfYVxKRiVfSQpNL1DKdq2DW65BQ491Nf5O3bUw1xFkvdQT7X8JWaU/EvRE0/Ahx/6L4BTTok6mrKioZ6SFKoTlKLf/Ab22ANOOinqSMqOyj6SFEr+pWbdOj+6Z/x46NQp6mjKjsb5S1Io+Zea6dP9cJMJE6KOpCxpqKckhWr+peD99+Gqq6Cx0S/QMmiQ7+yVolPLX5JCyb8U/OIXMHUq9O/vSz2XXw5mUUdVllTzl6RQ8o+7pia4/34YOdIvziKR0gLukhSq+cfdE0/4+fnPPjvqSIQQFnBX2UdiQsk/7u67z6/De+KJUUciqOwjyaHkH2effOJLPePGaVhnTOSd/NNTbyj5S0yo5h9Hn34KmzfDzJmwaROcdVbUEUlK3kM9zQJ8c4iET8k/bt57D2prffIH2G8/OPLIaGOSfwqUv/P+5hAJn5J/3Nx3n3+I68Yb/YRtRx+tYZ0xEij5q+UvMaLkHyfO+WGdw4bBD34QdTTSiryHeoJa/hIr6vCNk0WL4O23/bw9Ekt5D/WEgN8cIuFS8o+T3/7Wl3rGjIk6EmlD4Jq/yj4SE0r+cdHU5CdtO+kk2H33qKORNqSf1WpuzuNglX0kRlTzj9qCBbBkCbzzjn+SVyWfWEsn/6Ym/0Nazger5S8xoeQfpfXrfeduelhnr156kjfm0g/q5v2Ur1r+EhNK/lF66CGf+B9+2I/t79ULqqujjkp2ItAUPerwlRhR8o/StGkwYIBfhF1j+UtCZtknZ+rwlRhRh29UPvoInnwSzjhDib+EqOwjSaHkH5WZM/2QkTPOiDoSyUHgso9a/hITSv5ReeABOOQQOPDAqCORHARK/mr5S4yo5l9MzzwDjz8OW7fCX/4CN9wQdUSSo0A1f3X4Sowo+ReLc/D1r/vF2Csr/cgelXxKTuCav8o+EhMq+xTLCy/AsmVwzz3Q2AirVkFNTdRRSY5U9pGkUPIvlmnT/Bj+0aOjjkQCCFz2UctfYkLJvxiamvzonq9+Fbp2jToaCUBDPSUplPyL4emn/bw9p58edSQSkJ7wlaRQ8i+G6dOhSxfN25MAgWv+KvtITGi0T6G8/DLccYcf5fPgg77Wv8suUUclAQWe3kEtf4kJJf9CueoqmDsXevaEPfaACy+MOiIJQaCavzp8JUZU9imEdevgscfg4oth+XJ47z344hejjkpCoKGekhRK/oXw8MN+LP+4cVFHIiHTE76SFEr+hTB9Ouy3Hxx2WNSRSMj0hK8kRSjJ38yON7O3zGyJmU1q5fOjzewlM2sys2SvTr5qlZ+qedw4TdWcQCr7SFIETv5mVgncBpwA1AKnm1lti93eByYADwS9XuzNmuWnalbJJ5E0pbMkRRijfYYCS5xz7wKY2XRgFPB6egfn3NLUZ80hXC9+PvoIrrgCtmyBZ5+Fgw+Ggw6KOiopAA31lKQII/l/DvggY7seOCKfE5nZRGAiQN++fYNHViyTJ8Ndd0H//lBRAT/4QdQRSYFoegdJiliN83fOTQGmANTV1bmIw8mOc76D95hjYN68qKORAlPZR5IijA7f5cA+Gds1qffKw6uvwptvwmmnRR2JFIHKPpIUYST/F4GBZtbfzDoC44DZIZy3NMyY4TPCKadEHYkUgZ7wlaQInPydc03ARcBc4A1gpnNusZldY2YjAczscDOrB04FJpvZ4qDXjYV0yWf4cL8ylySehnpKUoRS83fOzQHmtHjvyozXL+LLQcmycKGfuuGKK6KORIpEyV+SIlYdviWhqQl+/GNYvRpeew2qqrQ6VxkJPL0D+OdAKvRwvURLyT9Xjz8ON9wAffr4xP+tb0H37lFHJUUSeKgn+G+Ojh1Di0kkH0r+uZoxA7p184uxd+oUdTRSZIGHeuZ9sEi49LNnLjZvhkce8WUeJf6yFHioZ94Hi4RLyT8Xc+fChg0a01/GQiv7iERMyT8XM2f6VblGjIg6EomIyj6SFEr+2dq0CWbPhq99zXf0SlkKPNQT1PKXWFCHb3vuuQdeeQU+/BA2blTJp8yFMtRTyV9iQMl/Z1avhnPP9cPyOnaEww/3E7hJ2Qql5q+yj8SAkv/OPPSQ/4/6/PMwZEjU0UgMqOwjSaGa/87MnAn77w+HHBJ1JBIT6QdzA5V91PKXGFDyb8vKlfD00zB2rNbilR3kvQ67Wv4SI0r+bUmvxasOXmkh75mZ1eErMaLk35aZM6G21q/HK5Ih7+SvDl+JEXX4ZlqwAB57zLfMnn0Wrr466ogkhiorNb2DlD4l/0zf+Aa8/rp/3bkzjB8fbTwSS3nX/NXhKzGisk/a4sU+8d9yi//PuWEDfP7zUUclMRS47KOWv8SAWv5pM2b4cXynnqqFNmSn1OErSaAsB34t3hkzYNgw2HPPqKORmMt7NUZ1+EqMKPkDvPoqvP22hnVKVlT2kSRQ8gff6q+s9DN2irQjcNlHLX+JgfKt+a9a5efuaW6G+++Hf/936NUr6qikBGiopyRB+Sb/K6+EyZO3b193XXSxSEkJPNRTyV9ioDyTf1OTn75hzBi47Tb/v3mPPaKOSkqEnvCVJCjP5P/UU9DQ4B/i6t076mikxKjDV5KgPDt8Z86ELl3g+OOjjkRKUN5DPdXhKzFSfsl/61bf0TtqFFRXRx2NlCC1/CUJyi/5P/EErFunMf2SNz3hK0lQHjX/bdtgzhzYvBnuvht23x2OOy7qqKREBR7qqbKPxEB5JP/p0+HMM7dvT5zoF2QXyYNW8pIkKI/kP20a7LMPPPqo3x44MNp4pKRVVsKWLXkeCGr5SywkP/mvWwePPw4XXwwHHRR1NJIA6vCVJEh+h+8jj0Bjozp4JTSBZ/VU8pcYSH7ynzEDBgyAurqoI5GE0MRukgTJTv5r1vihnWPHglnU0UhCqOwjSZC8mv/WrTB/vn/91FP+f+nYsdHGJImS91DP9ApxavlLDCQv+a9fD8ceu337wANhyJDo4pHEyXuop1mAbw6RcCUv+e++Ozz77PbtgQNV8pFQ5V32gQC9xSLhCiX5m9nxwC+ASuAu59z1LT7vBEwFvgA0AKc555aGce3PqKqCL3+5IKcWgYDJP9DBIuEJ3OFrZpXAbcAJQC1wupnVttjtHGCdc24/4GbghqDXFYlKoMa7Wv4SE2GM9hkKLHHOveuc2wpMB0a12GcUcG/q9SxguJlqMVKa1PKXJAgj+X8O+CBjuz71Xqv7OOeagPVAj5YnMrOJZrbQzBauXr06hNBEwqeavyRBrDp8nXNTgCkAdXV1LuJwRFoVaMCOkr+05c03Ye1a/7pzZxg8uKCXCyP5Lwf2ydiuSb3X2j71ZtYB6Ibv+BUpOXkP9QSVfaR1b7/t5x5rbvbbRxwBL7xQ0EuGkfxfBAaaWX98kh8HnNFin9nA14E/A2OAp5xzatlLSVLZR0I3bRo451cZ3G036Nat4JcMnPydc01mdhEwFz/U827n3GIzuwZY6JybDfwauM/MlgBr8V8QIiUpcIevkr9kcs6vOXL00fAf/1G0y4ZS83fOzQHmtHjvyozXm4FTw7iWSNQCD/VU2Ucy/e1vvt7/3e8W9bLJnthNpABU9pFQTZ/u/1GdckpRL6vkL5KjdPLPq9dKHb6SyTk/7fzw4dCrV1EvHauhniKlID0tf3Pz9tdZU8tfAOrrYdkyeO89ePdd+PGPix6Ckr9IjtLT8m/blkfyV4evbN0Khx0G6QdZq6th9Oiih6HkL5KjQAtyqcNX/vd/feK/8UY45BCoqYHu3YsehpK/SI4CJ3+1/MvbtGk+2V98MXTsGFkY6vAVyVGg1RjV4VvePv0Ufvc7P7InwsQPSv4iOVPLX/L2xz/Cxo1w+ulRR6LkL5KrQMlfHb7lbfp06NMHhg2LOhLV/EVylU7+eeVwdfiWnxUrYNEi/w/mj3+EiRPzGCYWPiV/kRxlDvXM62C1/MvLaaftuK74mWdGF0sGJX+RHAUu+6jlXz6WLfOJ/5JLYPx46NIF/uVfoo4KUPIXyZk6fCVr06f73y++GPr3jzaWFtThK5KjwEM9lfzLxwMPwBe/GLvED0r+IjnTE76SlcWL4dVXYzGsszVK/iI5UtlHsjJtGlRUwNixUUfSKtX8RXIUaKinOnyT7cMPfZ2/uRmmToURI/y4/hhS8hfJkYZ6SpuuuALuvnv79k03RRdLO1T2EcmRnvCVVm3eDLNmwVlnwSef+Hl8xoyJOqo2qeUvkiN1+Eqr/vAH2LABzj4bOneOOpp2qeUvkqNAQz1V9kmu+++HvfaCf/u3qCPJipK/SI70hK98xtq1ft6e00+Pxbw92VDyF8mRhnrKZ8yaBY2NfgqHEqGav0iOAg/1VPJPhvp6uOAC39G7eDEccAAcemjUUWVNLX+RHAUe6gl+HLiUtjvvhEcf9cn/85+Hn/wEzKKOKmtq+YvkKHDZB3zrP+Jl/CSA5ma47z449lh47LGoo8mLWv4iOQrc4Zv3wRIbzz4L77/vh3WWKCV/kRwFHuqZ98ESG1On+rH8o0dHHUnelPxFchRa2UdK06efwoMPwqmnwq67Rh1N3lTzF8mRyj5l6v33/Xj+efP89A0lXPIBJX+RnAVewD3vgyUyy5bBwIF+LD/AvvvC0UdHG1NASv4iOQo01FMt/9J0770+8d93n6/1H3ywn6u/hCn5i+RINf8y09wMv/kNDB8OZ54ZdTShKe2vLpEIKPmXmWeegaVL4RvfiDqSUCn5i+Qo8PQOoLJPKfnNb6BbN/ja16KOJFRK/iI5CmV6h2eegccf96NGJL7Wr/eTto0bB7vsEnU0oVLyF8lRoMZ79+7+94kT4Stf8fPBSPycfjr07Qv77w+bNsE3vxl1RKEL1OFrZnsAM4B+wFJgrHNuXSv7PQYcCTznnDspyDVFohao7DNiBLz8sk8oZ5zhx45LvLz8sl+Effhw/wXQrx8cfnjUUYUu6GifScCTzrnrzWxSavvSVva7EdgVOD/g9UQiF6jsYwZDhvjXffvCypWhxSUhmTwZqqv9U7zpn9QSKGjZZxRwb+r1vUCrE104554EVNyURAitz7Z3b1i1KnA8EqKNG/1yjKedlujED8GTfx/n3IrU64+APgHPJxJ7oSX/Pn2U/ONm2jT/BXB+8osU7ZZ9zOwJYM9WPro8c8M558zMBQnGzCYCEwH69u0b5FQiBZN+sDPwUP3evf1cMY2NUFUVOC4JweTJMGgQHHlk1JEUXLvJ3zk3oq3PzGylme3lnFthZnsBgZoxzrkpwBSAurq6QF8kIoXUoUNILX+A1ath770DxyR5uuoquP12cA4aGuCXvyypFbnyFbTDdzbwdeD61O+/CxyRSAmorAyp5g++01fJPxoNDXDjjb61X1cHu+2WuCd52xI0+V8PzDSzc4BlwFgAM6sDLnDOnZvafhY4AOhsZvXAOc65uQGvLRKZUNZhTyd/1f2j86tf+WG3d93lvwDKSKDk75xrAIa38v5C4NyM7X8Nch2RuAm17KPkH43GRrj1Vj+ev8wSP2hWT5G8hF72keJ76CFYvhzuvDPqSCKh5C+Sh1CSf5cu/mEitfyLZ9UqmDnT/+X9+tew335w4olRRxUJJX+RPIRS8zfzrX+1/Ivn+9/3D3GlTZlS8ouy5Ks871okoFBq/qCnfIvprbf8Q1zf+55/vmL9ejjvvKijioxa/iJ5CKXsA77Td8WK9veT4K691pfZJk1K/NQN2VDLXyQPoZR9QGWfYnn7bXjgAbjwwu0d7WVOLX+RPIRW9knP7+NcWTxVWlQffwz/9V+wbh288gp06gQ//GHUUcWGkr9IHkIr+/Tu7cebf/yxShFhu/RS//DWXnv57Suv3P5shSj5i+Qj1OQPvvWv5B+eBQv807uXXAI33RR1NLGkmr9IHkKr+esp3/Bt2+Zr+3vuCVdfHXU0saWWv0geQh3qCer0Daqx0Sf6ZctgzRp46SU/rLNr16gjiy0lf5E8hDrUE9TyD+qSS/y0zAMG+I7z88/3q3FJm5T8RfIQWtmnRw+frNTyz9+vfuUT/w9+4Kdnlqwo+YvkIbSyT4cGQuYSAAAFiklEQVQO/gtALf/sNTbCddf54ZvOwZw5cNxxcP31UUdWUpT8RfIQWtkHtJZvLtasgbFjYd48qK31X55f+QpMnbp9cWXJipK/SB4qK2HLlpBOpqd827Z8uV9WccECv/3mm35ennvvhbPPjja2EqfkL5KH0Gr+4Fv+ixaFdLIS4Rxs3OiXUWxu9ttLlsBf/uITfPrzxx7znw8d6he5HzwYrrnGb0sgSv4ieejQAf76VzjooODnmvRRb8avXUJD1V5t7mPO7eQMO/sMbCef5/tZu8fuNF7o5DZR7TZ/5v1mjOVV/WiyKpqp4M9dL2TqHt9j+Yb+focGoAyW2B082I9ULSQlf5E8nH++X+s7DK987lz6vLNlpwnTtTvvT9ufu5181v652zl2Z9fdyXmbKjqxvlNvPunYg20VPg2t3aWGd7rXsalqx7H53VK/ykn//oW/hrl2vqGjUldX5xYuXBh1GCIiJcXMFjnn6trbT9M7iIiUISV/EZEypOQvIlKGlPxFRMqQkr+ISBlS8hcRKUNK/iIiZUjJX0SkDMX2IS8zWw0sy+PQnsCakMMpBeV43+V4z1Ce9617zt6+zrle7e0U2+SfLzNbmM3TbUlTjvddjvcM5XnfuufwqewjIlKGlPxFRMpQEpP/lKgDiEg53nc53jOU533rnkOWuJq/iIi0L4ktfxERaYeSv4hIGSrZ5G9mx5vZW2a2xMwmtfJ5JzObkfr8L2bWr/hRhi+L+/6+mb1uZq+a2ZNmtm8UcYapvXvO2O8UM3NmVvJDArO5ZzMbm/q7XmxmDxQ7xkLI4t93XzObZ2Yvp/6NnxhFnGExs7vNbJWZvdbG52Zmt6T+PF41s8NCu7hzruR+AZXAO8AAoCPwV6C2xT7fAu5MvR4HzIg67iLd978Bu6ZeX1jq953NPaf26wLMB14A6qKOuwh/zwOBl4Huqe3eUcddpPueAlyYel0LLI067oD3fDRwGPBaG5+fCDyKX0/zSOAvYV27VFv+Q4Elzrl3nXNbgenAqBb7jALuTb2eBQw3a3ch1Lhr976dc/Occ5+mNl8AaoocY9iy+bsG+ClwA/DZVcFLTzb3fB5wm3NuHYBzblWRYyyEbO7bAelFfrsBHxYxvtA55+YDa3eyyyhgqvNeAHY3s73CuHapJv/PAR9kbNen3mt1H+dcE7Ae6FGU6Aonm/vOdA6+1VDK2r3n1I/C+zjn/ljMwAoom7/n/YH9zexPZvaCmR1ftOgKJ5v7vho408zqgTnAd4oTWmRy/T+ftQ5hnETix8zOBOqAYVHHUkhmVgHcBEyIOJRi64Av/RyD/+luvpkNcs59HGlUhXc6cI9z7n/M7IvAfWZ2sHOuOerASk2ptvyXA/tkbNek3mt1HzPrgP8RsaEo0RVONveNmY0ALgdGOue2FCm2QmnvnrsABwNPm9lSfF10dol3+mbz91wPzHbONTrn3gPexn8ZlLJs7vscYCaAc+7PQDV+ArSkyur/fD5KNfm/CAw0s/5m1hHfoTu7xT6zga+nXo8BnnKpHpQS1u59m9mhwGR84k9CHXin9+ycW++c6+mc6+ec64fv5xjpnFsYTbihyObf9yP4Vj9m1hNfBnq3mEEWQDb3/T4wHMDMDsQn/9VFjbK4ZgNnp0b9HAmsd86tCOPEJVn2cc41mdlFwFz8CIG7nXOLzewaYKFzbjbwa/yPhEvwHSrjoos4HFne941AZ+DBVP/2+865kZEFHVCW95woWd7zXOA4M3sd2Ab80DlX0j/ZZnnf/wn8ysy+h+/8nVDKjTozm4b/Eu+Z6se4CqgCcM7die/XOBFYAnwKfCO0a5fwn5uIiOSpVMs+IiISgJK/iEgZUvIXESlDSv4iImVIyV9EpAwp+YuIlCElfxGRMvR/6saJ1TyHsdYAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "def equilibrium(m0):\n", " return .5*m0**2\n", "\n", "def initialize(mesh, la):\n", " ug, ud = 0.25, -0.15\n", " xmin, xmax = .5*np.sum(mesh[:2]), .5*np.sum(mesh[-2:])\n", " xc = xmin + .5*(xmax-xmin)\n", " m0 = ug*(meshxc) + .5*(ug+ud)*(mesh==xc)\n", " m1 = equilibrium(m0)\n", " f0 = np.empty(m0.shape)\n", " f1 = np.empty(m0.shape)\n", " return f0, f1, m0, m1\n", "\n", "def relaxation(m0, m1, s):\n", " m1[:] = (1-s)*m1 + s*equilibrium(m0)\n", " \n", "# parameters\n", "Tf = 1. # final time\n", "N = 128 # number of points in space\n", "la = 1. # scheme velocity\n", "s = 1.8 # relaxation parameter\n", "# initialization\n", "x = mesh(N) # mesh\n", "dx = x[1]-x[0] # space step\n", "dt = dx/la # time step\n", "f0, f1, m0, m1 = initialize(x, la)\n", "plt.figure(1)\n", "plt.plot(x[1:-1], m0[1:-1], 'b', label='initial')\n", "# time loops\n", "t = 0.\n", "while (t" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "def transport(f0, f1):\n", " # Neumann boundary conditions\n", " f0[-1] = f0[-2]\n", " f1[0] = f1[1]\n", " # transport\n", " f0[1:-1] = f0[2:]\n", " f1[1:-1] = f1[:-2]\n", " \n", "# parameters\n", "Tf = 1. # final time\n", "N = 128 # number of points in space\n", "la = 1. # scheme velocity\n", "s = 1.8 # relaxation parameter\n", "\n", "# initialization\n", "x = mesh(N) # mesh\n", "dx = x[1]-x[0] # space step\n", "dt = dx/la # time step\n", "f0, f1, m0, m1 = initialize(x, la)\n", "plt.figure(1)\n", "plt.plot(x[1:-1], m0[1:-1], 'b', label='initial')\n", "# time loops\n", "t = 0.\n", "while (t