{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Finite Volume Burgers\n", "\n", "Copyright (C) 2010-2020 Luke Olson
\n", "Copyright (C) 2020 Andreas Kloeckner\n", "\n", "
\n", "MIT License\n", "Permission is hereby granted, free of charge, to any person obtaining a copy\n", "of this software and associated documentation files (the \"Software\"), to deal\n", "in the Software without restriction, including without limitation the rights\n", "to use, copy, modify, merge, publish, distribute, sublicense, and/or sell\n", "copies of the Software, and to permit persons to whom the Software is\n", "furnished to do so, subject to the following conditions:\n", "\n", "The above copyright notice and this permission notice shall be included in\n", "all copies or substantial portions of the Software.\n", "\n", "THE SOFTWARE IS PROVIDED \"AS IS\", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR\n", "IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,\n", "FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE\n", "AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER\n", "LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,\n", "OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN\n", "THE SOFTWARE.\n", "
" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "import scipy as sp\n", "import matplotlib.pyplot as plt\n", "\n", "from matplotlib import animation\n", "from IPython.display import HTML" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "def gaussian(x):\n", " u = sp.exp(-100 * (x - 0.25)**2)\n", " return u\n", "\n", "\n", "def step(x):\n", " u = np.zeros(x.shape)\n", " for j in range(len(x)):\n", " if (x[j] >= 0.6) and (x[j] <= 0.8):\n", " u[j] = 1.0\n", "\n", " return u\n", "\n", "def g1(x):\n", " return 1+gaussian(x)\n", "\n", "def g2(x):\n", " return 1+gaussian(x) + step(x)\n", "\n", "g = g1" ] }, { "cell_type": "code", "execution_count": 20, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "tsteps = 250\n", " dx = 0.00609756\n", "lambda = 0.95\n" ] }, { "data": { "text/plain": [ "[]" ] }, "execution_count": 20, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD4CAYAAAD8Zh1EAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3deXRdZ33u8e/v6GiwrHmyNVm240nybMs4hAAhCZCkTE2BNrchbVYgDaGsTveWdrWFe0tX7+2l0NuZpjSktBCgTS71pUkhEEJwZtmx5UF2YseTLMmSNVuzdN77xznHuIplyfI5Z+999HzW0rKsvbX3b+vYj/Z533e/rznnEBGR4At5XYCIiCSGAl1EJE0o0EVE0oQCXUQkTSjQRUTSRNirE5eVlbnly5d7dXoRkUDas2fPeedc+eW2eRboy5cvp6mpyavTi4gEkpmdmmmbmlxERNKEAl1EJE0o0EVE0oQCXUQkTSjQRUTSxKyBbma1ZvYjMztsZofM7Ncus4+Z2V+Y2TEzazazbckpV0REZjKXYYuTwG855/aaWT6wx8yecs4dvmSf24HVsY+dwN/G/hQRkRSZ9Q7dOdfunNsb+3wQaAGqp+32QeBrLupFoMjMKhNerVy14fFJ/u+rreza38boxJTX5YhIEl3Vg0VmthzYCrw0bVM1cOaSv7fGvtY+7fvvB+4HWLZs2dVVKlftH58/yRe+d5QLY5MAFOSE+f2faeCjO2o9rkxEkmHOnaJmlgc8Bvy6c25gPidzzj3knGt0zjWWl1/2yVVJkO82t/G5XYfYuqyIb//KW/nGJ3ZSX1nA7zzezNNHznldnogkwZwC3cwyiYb5151zj19ml7PApbd9NbGviQf2n+njt769n8a6Yv7+nkbesqKEG64r46v37qChqoBPf+NVjnYMel2miCTYXEa5GPAPQItz7ksz7LYLuCc22uV6oN851z7DvpJEzjk+u+sQxblZPHRPIzmZGRe35WaF+Ydf2kF2Zgaf/+7hKxxFRIJoLnfobwM+BtxsZvtiH3eY2QNm9kBsnyeAN4BjwN8DDyanXJnND1o62X+mj99492pKFme9afuSghwevOk6dh87z/PHz3tQoYgky6ydos653YDNso8DPpWoomR+IhHHF79/lBVli/m5bTUz7nf39XV85Scn+NPvHeWxT5YSfRMmIkGnJ0XTyPcOdXCkY5Bfv3U14YyZX9qczAw+fcsq9p7u4yev6y5dJF0o0NPIo6+coaowh/dtqpp1349sr6VkcRbfeuXMrPuKSDAo0NNEe/8IP3m9iw9vryEjNHsTSlY4xAe3VPHU4XP0Do2noEIRSTYFepp4fO9ZnIMPb5/7Q0Mf2V7L+FSEf9unEaYi6UCBngacc/xL0xl2rihhWWnunL+voaqADdUF/Mue1iRWJyKpokBPA3tP93Kye5iPNF79I/0f2V7LobYBjnTM6+FfEfERBXoa+P7hc4RDxnvWL7nq771941IAnjqk6QBEgk6BngZ+cPgc168spSAn86q/tyI/hy21RfygRYEuEnQK9IB7o+sCx7uGuLW+Yt7HeHfDEva39nNuYDSBlYlIqinQA+6HLZ0A3FJ/9c0tcbfGvjd+LBEJJgV6wD3Vco51S/OpLZn76Jbp1izJo7ZkkZpdRAJOgR5gfcPjNJ3suXiHPV9mxq31S9h97Dwj41rVSCSoFOgB9uIb3UQc3LT22hcLeeeacsYnIzSd6klAZSLiBQV6gD1/vJvcrAw21RRd87F2LC8hHDKeO9adgMpExAsK9AB77th53rKihKzwtb+Mi7PDbKkt4gXNkS4SWAr0gDo3MMrxriFuuK40Yce8YVUZB8720z8ykbBjikjqKNAD6oXj0aaRG64rS9gxb7iulIiDl0+oHV0kiBToAfXcsfMULsqkobIgYcfcuqyI7HCI546p2UUkiBToAfX88W7eurKU0BzmPp+r7HAGO5aXaK1RkYBSoAdQW98IZ/tG2LmyJOHH3rmihNfOXaBvWIteiASNAj2Amk71AtBYl/hA3768GIhOySsiwaJAD6C9p3pZlJlBfWV+wo+9pbaIjJCx55QCXSRoFOgB1HSqhy21RYQzEv/y5WaFaagsoOmkAl0kaBToATM0NklL+yCNsaaRZNheV8z+1j4mpiJJO4eIJJ4CPWD2n+ljKuLYXpe8QG9cXszoRITDbVqWTiRIFOgBs+dUL2awdVly79Dj5xKR4FCgB0zTqV7WVORTuOjql5ubq8rCRVQXLVKgiwSMAj1AnHPsO9PH1mXXPrvibLYuK+JVDV0UCRQFeoCc7B6mf2SCLbXJD/QttUW09Y/SOah1RkWCQoEeIM2tfQAJmf98NptjvzSaz/Qn/VwikhgK9ADZd6aPnMwQa5bkJf1c66sKCNlPf4mIiP/NGuhm9rCZdZrZwRm2F5rZ/zOz/WZ2yMzuTXyZAtDc2s+GqsKkPFA0XW5WmDVL8tnXqjt0kaCYSzI8Atx2he2fAg475zYDNwFfNLOsay9NLjUxFeHg2f6LTSGpsLmmiObWPpxzKTuniMzfrIHunHsWuNKKBw7INzMD8mL7TiamPIl77dwgY5MRNtUUpuycm2uL6Bue4HTPcMrOKSLzl4j37n8F1ANtwAHg15xzl31m3MzuN7MmM2vq6upKwKkXjv2xzslUjHCJ21wb/eWxX80uIoGQiEB/L7APqAK2AH9lZpddRsc595BzrtE511heXp6AUy8cza19FOVmsqwkN2XnXLMkn+xwiP1n1DEqEgSJCPR7gcdd1DHgBLAuAceVSzS39rOxupBoy1ZqZGaEaKgq4MBZ3aGLBEEiAv00cAuAmS0B1gJvJOC4EjM2OcVr5wbZUJ269vO49VUFtLQNEImoY1TE7+YybPFR4AVgrZm1mtl9ZvaAmT0Q2+XzwA1mdgD4IfAZ55wWpUyg1zouMBlxbKjyItALGRyb5EyvOkZF/C482w7Oubtm2d4GvCdhFcmbHGyLNnlsqL5s10RSra+KnvNQ2wB1pYtTfn4RmTs9KRoAB8/2k58TTmmHaNyaJflkhIxDbWpHF/E7BXoAHGwbYENVajtE43IyM1hdkcchLXYh4nsKdJ+bmIrQ0j7gSXNLXENVgQJdJAAU6D53vOsC45MRT0a4xK2vKqRrcExT6Yr4nALd5w6ejd4Zr/dghEvcpR2jIuJfCnSfO3i2n9ysDFaUeTfCpCEW6Fo0WsTfFOg+19I+wLql0ZEmXinIiU45oJEuIv6mQPcx51w00Cu96xCNW6+OURHfU6D7WFv/KAOjk9T7JNBPdQ8zMDrhdSkiMgMFuo+1xO6IGyrzPa7kp52yLbpLF/EtBbqPtbRHw3PtUu/v0Bs00kXE9xToPtbSMUBdaS552bNOuZN0FfnZlOVlKdBFfEyB7mNH2gep98HdOYCZ0VBVqJEuIj6mQPep4fFJTnQPsc4H7edx66sKONZ5gbHJKa9LEZHLUKD71NGOQZzDFyNc4tZXFTAZcbzWccHrUkTkMhToPtXSPghAg68CPTrSRc0uIv6kQPeplvYB8rPD1BQv8rqUi+pKclmclXFx9I2I+IsC3aeiT4jmezIH+kxCIWPt0vyL7x5ExF8U6D7knONIx6Cv2s/j6isLaOkYwDktGi3iNwp0H2rtHeHCmD8e+Z+uoaqAwdFJWntHvC5FRKZRoPvQ4Vgb9bql/hmyGBf/JXNY7egivqNA96GW9gHMYK0PA33d0nzMUMeoiA8p0H2opX2AFaWLyc3y/pH/6XKzwiwvXaxAF/EhBboPtbT7s0M0rr5SI11E/EiB7jODoxOc7hmm3keP/E9Xv7SA0z3DDGpudBFfUaD7zGvnone+fr5Dj0+le7RDd+kifqJA95nDsaYMPyw7NxONdBHxJwW6z7S0D1CQE6aqMMfrUmZUWZhD4aJMdYyK+IwC3Wda2georyzw1SP/05kZ9ZX5F99NiIg/KNB9JBJxHPXpI//T1VcWcLRjgKmIpgAQ8YtZA93MHjazTjM7eIV9bjKzfWZ2yMx+nNgSF45TPcMMj0/5asrcmTRUFjA6EeFk95DXpYhIzFzu0B8Bbptpo5kVAX8DfMA5tx74SGJKW3iOxNqkg3KHDnBYa4yK+Masge6cexboucIu/wV43Dl3OrZ/Z4JqW3Ba2gcIGaxekud1KbNavSSPcMjUMSriI4loQ18DFJvZM2a2x8zumWlHM7vfzJrMrKmrqysBp04vh9sHWVmeR05mhtelzCo7nMF15XkKdBEfSUSgh4HtwM8A7wX+wMzWXG5H59xDzrlG51xjeXl5Ak6dXuIjXIJCUwCI+EsiAr0V+J5zbsg5dx54FticgOMuKP0jE5ztG/H1I//TNVQV0DEwSu/QuNeliAiJCfR/A240s7CZ5QI7gZYEHHdBCVKHaFy8VjW7iPjDrPOzmtmjwE1AmZm1Ap8DMgGcc192zrWY2X8AzUAE+IpzbsYhjnJ5R2LzogRhyGLcpVMA3LCqzONqRGTWQHfO3TWHfb4AfCEhFS1QLe0DlCzOoiI/2+tS5qwsL5vy/GzN6SLiE3pS1Cda2gdiqwH595H/y6mvLFDHqIhPKNB9YCriOHouGI/8T1dfmc+xzkHGJyNelyKy4CnQfeDE+SFGJyKBDPSGygImphzHuy54XYrIgqdA94GWiyNcgjNkMa5BUwCI+IYC3QeOdAwQDhmrKvz/yP90K8oWkxUOaeiiiA8o0H2gpX2QVRV5ZIf9/8j/dOGMEGuX5NPSoUAX8ZoC3QfiI1yCKj4FgHOaG13ESwp0j/UOjdPePxrIDtG4hsoCeobG6Rwc87oUkQVNge6xeNvz+qpCjyuZPy0aLeIPCnSPHQ7wCJe4dRrpIuILCnSPHWobYGlBDqV5wXnkf7rCRZlUFy3SSBcRjynQPXa4bYCGquC2n8dFpwBQoIt4SYHuodGJKY51XWB9GgR6Q1VB7InXKa9LEVmwFOgeev3cBaYiLlBT5s6koTKfiIOjHZqoS8QrCnQPHWrrB0ibJhfQSBcRLynQPXS4fYC87DC1xblel3LNaotzWZyVoXZ0EQ8p0D10uG2A+sp8QqFgzYF+OaGQsU4doyKeUqB7JBJxtLQPBPqBoukaKgs4oikARDyjQPfIqZ5hhsan0qJDNK6+soDBsUlae0e8LkVkQVKgeyT+VGU6dIjGxZ92VceoiDcU6B453N4f2DnQZ7J2aT5mmgJAxCsKdI8cbhtgVUUeOZnBmwN9JrlZYVaULlbHqIhHFOgeOZQmj/xPV19ZoMUuRDyiQPdA1+AYnYNjadUhGtdQVcCZnhEGRye8LkVkwVGgeyDeJJGed+jRjtEjmgJAJOUU6B6IjwJJxzv0es2NLuIZBboHDrUNUF20iKLcLK9LSbilBTkU5WaqY1TEAwp0Dxw625+WzS0AZsb6qgIO6Q5dJOUU6CnWPzLBG+eH2FyTPo/8T7exuogjHQOMTWpudJFUUqCn2KGz0SlzN9YUeVxJ8myqKWRiynGkXR2jIqmkQE+x/a3RQN9Unb536Jti7z6aW/s8rkRkYZk10M3sYTPrNLODs+y3w8wmzezDiSsv/TS39rGsJJfixenXIRpXXbSI0sVZF395iUhqzOUO/RHgtivtYGYZwJ8A309ATWmtubX/4h1sujIzNtYUckCBLpJSswa6c+5ZoGeW3T4NPAZ0JqKodNV9YYyzfSNpH+gAm2qKeL1zkOHxSa9LEVkwrrkN3cyqgZ8F/nYO+95vZk1m1tTV1XWtpw6c5nj7eRp3iMZtrikk4uDgWQ1fFEmVRHSK/h/gM865yGw7Oucecs41Oucay8vLE3DqYGlu7ccMNqRxh2jcRnWMiqRcOAHHaAS+aWYAZcAdZjbpnPtOAo6dVppb+7iuPI+87ET82P2tIj+HysKci+9KRCT5rjlZnHMr4p+b2SPAdxXmb+acY39rP+9YU+Z1KSmzqaZQd+giKTSXYYuPAi8Aa82s1czuM7MHzOyB5JeXPtr7Rzl/YYzNC6D9PG5TTREnu4fpH9ZUuiKpMOsdunPurrkezDn3y9dUTRr7aYdo+refx8Wv9cDZfm5cvXDemYh4RU+Kpkhzax/hkF2cXnYh2FQdfTeyX80uIimhQE+R5tZ+1i7NT6s1RGdTmJtJXWmu2tFFUkSBngLOOZpb+xbE+PPpNtUU6YlRkRRRoKfAqe5hBkYn03rK3JlsrimkrX+UrsExr0sRSXsK9BSItyFvXICBvrFaDxiJpIoCPQX2n+knOxxizZJ8r0tJuQ3VhYQM9p1RoIskmwI9Bfac7mVzbRGZGQvvx704O8y6pQXsPd3rdSkiaW/hJUyKjU5McehsP9vrir0uxTPb64rZd7qPyalZp/sRkWugQE+y5tZ+JiOO7csWdqAPjU9x9JyWpBNJJgV6ku05FW1q2LbA79AB9p5Ss4tIMinQk2zPqV5Wli2mJI2XnJtNTfEiyvOzL/5yE5HkUKAnkXOOvad7F/TdOUSXpNu+rJg96hgVSSoFehKd7B6mZ2h8QXeIxm2vK+ZMzwidA6NelyKSthToSdR0MroUqwL9p30IanYRSR4FehK9fKKHotxMVpXneV2K5zZWF5KTGeKlE7OtNy4i86VAT6KXTvTwluUlhELmdSmeywqH2LasWIEukkQK9CRp7x/hdM8wO1eWel2Kb+xcUcqRjgGtYCSSJAr0JHnpjeid6M4VJR5X4h87V5bgHLxyUnfpIsmgQE+Sl070kJ8TXlArFM1mS20RWRkhXjrR7XUpImlJgZ4kL53oZsfyEjLUfn5RTmYGW2qLeFnt6CJJoUBPgs7BUd7oGlJzy2XsXFnCwbYBLoxNel2KSNpRoCfBi/H2c3WIvsn1K0uZijheVrOLSMIp0JNg9+tdFOSEL67WIz+1va6Y7HCI3a8r0EUSTYGeYM45dr9+nhuuK1P7+WXkZGawY3kJzx0773UpImlHgZ5gJ7uHaesf5cbVZV6X4ls3ri7j6LlBzesikmAK9ATb/XoXADeuUqDPJP6z2a27dJGEUqAn2O5j56kpXkRdaa7XpfhWQ2UBxbmZCnSRBFOgJ9DkVITnj3dz46oyzNR+PpNQyLhhVRm7Xz+Pc87rckTShgI9gfad6WNwdFLt53PwjtVldA6OcaRD64yKJIoCPYGePtJJOGS8fXW516X43rvWVgDRn5mIJMasgW5mD5tZp5kdnGH7L5pZs5kdMLPnzWxz4ssMhqePdLJjeQmFizK9LsX3Kgpy2FhdqEAXSaC53KE/Atx2he0ngHc65zYCnwceSkBdgdPaO8yRjkFuXlfhdSmBcfO6Cvae7qVnaNzrUkTSwqyB7px7FphxNiXn3PPOufi6Yi8CNQmqLVB+FLvTvLlegT5Xt9RX4Bw8c1R36SKJkOg29PuAJ2faaGb3m1mTmTV1dXUl+NTeevpIJ3WluawsW+x1KYGxoaqQsrxsNbuIJEjCAt3M3kU00D8z0z7OuYecc43Oucby8vTpOBwam+T5493cvK5CwxWvQihk3LyunB8f7WJscsrrckQCLyGBbmabgK8AH3TOLbhZl54+0snYZITb1i/1upTAuW3DUgbHJnn+2IL7ZyOScNcc6Ga2DHgc+Jhz7rVrLyl4njzYTlleNo3LNf/51XrbqjLys8M8caDd61JEAm8uwxYfBV4A1ppZq5ndZ2YPmNkDsV0+C5QCf2Nm+8ysKYn1+s7w+CRPH+nk9g1LNbviPGSHM3h3wxK+f/gcE1MRr8sRCbTwbDs45+6aZfvHgY8nrKKAeeZoF6MTEW7fqOaW+bp9YyWPv3qW549388416dO3IpJqelL0Gj1xoJ3SxVm8Rc0t8/b21WXkZYd5olnNLiLXQoF+DQZHJ/hByzlu27CUcIZ+lPOVk5nBrfUVPHmwndEJjXYRmS+l0DV48mAHoxMR7ty2IJ+lSqg7t9UwMDqpMeki10CBfg0e29PKirLFbFtW5HUpgfe2VWUsKcjmsT2tXpciElgK9Hk60zPMSyd6uHNrtR4mSoCMkPGhrdU881oXXYNjXpcjEkgK9Hl6fO9ZAH52W7XHlaSPD2+rYSri+Ld9Z70uRSSQFOjzMBVxfLvpDDdcV0pNsZaaS5TVS/LZXFvEN185o5WMROZBgT4PP2w5x9m+Ee55a53XpaSdj11fx7HOC7xwXFMBiFwtBfo8fO2FU1QW5nBr/RKvS0k779tUSXFuJv/4wkmvSxEJHAX6VTrWOcjuY+e5+/o6jT1PgpzMDH5+xzKeOhx9FyQic6dEukqPPH+SrIwQP7+j1utS0tYv7lwGwD+9cMrjSkSCRYF+FToHRvl2Uyt3bqumLC/b63LSVm1JLndsrOSfXzxF//CE1+WIBIYC/Sr8/U/eYHIqwidvus7rUtLep961igtjkzzy/EmvSxEJDAX6HPUOjfP1l07zgc1V1JVqmblkq68s4Nb6JTz83AkujE16XY5IICjQ5+jvnn2D4fEpHnzXKq9LWTB+9eZV9I9M8NXdJ7wuRSQQFOhzcKZnmIefO8GdW6tZsyTf63IWjC21RbynYQlf/vFxOgdHvS5HxPcU6HPwhe8dJWTwX9+71utSFpzfvaOesckIf/bU616XIuJ7CvRZ7DnVy679bXzi7SupKlrkdTkLzoqyxXzsrXV865XTHG4b8LocEV9ToF/B2OQUn3msmarCHH7lnRrZ4pVfu2U1JYuz+MxjzUxq3VGRGSnQr+Cvnz7Gsc4L/PGdG8nLnnX5VUmSotws/scHNnDgbD//oA5SkRkp0Gew/0wff/PMce7cVs1Nayu8LmfBu2PjUt67fglfeuo1jnYMel2OiC8p0C+jd2icB7++lyUFOXz2fQ1elyOAmfFHH9pIwaJMPvnPexgc1ROkItMp0KeZijh+49v76Boc42/v3kZRbpbXJUlMeX42f3XXVk71DPPb/9pMJKI500UupUC/hHOO3//OQZ452sVn39/AphqtFeo3O1eW8pnb1vLkwQ7++IkWLYQhcgn19MU45/jT7x/l0ZdP8+BN13H39Vq8wq8+8faVnO0d4Su7T1CSl8WDN+npXRFQoAMQiTg+/++H+epzJ/mFHbX8Nz1A5Gtmxufev57e4Qn+938cZWR8it989xot1i0L3oIP9KGxSX77sWb+vbmde9+2nD/4mQYFQwCEQsaXPrqZ3KwM/vLpY3T0j/L5D20gJzPD69JEPLOgA72lfYBf/cZeTpwf4nduX8evvGOlwjxAwhkh/uedG6koyOEvfvg6B9sG+Mu7trKqIs/r0kQ8sSA7RYfHJ/lfTx7h/X+5m/6RSf754zt54J3XKcwDyMz4zXev4av37qCjf4Q7/vwn/NlTrzE6MeV1aSIpZ16NEmhsbHRNTU0pPeeFsUm++fJpvvzj45y/MM5HG2v43dvrKV6soYnpoHNwlD/6bgu79rexpCCbB29axUcba1mUpWYYSR9mtsc513jZbeke6FMRR9PJHnbtb+M7r55laHyKG1eV8RvvXsP2uuKkn19S78U3uvnS91/j5ZM95OeE+bltNbx/cyVba4sJhfQuTILtmgLdzB4G3gd0Ouc2XGa7AX8O3AEMA7/snNs7W1HJCvTxyQiH2vrZc6qXppO9vHKyh+6hcbLDId63qYq7r1/G1mUK8nTnnOOVk718/aVTPHmgg/GpCOX52exYXsz2uhIa64ppqCogM2NBtjpKgF1roL8DuAB8bYZAvwP4NNFA3wn8uXNu52xFzTfQ+4cneL1zkO6hcXpiHx39o5zsHuJU9zBn+0aYij1BuKwkl8a6Ym6ur+BdaytYrAm2FqSB0QmebunkR0c7aTrZy9m+EQAyQkZN8SKWly6mrjSXJQU5lC7OomRxFqV5WeRmhcnNymBRVgaLMjPIzQqToTt88diVAn3WhHPOPWtmy6+wyweJhr0DXjSzIjOrdM61z6vaWTz7eheffvTV//S1wkWZLC/NZUttER/aUsW6ygIa64qpKMhJRgkSMAU5mXxoazUf2loNQHv/CE0neznaMXjxRmDv6V4GR2dfuzQcMjLiH2ZkZET/DIWMcMgIWXTb9P71y/0amN4J/6Z9EnEM8aWf31HLx9++MuHHTcQtazVw5pK/t8a+9qZAN7P7gfsBli1bNq+T7VxZwiP37qAsL5uS2N2Uxh7L1agsXMT7Ny/i/Zv/89dHJ6boHhqn+8IYPUPjDI9PMTI+xfDEFCPjk4yMRxibnGLKOaamHFPOEYk4JiOOiHNMxT+fNsfM5d4DT39jPH2f6e+cL/s++k3H0DQIQVGWl52U46a0DcI59xDwEESbXOZzjIr8HCrW6s5bEi8nM4PqokVUa2UqCahE9AidBWov+XtN7GsiIpJCiQj0XcA9FnU90J+s9nMREZnZrE0uZvYocBNQZmatwOeATADn3JeBJ4iOcDlGdNjivckqVkREZjaXUS53zbLdAZ9KWEUiIjIveqpCRCRNKNBFRNKEAl1EJE0o0EVE0oRnsy2aWRdwap7fXgacT2A5QbEQr1vXvDAsxGuG+V13nXOu/HIbPAv0a2FmTTNNTpPOFuJ165oXhoV4zZD461aTi4hImlCgi4ikiaAG+kNeF+CRhXjduuaFYSFeMyT4ugPZhi4iIm8W1Dt0ERGZRoEuIpImfB3oZnabmR01s2Nm9juX2Z5tZt+KbX9plqXyAmEO1/ybZnbYzJrN7IdmVudFnYk223Vfst/PmZkzs8APcZvLNZvZR2Ov9yEz+0aqa0y0Ofz7XmZmPzKzV2P/xu/wos5EMrOHzazTzA7OsN3M7C9iP5NmM9s275M553z5AWQAx4GVQBawH2iYts+DwJdjn/8C8C2v607BNb8LyI19/smgX/Ncrzu2Xz7wLPAi0Oh13Sl4rVcDrwLFsb9XeF13Cq75IeCTsc8bgJNe152A634HsA04OMP2O4AniS4Jez3w0nzP5ec79LcAx5xzbzjnxoFvEl2Q+lIfBP4x9vm/ArfY9FVzg2XWa3bO/cg5Nxz764tEV4gKurm81gCfB/4EGE1lcUkyl2v+BPDXzrleAOdcZ4prTLS5XLMDCmKfFwJtKawvKZxzzwI9V9jlg8DXXNSLQJGZVc7nXH4O9JkWn77sPs65SaAfKE1Jdckxl2u+1Fv9J6QAAAIESURBVH1Ef7MH3azXHXsbWuuc+/dUFpZEc3mt1wBrzOw5M3vRzG5LWXXJMZdr/u/A3bHFdJ4APp2a0jx1tf/vZ5TSRaIlcczsbqAReKfXtSSbmYWALwG/7HEpqRYm2uxyE9F3Ys+a2UbnXJ+nVSXXXcAjzrkvmtlbgX8ysw3OuYjXhQWBn+/Q57L49MV9zCxM9C1ad0qqS445LbhtZrcCvwd8wDk3lqLakmm2684HNgDPmNlJou2MuwLeMTqX17oV2OWcm3DOnQBeIxrwQTWXa74P+DaAc+4FIIfoBFbpbE7/7+fCz4H+CrDazFaYWRbRTs9d0/bZBfxS7PMPA0+7WC9DQM16zWa2Ffg7omEe9DbVuCtet3Ou3zlX5pxb7pxbTrTv4APOuSZvyk2Iufz7/g7Ru3PMrIxoE8wbqSwyweZyzaeBWwDMrJ5ooHeltMrU2wXcExvtcj3Q75xrn9eRvO4BnqV3+A6idyXHgd+Lfe0Pif5nhuiL/S9EF6h+GVjpdc0puOYfAOeAfbGPXV7XnIrrnrbvMwR8lMscX2sj2tR0GDgA/ILXNafgmhuA54iOgNkHvMfrmhNwzY8C7cAE0Xdd9wEPAA9c8jr/dexncuBa/m3r0X8RkTTh5yYXERG5Cgp0EZE0oUAXEUkTCnQRkTShQBcRSRMKdBGRNKFAFxFJE/8fP58+feyTWKEAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "nx = 164\n", "x = np.linspace(0, 1, nx, endpoint=False)\n", "dx = x[1] - x[0]\n", "xx = np.linspace(0, 1, 1000, endpoint=False)\n", "\n", "lmbda = 0.95\n", "nt = 250\n", "print('tsteps = %d' % nt)\n", "print(' dx = %g' % dx)\n", "print('lambda = %g' % lmbda)\n", "\n", "J = np.arange(0, nx) # all vertices\n", "Jm1 = np.roll(J, 1)\n", "Jp1 = np.roll(J, -1)\n", "\n", "plt.plot(x, g(x))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Plot the solution:" ] }, { "cell_type": "code", "execution_count": 21, "metadata": {}, "outputs": [], "source": [ "if 1:\n", " # Burgers\n", " def f(u):\n", " return u**2/2\n", " def fprime(u):\n", " return u\n", "else:\n", " # advection\n", " def f(u):\n", " return u\n", " def fprime(u):\n", " return 1+0*u\n", "\n", "steps_per_frame = 2" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Part I: Lax-Friedrichs\n", "\n", "Implement `rhs` for a Lax-Friedrichs flux:\n", "\n", "Recall (local) Lax-Friedrichs:\n", "$$\n", "f^{\\ast} (u^{_-}, u^+) = \\frac{f (u^-) + f (u^+)}{2} - \\frac{\\alpha}{2} (u^+ - u^-) \n", "\\quad\\text{with}\\quad\n", "\\alpha = \\max \\left( |f' (u^-)|, |f' (u^+)| \\right).$$\n", "Recall FV:\n", "$$ \\bar{u}_{j,\\ell+1} = \\bar{u}_{j,\\ell} - \\frac{h_t}{h_x} (f (u_{j + 1 / 2,\\ell}) -\n", " f (u_{j - 1 / 2,\\ell})) . $$\n" ] }, { "cell_type": "code", "execution_count": 29, "metadata": {}, "outputs": [], "source": [ "#clear\n", "def rhs(u):\n", " uplus = u[Jp1]\n", " uminus = u[J]\n", " alpha = np.maximum(np.abs(fprime(uplus)), np.abs(fprime(uminus)))\n", " # right-looking, between J and Jp1\n", " fluxes = (\n", " (f(uplus)+f(uminus))/2\n", " - alpha/2*(uplus-uminus)\n", " )\n", " return - 1/dx*(fluxes[J]-fluxes[Jm1])" ] }, { "cell_type": "code", "execution_count": 30, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\n", "\n", "\n", "\n", "\n", "\n", "
\n", " \n", "
\n", " \n", "
\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
\n", "
\n", " \n", " \n", " \n", " \n", " \n", " \n", "
\n", "
\n", "
\n", "\n", "\n", "\n" ], "text/plain": [ "" ] }, "execution_count": 30, "metadata": {}, "output_type": "execute_result" }, { "data": { "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "u = g(x)\n", "\n", "def timestepper(n):\n", " for i in range(steps_per_frame):\n", " dt = dx*lmbda/np.max(np.abs(u))\n", " u[:] = u + dt*rhs(u)\n", " \n", " line.set_data(x, u)\n", " return line\n", "\n", "fig = plt.figure(figsize=(5,5))\n", "line, = plt.plot(x, u)\n", "\n", "ani = animation.FuncAnimation(\n", " fig, timestepper,\n", " frames=nt//steps_per_frame,\n", " interval=30)\n", "html = HTML(ani.to_jshtml())\n", "plt.clf()\n", "html" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Part II: Second-Order Reconstruction\n", "\n", "First, need a second-order time integrator:" ] }, { "cell_type": "code", "execution_count": 24, "metadata": {}, "outputs": [], "source": [ "def rk2_step(dt, u, rhs):\n", " k1 = rhs(u)\n", " k2 = rhs(u+dt*k1)\n", " return u+0.5*dt*(k1+k2)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now upgrade the reconstruction to second order. \n", "\n", "**NOTE:** It's OK to assume (here!) that the wind blows from the right to simplify upwinding." ] }, { "cell_type": "code", "execution_count": 25, "metadata": {}, "outputs": [], "source": [ "#clear\n", "def rhs(u):\n", " # right-looking, between J and Jp1\n", " fluxes = f(u[J] + 1/2*(u[J]-u[Jm1]))\n", " return - 1/dx*(fluxes[J]-fluxes[Jm1])" ] }, { "cell_type": "code", "execution_count": 26, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\n", "\n", "\n", "\n", "\n", "\n", "
\n", " \n", "
\n", " \n", "
\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
\n", "
\n", " \n", " \n", " \n", " \n", " \n", " \n", "
\n", "
\n", "
\n", "\n", "\n", "\n" ], "text/plain": [ "" ] }, "execution_count": 26, "metadata": {}, "output_type": "execute_result" }, { "data": { "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "u = g(x)\n", "\n", "def timestepper(n):\n", " # to simplify upwinding\n", " assert np.min(u) >= 0\n", " \n", " for i in range(steps_per_frame):\n", " dt = 0.7*dx*lmbda/np.max(np.abs(u))\n", " u[:] = rk2_step(dt, u, rhs)\n", " \n", " line.set_data(x, u)\n", " return line\n", "\n", "fig = plt.figure(figsize=(5,5))\n", "line, = plt.plot(x, u)\n", "\n", "ani = animation.FuncAnimation(\n", " fig, timestepper,\n", " frames=nt//steps_per_frame,\n", " interval=30)\n", "html = HTML(ani.to_jshtml())\n", "plt.clf()\n", "html" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.8.3" } }, "nbformat": 4, "nbformat_minor": 4 }