Modeling traffic flow with the Lighthill-Whitham-Richards model in PyClaw.

pyclaw_example_traffic.ipynb

{
 "metadata": {
  "name": "",
  "signature": "sha256:043169b5892cf2f652ea9fbc33d62fe344f2a65754105c6b2ff7b08a70c7ab1d"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import matplotlib.pyplot as plt\n",
      "import time\n",
      "from IPython.display import clear_output\n",
      "import numpy as np"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Here is an example python script for solving the LWR traffic model:\n",
      "\n",
      "$$q_t + u_\\text{max}( q(1-q) )_x = 0. $$"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def traffic(use_petsc=0,solver_type='classic'):\n",
      "    from clawpack import riemann\n",
      "\n",
      "    if use_petsc:\n",
      "        import clawpack.petclaw as pyclaw\n",
      "    else:\n",
      "        from clawpack import pyclaw\n",
      "\n",
      "    #===========================================================================\n",
      "    # Set up solver and solver parameters\n",
      "    #===========================================================================\n",
      "    if solver_type=='sharpclaw':\n",
      "        solver = pyclaw.SharpClawSolver1D(riemann.traffic_1D)\n",
      "    else:\n",
      "        solver = pyclaw.ClawSolver1D(riemann.traffic_1D)\n",
      "        \n",
      "    solver.num_waves = 1\n",
      "    solver.bc_lower[0] = pyclaw.BC.periodic\n",
      "    solver.bc_upper[0] = pyclaw.BC.periodic\n",
      "\n",
      "    #======================================\n",
      "    # Set up domain and initialize solution\n",
      "    #======================================\n",
      "    x = pyclaw.Dimension('x',-30.0,30.0,500)\n",
      "    domain = pyclaw.Domain(x)\n",
      "    num_eqn = 1\n",
      "    state = pyclaw.State(domain,num_eqn)\n",
      "\n",
      "    grid = state.grid\n",
      "    xc=grid.x.centers\n",
      "\n",
      "    # Stopped traffic at a light:\n",
      "    #state.q[0,:] = 0.75*(xc<0) + 0.1*(xc>0.)\n",
      "    \n",
      "    # Formation of a traffic jam\n",
      "    state.q[0,:] = 0.25 + 0.25*np.exp(-0.01*xc**2)\n",
      "\n",
      "    state.problem_data['efix']=True\n",
      "    state.problem_data['umax']=1.\n",
      "\n",
      "    #===========================================================================\n",
      "    # Setup controller and controller parameters. Then solve the problem\n",
      "    #===========================================================================\n",
      "    claw = pyclaw.Controller()\n",
      "    claw.tfinal = 50.0\n",
      "    claw.solution = pyclaw.Solution(state,domain)\n",
      "    claw.solver = solver\n",
      "    claw.num_output_times = 25\n",
      "    claw.keep_copy = True\n",
      "\n",
      "    status = claw.run()\n",
      "\n",
      "    return claw"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "claw = traffic()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%pylab inline\n",
      "from matplotlib import animation\n",
      "import matplotlib.pyplot as plt\n",
      "from clawpack.visclaw.JSAnimation import IPython_display\n",
      "\n",
      "xc=claw.solution.state.grid.x.centers\n",
      "\n",
      "fig = plt.figure(figsize=[8,4])\n",
      "ax = plt.axes(xlim=(xc[0], xc[-1]), ylim=(0, 1.0))\n",
      "line, = ax.plot([], [], lw=1)\n",
      "\n",
      "def fplot(i):\n",
      "    frame = claw.frames[i]\n",
      "    q = frame.q[0,:]\n",
      "    line.set_data(xc, q)\n",
      "    return line,\n",
      "\n",
      "animation.FuncAnimation(fig, fplot, frames=len(claw.frames))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    }
   ],
   "metadata": {}
  }
 ]
}