{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Hughes model for pedestrian flow\n", "In 1D we can rewrite Hughes model as \n", "\\begin{align*}\n", "\\rho_t + \\nabla \\cdot (\\rho f(\\rho) \\textrm{sign} \\phi_x)_x &= 0\\\\\n", "\\lvert \\phi_x\\rvert &= \\frac{1}{f(\\rho)}\n", "\\end{align*}\n", "Note that the flux is discontinuous and we expect the formation of shocks at the discontinuity. \n", "\n", "### The 1D problem \n", "\n", "Let's consider the eikonal equation first\n", "\\begin{align*}\n", "\\lvert \\phi'\\rvert= 1 \\quad \\phi(0) = \\phi(1) = 0\n", "\\end{align*}\n", "\n", "* Let $\\phi:[-1,1] \\rightarrow \\mathbb{R}$ be a classic solution. Which theorem immediately gives a contradiction ?\n", "\n", "Numerical methods\n", "\n", "* Fast sweeping\n", "* Level set formulation\n", "* Fast marching\n", "\n", "We use a fast sweeping method:\n", "* H Zhao, *A fast sweeping method for eikonal equations*, Mathematics of Computation, 2005\n" ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "collapsed": true }, "outputs": [], "source": [ "import numpy as np\n", "import matplotlib.pyplot as plt\n", "from matplotlib import animation\n", "from IPython.display import HTML" ] }, { "cell_type": "code", "execution_count": 6, "metadata": { "collapsed": true }, "outputs": [], "source": [ "# Initialize solver\n", "npoints = 101 # number of gridpoints\n", "x = np.linspace(0,1.0,npoints) # mesh\n", "T = 1.0 # final time\n", "nt = 1000 # number of time steps\n", "dt = T/nt # discrete time step\n", "dx = 1.0/npoints # interval size\n", "lam = dt/dx\n", "beta = 1.0 # outflow rate\n", "rhomin = 0.5\n", "frhomin = rhomin * (1.0-rhomin) # maximum outflux\n", "flux = np.zeros(npoints+1)" ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "collapsed": true, "scrolled": true }, "outputs": [], "source": [ "def eikonal(npoints, costs):\n", " \n", " # Set the values at the boundary to zero, and inside the domain to a larger value\n", " phi = np.zeros(npoints)\n", " phi[1:-1] = 1.0\n", " it = 0\n", " dx = 1.0/npoints\n", " \n", " # Sweep from the left\n", " for i in np.arange(npoints-2):\n", " it = 1+i\n", " phiold = phi[it]\n", " \n", " phimin = np.minimum(phi[it-1], phi[it+1])\n", " phi[it] = np.minimum(phimin + costs[it] * dx, phiold)\n", "\n", " # Sweep from the right\n", " for i in np.arange(npoints-2):\n", " it = 1+i\n", " phiold = phi[npoints-1-it]\n", " phimin = np.minimum(phi[npoints-1 -it+1], phi[npoints - 1- it-1])\n", " phi[npoints-1 -it] = np.minimum(phimin + costs[npoints - 1 - it] * dx, phiold)\n", "\n", " phix = -np.gradient(phi,dx)\n", " return phi, phix\n" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAD8CAYAAACb4nSYAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xl8lNd97/HPbyQkAVpBKwKxSmIzO8Yri21s4rYQN2lj\nO06dxrdu0jpukzb3Jm1v0jq53dLGSRo3iZu4adLsbhbaOJeCzWI7XsAGbDYtLEZCoAWhDQHaTv+Y\nGVmWhRmkmXlm+b5fL17WzDyj+T0gf3XOec5zjjnnEBGR5ODzugAREYkehb6ISBJR6IuIJBGFvohI\nElHoi4gkEYW+iEgSUeiLiCQRhb6ISBJR6IuIJJFUrwsYLj8/382YMcPrMkRE4sorr7zS4pwruNJx\nMRf6M2bMYM+ePV6XISISV8zsjVCO0/COiEgSUeiLiCQRhb6ISBJR6IuIJBGFvohIElHoi4gkEYW+\niEgSUegnCOccP3j5JFVnOr0uRURiWMzdnCWjs7eujU/+5HUAVlcU8ODNs7ipPN/jqkQk1qilnyCq\nAy38D904k0MNHdz3zZfYdqjR46pEJNYo9BNEdWMXGeN8/MWvzeO5/7OOtFQfu0+0el2WiMQYDe8k\niJqmTuYUZuLzGRm+FOYUZHJE4/siMoxa+gmiurGTisKswceVxVlUNyr0ReStFPoJoP1CL40dl6go\nfjP0K4qyON1+kfYLvR5WJiKxRqGfAGoCLfqKoszB5+YGfgHUqLUvIkMo9BNAdWMXAOVDhneCrX6N\n64vIUAr9BFDd2MmEtBRKc8cPPjclJ4PM9FSN64vIWyj0E0BNUyflgZk7QWZGRVGm7tAVkbdQ6CeA\n6sYuyouy3vZ8ZXE21Y2dOOc8qEpEYpFCP861dffQ3HnpLRdxgyqLMjnX3Utz5yUPKhORWKTQj3OD\nF3FHaOkHL+ZWaVxfRAIU+nGuanC65gjDO4HnNK4vIkEK/ThX09hJZnoqU3Iy3vba5Mx08jPTNYNH\nRAYp9ONcdaN/zR0zG/H1ymLN4BGRNyn041xNY9eIF3GDKoqyqG7sYmBAM3hERKEf1852XeLs+Z4R\nx/ODKouyuNDbT/25C1GsTERiVUihb2YbzKzKzGrN7JMjvP5xMztkZq+Z2dNmNn3Ia/ebWU3gz/3h\nLD7Z1TRdfuZOUOXgcgwdUalJRGLbFUPfzFKAx4B3AfOBe8xs/rDD9gIrnHOLgCeBvw+8dxLwGWAV\ncC3wGTPLC1/5yW2khdaGq9AMHhEZIpSW/rVArXPumHOuB/gBsGnoAc657c657sDDF4Gpga/vALY6\n51qdc+eArcCG8JQuNU1dZKWnUpz99pk7QRPTU5k+eYIWXhMRILTQLwXqhjyuDzx3OQ8Avxzle+Uq\nVDd2Mqfo8jN3guYWZ3H4tIZ3RCS00B8pUUacCmJm9wErgM9fzXvN7EEz22Nme5qbm0MoScA/c6e8\n8PJDO0Fzi7M5fvY8F3r6o1CViMSyUEK/Hpg25PFUoGH4QWZ2G/DnwEbn3KWrea9z7nHn3Arn3IqC\ngoJQa09qoczcCZpXko1z6CYtEQkp9HcD5WY208zSgLuBzUMPMLOlwNfxB37TkJe2ALebWV7gAu7t\ngedkjEKZuRM0r8R/jIZ4RCT1Sgc45/rM7CH8YZ0CPOGcO2hmjwB7nHOb8Q/nZAI/Dowvn3TObXTO\ntZrZZ/H/4gB4xDnXGpEzSTLBmTuhDO9My5vAxLQUXcwVkSuHPoBz7ingqWHPfXrI17e9w3ufAJ4Y\nbYEyspqmLjLTUykZYc2d4Xw+o7I4i0Nq6YskPd2RG6eutObOcHNLsjlyukMbqogkOYV+nKpteuc1\nd4abV5JNx8U+GtovRrAqEYl1Cv041Hq+h5au0GbuBM0LLsegIR6RpKbQj0PBqZdzQriIGxRcg0cz\neESSm0I/DgWna15NSz8rYxzTJo3nsGbwiCQ1hX4cCu6WFcrMnaHmFmdreEckySn041BNY9dVzdwJ\nmleSzfGW81zs1XIMIslKoR+Hapo6r2rmTtC84iwGtByDSFJT6MeZ4Myd8sLQx/OD5pZkA7qYK5LM\nFPpxZnD5hVG09KdP8i/HcKhBoS+SrBT6caZ6FDN3gnw+Y/6UbA4o9EWSlkI/ztQ0dpI1ipk7QQum\n5HD4dAf9A1qOQSQZKfTjTKi7ZV3OginZdPf0c7zlfJgrE5F4oNCPMzWNXVSM4iJu0IIpOQAcbGgP\nV0kiEkcU+nEkuFvWaC7iBpUXZZKW4uOgxvVFkpJCP45UN47+Im7QuBQfc0uyOHBKLX2RZKTQjyM1\nTf7pmmMJffCP6x9s0Nr6IslIoR9Hqhs7ycpIpSg7fUzfZ8GUHNov9FJ/7kKYKhOReKHQjyPVjV1U\nFGWNeuZO0MJSXcwVSVYK/TjhnKOmcXRr7gw3tziLFJ/pYq5IElLox4mWrh7OdfeOas2d4TLGpTCn\nIFMXc0WSkEI/TgRXxhzrRdygBVqOQSQpKfTjxJuhP/bhHYAFpTk0d16iqUMbpYskE4V+nKhu7CJn\n/DgKssY2cydo4RT/Mssa1xdJLgr9OFHT2EllGGbuBM0PhL7G9UWSi0I/DjjnqG7sHNPyC8NlZYxj\nVsFE9tcr9EWSiUI/DjR1XqLjYl/YLuIGLZmay766Nt2ZK5JEFPpxoHoMu2W9k8XTcmnpusTpdl3M\nFUkWCv04EFxorTLMLf3F03IB2F/XFtbvKyKxS6EfB6rOdJCfmcbkzPDM3AmaV5LFuBRjX71CXyRZ\nKPTjQFVgzZ1wS09NYX5Jtlr6IklEoR/jBgb8a+5UFoc/9ME/xPN6fbv2zBVJEgr9GHeq7QLdPf1h\nH88PWjw1l/M9/Rxt7orI9xeR2KLQj3FHzgSWX4hgSx90MVckWSj0Y1y4F1obblb+RLLSU9mvi7ki\nSSGk0DezDWZWZWa1ZvbJEV5fbWavmlmfmb132Gv9ZrYv8GdzuApPFlVnOpmaN57M9NSIfH+fz1g0\nLYf9dbozVyQZXDH0zSwFeAx4FzAfuMfM5g877CTwQeB7I3yLC865JYE/G8dYb9KpOtMZsfH8oMVT\nczl8uoOLvf0R/RwR8V4oLf1rgVrn3DHnXA/wA2DT0AOccyecc68BAxGoMWn19A1wtLkrYuP5QYun\n5dI34Dh0WituiiS6UEK/FKgb8rg+8FyoMsxsj5m9aGbvvqrqktyJs+fpG3ARb+kvCVzM3XdS4/oi\niS6UgeKR1vK9mkndZc65BjObBTxjZq87546+5QPMHgQeBCgrK7uKb53YqgIzdyI1Rz+oKDuDkpwM\n9moGj0jCC6WlXw9MG/J4KtAQ6gc45xoC/z0G7ACWjnDM4865Fc65FQUFBaF+64RXdaaTFJ8xq2Bi\nxD9r+fQ8dh9v1YqbIgkulNDfDZSb2UwzSwPuBkKahWNmeWaWHvg6H7gRODTaYpNNVWMnM/Mnkp6a\nEvHPWjljEmc6LnKq7ULEP0tEvHPF0HfO9QEPAVuAw8CPnHMHzewRM9sIYGYrzawe+C3g62Z2MPD2\necAeM9sPbAf+1jmn0A9RdWPkZ+4ErZiRB8CeE+ei8nki4o2QJn87554Cnhr23KeHfL0b/7DP8Pf9\nCrhmjDUmpe6ePk62dvOeZW/7a42IucXZZKansvtEK+9eejXX6UUknuiO3BhV09iFc5G7E3e4FJ+x\nbHoer7yhlr5IIlPox6hozdwZauX0PKoaO2nv7o3aZ4pIdCn0Y9ThMx2MH5dC2aQJUfvM5TPycA5e\nPanWvkiiUujHqCOn/Wvop/hGuk0iMpZMyyXVZ+w+0Rq1zxSR6FLoxyDnHIfPdDCvJHpDOwAT0lJZ\nUJqjGTwiCUyhH4MaOy7R1t3LvJLsqH/2yul57Ktv41KfFl8TSUQK/Rh0OLDw2dzi6If+ihmT6Okb\n4MApLbUskogU+jHo8Bl/6Edz5k7Q8un+m7R2a4hHJCEp9GPQ4dOdlOaOJ2f8uKh/dkFWOrMKJvLS\nsbNR/2wRiTyFfgw6cjr6F3GHunF2Pi8db6WnT9sjiCQahX6Mudjbz7GW856M5wfdOGcy3T392jdX\nJAEp9GNMbVMX/QPOk5k7QdfPyscMnqtp8awGEYkMhX6MGZy54+HwTs6EcVxTmsPztQp9kUSj0I8x\nh093kjHOx4zJkd845Z3cOCeffXVtdF3q87QOEQkvhX6MOXKmg8qi6C6/MJKb5uTTN+B4+bhm8Ygk\nEoV+DHHOcfh0h6cXcYOWT88jLdXH87UKfZFEotCPIU2dlzjX3evpdM2gjHEprJyRp3F9kQSj0I8h\nb17E9b6lD3DD7HyOnOmkufOS16WISJgo9GPIoUDoz4uB4R3wj+sD/OqoWvsiiUKhH0MOnupg2qTx\n5EyI/vILI1lYmkN2Rqrm64skEIV+DDnQ0M6CkhyvyxiU4jNWVxSwvaqZgQHndTkiEgYK/RjRcbGX\nN852s7A0NoZ2gm6bV0RL1yUtySCSIBT6MeJQg388f0Fp7LT0AdZWFpDiM7YdbvS6FBEJA4V+jDgY\nDP0psdXSz52QxorpeWw71OR1KSISBgr9GHHwVDuFWekUZmV4XcrbrJ9fRFVjJ3Wt3V6XIiJjpNCP\nEQca2lkYY0M7QbfOKwLQEI9IAlDox4ALPf3UNnXF3NBO0Mz8icwumMjThzXEIxLvFPox4MiZDgYc\nLJgSmy19gNvmF/HisbN0XOz1uhQRGQOFfgw4EKMXcYe6bV4RfQOOXdXNXpciImOg0I8BhxrayRk/\njql5470u5bKWleUxaWIa/31Q4/oi8UyhHwMOnOpgYWk2Zt6uof9OUnzGHQuK2Ha4kQs9/V6XIyKj\npND3WG//AFVnOmN6PD9o4+JSunv62apZPCJxS6HvsZrGLnr6B2J6PD9o1cxJFGdnsHnfKa9LEZFR\nUuh77EBDOxDbM3eCfD5j45Ip7Khq5tz5Hq/LEZFRUOh77LX6NjLTU5mZ7+1G6KHauHgKfQOOpw6c\n9roUERmFkELfzDaYWZWZ1ZrZJ0d4fbWZvWpmfWb23mGv3W9mNYE/94er8ETxWn0715TmeL4ReqgW\nTMlmTmEmP9/X4HUpIjIKVwx9M0sBHgPeBcwH7jGz+cMOOwl8EPjesPdOAj4DrAKuBT5jZnljLzsx\nXOzt5/DpDhZPy/W6lJCZGZsWT+Hl4600tF3wuhwRuUqhtPSvBWqdc8eccz3AD4BNQw9wzp1wzr0G\nDAx77x3AVudcq3PuHLAV2BCGuhPC4dMd9PY7lkyL/fH8oTYtKQVg83619kXiTSihXwrUDXlcH3gu\nFGN5b8LbX+ffmCSeWvoAZZMnsKwslx/vqcM57aglEk9CCf2RBptD/T89pPea2YNmtsfM9jQ3J89t\n/vvr/cspF2fH3nLKV3LvqukcbT7PC0fPel2KiFyFUEK/Hpg25PFUINR+fUjvdc497pxb4ZxbUVBQ\nEOK3jn/769pYPC03pu/EvZxfX1RC3oRxfPuFN7wuRUSuQiihvxsoN7OZZpYG3A1sDvH7bwFuN7O8\nwAXc2wPPJb327l6OtZxnSZwN7QRljEvht1dOY+vhRk6364KuSLy4Yug75/qAh/CH9WHgR865g2b2\niJltBDCzlWZWD/wW8HUzOxh4byvwWfy/OHYDjwSeS3qvnQqM50+Nz9AHuG/VdAac43svnfS6FBEJ\nUWooBznnngKeGvbcp4d8vRv/0M1I730CeGIMNSak4EXca6bG18ydoaZNmsCtcwv5/ssneeiWOaSn\npnhdkohcge7I9ci+unZmFUwkZ/w4r0sZkw9cP4OWrh7+/4EzXpciIiFQ6HvAOce+ujaWxPHQTtDN\nc/KZmT+RJ54/oembInFAoe+B0+0Xaem6FHfz80fi8xkPrp7F/ro2dmpXLZGYp9D3QLzelHU571k2\nldLc8Ty6rUatfZEYp9D3wN66NtJSfMwryfK6lLBIS/Xx0C1z2F/Xxo4qtfZFYplC3wO7T7SyeFpO\nQs12Cbb2v7itWq19kRim0I+yCz39HDjVzvLpk7wuJazSUn189JY57K9vZ3tVk9fliMhlKPSjbH99\nG739jpUzEm+F6d9cNpWpeeP5hy3V9A+otS8SixT6UbbnhP+G5OXTEy/001J9/O8Nczl0uoPvvqQ1\neURikUI/yva8cY6KokxyJ6R5XUpE/MaiEm6cM5nPb6miqfOi1+WIyDAK/SjqH3C88sY5VsxIrPH8\nocyMRzYt5FLvAH/z1BGvyxGRYRT6UVTd2EnnxT5WJODQzlCzCzL5/TWz+OneU7x4TOvti8QShX4U\nBcfzVyZwSz/oD9fNYdqk8fzZT17n/KU+r8sRkQCFfhTtPnGOoux0puaN97qUiMsYl8LfvWcRx8+e\n59M/P+h1OSISoNCPoj0nWlkxY1Jc7pQ1GjfMzufhW8r5j1frefKVeq/LEREU+lFzqu0CDe0XWZng\n4/nDPXxrOdfNmsT//dkBaho7vS5HJOkp9KMkOJ6fyDN3RpLiM75091ImpKXwke++Slt3j9cliSQ1\nhX6UvHS8lcz0VOYWJ8Yia1ejKDuDr9y7jJNnu3ng3/Zwoaff65JEkpZCP0qer23hulmTSE1Jzr/y\n62dP5kt3L+HVk+f4g+++Qm//gNcliSSl5EygKKtr7eaNs93cMDvf61I89a5rSvjcuxeyvaqZP/3x\nfvoU/CJRF9LG6DI2vzraAsBN5ckd+gDvXzWdtu5ePr+lio4LvTz2/mVMSNOPoUi0qKUfBc/XnqUg\nK53ywkyvS4kJf7huDv/vroXsrG7mnsdfpKXrktcliSQNhX6EOef41dEWbpw9OWnm54fi/aum8/UP\nrKCqsZO7/vn5wS0kRSSyFPoRVtXYSUtXDzfO0dDOcOvnF/H937uO/n7He776Kx7fdZQBrcMvElEK\n/Qh7rsY/nq/QH9nSsjye+qObuXVeIX/91BHu/9eXOXm22+uyRBKWQj/Cnq9tYVb+RKbkJv56O6OV\nOyGNr923nM+9eyGvvnGO9Y/u5EvbarjYq/n8IuGm0I+g3v4BXjreyg1zJntdSswzM+67bjrb/mQN\nt80v4tFt1dz+6C5+urdeWy+KhJFCP4L21bXR3dPPTRraCVlJzngeu3cZ//7AKiakpfCxH+7n9kd3\nsnl/g+b1i4SBQj+CnqtpwQyum6WW/tW6qTyfpx6+mX9+/zJ8Zjz8/b2s/YcdfOPZY3Rc7PW6PJG4\npbtiImh7VRNLpuUm7H64kebzGXdeU8IdC4rZeqiRJ547zud+cZhHt1azcckU7l5ZxqKpOZoKK3IV\nFPoR0thxkdfq2/nEHZVelxL3UnzGhoXFbFhYzOv17fzbCyf42d4Gvv9yHXOLs7hraSkbl0yhJEcX\ny0WuxJyLrYtkK1ascHv27PG6jDH73ksn+bOfvs6WP15NZRKurBlpnRd7+c/9p/nRnjr21bVhBqtm\nTuLXFk1hw4JiCrLSvS5RJKrM7BXn3IorHqfQj4wPfWs3NU2d7PrEOg0/RNiJlvP8fF8DP99/imPN\n5/EZXDtzEhsWFHP7gmJNl5WkoND3UHdPH0sf2cq9q8r4zG8s8LqcpOGco7qxi1+8fppfvn6amqYu\nAK4pzWH9/CJum1fEvJIs/RKWhBRq6GtMPwKeq2nhUt8At80r8rqUpGJmVBZnUVmcxcfXV3CsuYst\nBxv570NneHRbNV/YWk1p7nhumVvILfMKuX7WZDLGpXhdtkhUhRT6ZrYB+BKQAnzDOfe3w15PB74N\nLAfOAu9zzp0wsxnAYaAqcOiLzrkPh6f02PX04SayMlK5dmZybY0Ya2YVZPKRtZl8ZO1smjovsv1I\nE1sPNfHkK/V858U3GD8uhRvnTGbd3EJumVuoC8GSFK4Y+maWAjwGrAfqgd1mttk5d2jIYQ8A55xz\nc8zsbuDvgPcFXjvqnFsS5rpj1sCA4+kjjaypKGBcku6SFYsKszJ438oy3reyjIu9/bx47CxPH27i\nmSNNbDvcBMDc4izWzS1kXWUhy8pyk3aXM0lsobT0rwVqnXPHAMzsB8AmYGjobwL+MvD1k8BXLEkH\nTvfVt9HS1cP6+RraiVUZ41JYW1nI2spCHnGOmqYunjnSxI6qJv5l1zG+uuMo2Rmp3FxRwLrKQtZU\nFGg2kCSMUEK/FKgb8rgeWHW5Y5xzfWbWDgRvQ51pZnuBDuAvnHPPjq3k2LbtUCMpPmNtRaHXpUgI\nzIyKoiwqirL48JrZdFzs5bmaFnZUNbG9qplfvHYagEVTc1hbUcDauYUsnppLii8p2zSSAEIJ/ZF+\nuodP+bncMaeBMufcWTNbDvzMzBY45zre8mazB4EHAcrKykIoKTY55/iv105z/azJ5EwY53U5MgrZ\nGeO485oS7rymhIEBx6HTHYO/AL6yvZYvP1NL3oRx3FxewNrKAtZUFDA5U70AiR+hhH49MG3I46lA\nw2WOqTezVCAHaHX++aCXAJxzr5jZUaACeMucTOfc48Dj4J+yOYrziAn76to42drNR2+Z43UpEgY+\nn7GwNIeFpTk8dEs5bd097Ar0AnZVN7N5fwNmsKg0hzWVhayrLGCRegES40IJ/d1AuZnNBE4BdwP3\nDjtmM3A/8ALwXuAZ55wzswL84d9vZrOAcuBY2KqPMT/f10Baqo87FhZ7XYpEQO6ENDYunsLGxVMY\nGHAcaGhnR1UzO6qa+MozNXz56RryJoxjdYW/F7C6XL0AiT1XDP3AGP1DwBb8UzafcM4dNLNHgD3O\nuc3AN4HvmFkt0Ir/FwPAauARM+sD+oEPO+daI3EiXuvrH+C/XmvgtnmFZGdoaCfR+XzGoqm5LJqa\ny8O3lnPufA+7aprZWdXMzupmfr4v0AuYmsvaigLWzS1kUWkOPvUCxGO6IzdMdlU38ztPvMzX7lvO\nBrX0k9rAgOP1U4FeQHUT++racA4mTUxjdXk+aysLWV1RwKSJWn1Vwkd35EbZz/adIisjlbWVBV6X\nIh7z+YzF03JZPC2XP7rtzV7AjqpmdlU387NAL2Dx1FzWVhawtlK9AIkehX4YXOztZ8uBM/z6oim6\nrV/eJm9iGpuWlLJpSelbegHbq5r40tM1fHFbDZMmprFmyLWAPPUCJEIU+mHw9OEmzvf0s2nJFK9L\nkRg3vBfQer6HZ2ua2X6kiZ3Vzfx07ynMYMm0XNZWFLK2soBr1AuQMFLoh8F/vFpPUXY6q7Qtolyl\nSUN6Af2DvQD/fQFffLqaR7dVMznQC1ijXoCEgUJ/jOpau9le1cRHbynX/GwZkxSfsWRaLkum5fLH\nt1VwtusSzw7eHdzET/aewhfsBVT6ewELp6gXIFdHoT9G//7iG/jMuPfa+L2TWGLT5Mx03r20lHcv\n9fcCXqtvG7wvILhUdH5mWuC+gEJWl+drP2a5IoX+GFzs7eeHe+q4Y0ERxTkZXpcjCSzFZywty2Np\nWR4fW+/vBeyqaWb7kWaeOdLET15VL0BCo9Afg837G2jr7uUD183wuhRJMpMz07lr6VTuWjqV/gHH\n/kAvYGdVE1/YOqQXUO5fJE69AAlS6I+Sc47vvPAGFUWZXDdLm6WId1J8xrKyPJaV5fHx9RW0dF1i\nV7X/voBnhl0LWBdYUnrBlGz1ApKUQn+U9tW18fqpdj67aYH2XJWYkp+Zzm8um8pvLvP3AvbVtbGz\nqokd1c3849Zq/nFrNfmZ6ayuyGddZSGrywu0KmwSUeiP0rd+dYLM9FTuWjbV61JELivFZyyfnsfy\n6Xl8/PbKt/YChlwLWFqWN7hG0PwS9QISmUJ/FI42d/Gf+xv4XzfPIjNdf4USP0LtBaypKGDd3AJu\nnqNeQKJRYo3CV56pJT01hQdXz/K6FJFRG94LaO4M9AKqm9l2uJH/eLUen8GysrzBNYIWTMnWcGac\n0yqbV+locxfrv7CT37t5Fp+6c57X5YhERF//wOCMoB1Vzbx+qh2Agqz0wTWCbi4vIGe8egGxQqts\nRsiXn65RK18SXmqKj+XTJ7F8+iT+JNAL2FntvzFs66FGnnylPjBryH9fwJqKAvUC4oRa+lehtqmL\n9Y/u5MHVs/jUu9TKl+Q0tBewvaqJA6f8W16rF+AttfQj4NFt1Ywfl8KDN6uVL8lreC+gqfMiu6pb\n2F7VxH8fPKNeQIxTSz9Ez9e28P5vvMQf3VrOx9ZXeF2OSEzq6x9gX13b4K5hQ3sBawNrBN1Unq9e\nQASE2tJX6Iegp2+ADV/aRf+AY8sfr9ZGKSIhauq8yM4q/4ygZ6ub6bjY5581VJbHmkr/UND8EvUC\nwkGhH0aPba/l81uq+NffXcm6ykKvyxGJS8FewPaqJnZUNXOwwd8LKMxKH5wSelN5PtkZ6gWMhkI/\nTOpau1n/6E7WVRby1fuWe12OSMJo6rjonxFU7d87uDPYC5geuC+gopB5JVnqBYRIoR8Gzjl+91u7\nefl4K9s+voYpueO9LkkkIfX1D7C3ro0dw3oBRdmBu4MrC7lRvYB3pNk7YfDN546zo6qZv9q4QIEv\nEkGpKT5WzpjEyhmT+MQdc2kM9AJ2VjXzywNn+NGeelJ9xrLpeYGVQguYW6xewGiopX8Ze0+e47e+\n9gK3zSviq/ct0w+XiEd6+wfYe/LNXsCh0/5eQHF2xuAaQTfOyScryXsBGt4Zg/buXu788rOYwS8e\nvlnTy0RiSGNHcEZQE89Wt9B5qY/UwWsBydsLUOiPUl//AB/+91fYUdXMkx+5gSXTcj2rRUTeWW//\nAK++cY4d1c1sP9LEkTOdQHL2AhT6o+Cc41M/eZ0f7K7jkU0L+J3rZ3hSh4iMzpn2i+ys9g8DPVfz\n9l7AurkFVBYlZi9AoT8Kn99yhMe2H+WhdXP40zsqPalBRMKjt3+AV944F1gp9M1eQElOxuAaQYnU\nC1DoX6VvPHuMz/3iMPdcW8Zf37UwIVsCIsnscr2AFTMCvYDKQiqKMuP2/32Ffoicc3xhazX/9Ewt\nd15TzD/ds4wUbRUnktDeqRcQvDv4xjn5cbUznkI/BL39A3zqJ6/z5Cv13HPtND67aSGpKb6ofLaI\nxI7T7Rf8M4KqmnmutoWuS32MSzFWTJ80+Esg1nsBCv0raOm6xMd+uI9na1r42G0VPHzrnJj+BxWR\n6Aj2ArYQrKPQAAAHEUlEQVRXNbGzqnmwFzAlJ4M1gSmhsdgLUOi/g+dqWvjYj/bRfqGXz717Ib+9\nYlpEP09E4lewF7C9qonna88O9gJWznizF1Be6H0vQKE/ggs9/XxxWzWPP3uM2QWZ/NM9S5lXkh2R\nzxKRxNPTF7wW4L8gXNXo7wWU5o5n9ZAZQV70AhT6Qzjn2Hqokb/6z0OcarvAPdeW8elfn8/4NK2L\nLyKj19B2YfBi8PO1LZzv6fesFxDW0DezDcCXgBTgG865vx32ejrwbWA5cBZ4n3PuROC1TwEPAP3A\nw865Le/0WeEO/YMN7Xx+SxU7qpqpKMrks5sWsmrW5LB9fxER8PcC9rzROjgUVN3YBfh7AWsqC1hb\n4e8FTIxQLyBsoW9mKUA1sB6oB3YD9zjnDg055g+ARc65D5vZ3cBdzrn3mdl84PvAtcAUYBtQ4Zzr\nv9znhSv0q8508sVt1fzywBmyM1J5+NZy7r9hBuM0O0dEouBUW3BG0Nt7AcGVQueEsRcQztC/HvhL\n59wdgcefAnDO/c2QY7YEjnnBzFKBM0AB8Mmhxw497nKfN5bQd87xbE0L33zuODurm8lMT+VDN83k\ngZtmatE0EfHMO/UCgsNAN8yePKZeQDjX0y8F6oY8rgdWXe4Y51yfmbUDkwPPvzjsvaUhfOZVq2vt\n5kPf2k1NUxf5mel8fH0FH7huOnkT0yLxcSIiIUtL9XHD7HxumJ3Pp+6cN9gL2F7VxM/2nuK7L50k\nLcXH7QuK+Mq9yyJaSyihP1LfY3j34HLHhPJezOxB4EGAsrKyEEp6u5KcDKZNmsDvr5nNbywuIT1V\nF2lFJDaV5o7n3lVl3LuqzN8LONHKjupmUqOwGkAooV8PDJ3IPhVouMwx9YHhnRygNcT34px7HHgc\n/MM7oRY/VGqKjyc+uHI0bxUR8Uxaqo8b5uRzw5z8qHxeKFc1dwPlZjbTzNKAu4HNw47ZDNwf+Pq9\nwDPOf7FgM3C3maWb2UygHHg5PKWLiMjVumJLPzBG/xCwBf+UzSeccwfN7BFgj3NuM/BN4DtmVou/\nhX934L0HzexHwCGgD/jDd5q5IyIikZUUN2eJiCS6UGfvaNK6iEgSUeiLiCQRhb6ISBJR6IuIJBGF\nvohIEom52Ttm1gy8MYZvkQ+0hKmceJFs55xs5ws652QxlnOe7pwruNJBMRf6Y2Vme0KZtpRIku2c\nk+18QeecLKJxzhreERFJIgp9EZEkkoih/7jXBXgg2c452c4XdM7JIuLnnHBj+iIicnmJ2NIXEZHL\niMvQN7MNZlZlZrVm9skRXk83sx8GXn/JzGZEv8rwCuGcP25mh8zsNTN72syme1FnOF3pnIcc914z\nc2YW9zM9QjlnM/vtwL/1QTP7XrRrDLcQfrbLzGy7me0N/Hzf6UWd4WJmT5hZk5kduMzrZmZfDvx9\nvGZm4d1KyzkXV3/wL+98FJgFpAH7gfnDjvkD4GuBr+8Gfuh13VE453XAhMDXH0mGcw4clwXswr8t\n5wqv647Cv3M5sBfICzwu9LruKJzz48BHAl/PB054XfcYz3k1sAw4cJnX7wR+iX/nweuAl8L5+fHY\n0r8WqHXOHXPO9QA/ADYNO2YT8G+Br58EbrVwbTnvjSues3Nuu3OuO/DwRfy7lMWzUP6dAT4L/D1w\nMZrFRUgo5/x7wGPOuXMAzrmmKNcYbqGcswOyA1/nMMLue/HEObcL/74jl7MJ+LbzexHINbOScH1+\nPIb+SBu1D99s/S0btQPBjdrjVSjnPNQD+FsK8eyK52xmS4Fpzrn/imZhERTKv3MFUGFmz5vZi2a2\nIWrVRUYo5/yXwH1mVg88BXw0OqV55mr/f78qoeyRG2vGslF7vAr5fMzsPmAFsCaiFUXeO56zmfmA\nR4EPRqugKAjl3zkV/xDPWvy9uWfNbKFzri3CtUVKKOd8D/At59w/mtn1+HfpW+icG4h8eZ6IaH7F\nY0v/ajZqZ9hG7fEqpA3mzew24M+Bjc65S1GqLVKudM5ZwEJgh5mdwD/2uTnOL+aG+rP9c+dcr3Pu\nOFCF/5dAvArlnB8AfgTgnHsByMC/Rk2iCun/99GKx9Afy0bt8eqK5xwY6vg6/sCP93FeuMI5O+fa\nnXP5zrkZzrkZ+K9jbHTOxfNem6H8bP8M/0V7zCwf/3DPsahWGV6hnPNJ4FYAM5uHP/Sbo1pldG0G\nficwi+c6oN05dzpc3zzuhnfcGDZqj1chnvPngUzgx4Fr1iedcxs9K3qMQjznhBLiOW8BbjezQ0A/\n8Ann3Fnvqh6bEM/5T4B/MbOP4R/m+GA8N+LM7Pv4h+fyA9cpPgOMA3DOfQ3/dYs7gVqgG/jdsH5+\nHP/diYjIVYrH4R0RERklhb6ISBJR6IuIJBGFvohIElHoi4gkEYW+iEgSUeiLiCQRhb6ISBL5H0u7\nynr1p8hDAAAAAElFTkSuQmCC\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "cost = 0.1 + 2.0 * np.exp(-(x-0.2)**2/0.01)#np.ones(npoints) #\n", "phi,b = eikonal(npoints, cost)\n", "plt.plot(x,phi)\n", "plt.show()\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Nonlinear conservation laws with discontinuous flux\n", "\n", "Such problems were analysed and solved by Towers and co-workers in a series of papers. In particular they considered a 1d problem\n", "\n", "\\begin{align*}\n", "\\partial_t u (x,t) = (k(x) f(u))_x\n", "\\end{align*}\n", "where $k$ is stationary but can have jump discontinuities.\n", "The developed numerical schemes have the form\n", "\\begin{align*}\n", "u^{n+1}_j = u^n_j + \\lambda (k_{j + \\frac{1}{2}} h_{j + \\frac{1}{2}} - k_{j - \\frac{1}{2}} h_{j - \\frac{1}{2}}) \n", "\\end{align*}\n", "where $\\lambda = \\frac{\\Delta t}{\\Delta x}$ and the arguments of the numerical flux have to be transposed if the flux changes sign (to ensure monotonicity of the scheme):\n", "\\begin{align*}\n", "h_{j+\\frac{1}{2}} &=\n", "\\begin{cases}\n", "h(u_{j+1},u_j) \\text{ if } k_{j+\\frac{1}{2}} \\geq 0\\\\\n", "h(u_j,u_{j+1}) \\text{ if } k_{j+\\frac{1}{2}} < 0.\n", "\\end{cases}\n", "\\end{align*}\n", "Different numerical fluxes were analysed - Godunov or EO flux. For further information see\n", "* J Towers, *Convergence of a finite difference scheme for conservation laws with a discontinuous flux*, SIAM J Numer. Anal., 38(2), 2000" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "data": { "text/html": [ "" ], "text/plain": [ "" ] }, "execution_count": 10, "metadata": {}, "output_type": "execute_result" } ], "source": [ "def hughes(self, rho, rhoold, npoints):\n", " \n", " # Solve eikonal equation in every time step\n", " [phi2,b] = eikonal(npoints,1.0/(1.0-rho))\n", "\n", " # Determine where the flux is positive and negative\n", " bmid = 0.5*(b[1:] + b[:-1])\n", " bplus = np.array(1.0 * (bmid >0.0))\n", " bmin = np.array(1.0 * (bmid <0.0))\n", " bcalc = np.concatenate([[0],bmid,[0]])\n", " \n", " rho1 = np.maximum(rho,rhomin)\n", " rho2 = np.minimum(rho,rhomin)\n", " \n", " f = rho[:] * (1.0-rho[:])\n", " \n", " # Godunov flux - calculate the solution to the corresponding Riemann problem\n", " nextisbigger = 1.0 * (rho[1:]>rho[:-1])\n", " fmin = np.minimum(f[1:],f[:-1])\n", " fmax = np.maximum(f[1:],f[:-1])\n", "\n", " flux[1:-1] = bplus[:] * (nextisbigger[:] * fmin + (1.0-nextisbigger[:]) * fmax) + (1.0-bplus[:]) * (nextisbigger[:] * fmax + (1.0-nextisbigger[:]) * fmin)\n", " \n", " rho[0] = rho[0] - lam * beta * rho[0]\n", " rho[1:-1] = rho[1:-1] - lam * (flux[2:-1] * bcalc[2:-1] - flux[1:-2] * bcalc[1:-2])\n", " rho[npoints-1] = rho[npoints-1] - lam * beta * rho[npoints-1]\n", " line.set_data(x,rho)\n", " \n", " rhoold[:] = rho\n", " line.set_data(x,rho)\n", "\n", "# Set initial datum - solver works only for initial data smaller than 0.5 correctly, for larger values the boundary\n", "# conditions have to be implemented correctly.\n", "\n", "rhoinit = 0.5\n", "rho = rhoinit * np.ones(npoints)\n", "rhoprev = rhoinit * np.ones(npoints)\n", "# Set figure properties for the movie\n", "fig1, ax1 = plt.subplots(1,1)\n", "ax1.set_xlim(0, 1)\n", "ax1.set_ylim(0, 0.5)\n", "\n", "line, = ax1.plot([],[],lw=2,color=\"b\")\n", "anim2 = animation.FuncAnimation(fig1, hughes, fargs=(rho,rhoprev,npoints), frames=200, interval=20, blit=False)\n", "HTML(anim2.to_html5_video())" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.6.3" } }, "nbformat": 4, "nbformat_minor": 2 }