amr_color_Rgen_spread.ipynb

{
 "metadata": {
  "name": "",
  "signature": "sha256:9a142035df11e8faad0044e60d60b7636844afec60e1a074bfe965fada51127a"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "% pylab inline"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Populating the interactive namespace from numpy and matplotlib\n"
       ]
      }
     ],
     "prompt_number": 396
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def q_exact(x):\n",
      "    #qtrue = exp(-100. * (x - 0.5)**2) * sin(80. * x)\n",
      "    qtrue = where(x < 0.1,  1., 0.2)\n",
      "    return qtrue"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 397
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Step 1: Initialize coarse grid, domain [0,1]\n",
      "# We have mx cells, indexed by 0 to mx-1, and x[i] is the mid-point of cell i\n",
      "mxc = 50\n",
      "dx = 1./mxc\n",
      "xc = linspace(dx/2, 1.-dx/2, mxc)  # computational grid\n",
      "\n",
      "tfinal = 1\n",
      "nsteps = 40\n",
      "dt = float(tfinal)/nsteps\n",
      "\n",
      "# intital condition\n",
      "#q0 = where(abs(xc) < 0.5,  1., 0.)\n",
      "#q0 = where(abs(xc-0.5) < 0.1,  1., 0.)\n",
      "#q0 = exp(-(xc-0.5)**2/(2))\n",
      "q0 = q_exact(xc)\n",
      "\n",
      "# Coefficient function u\n",
      "def u_fnc(x):\n",
      "    #u = where(abs(x) > 0.5,  exp(x), 1.)\n",
      "    #u = x**2\n",
      "    #u = where(abs(x) < 0.5,  x**2, 0.1)\n",
      "    n = numpy.array([x]).size\n",
      "    u = zeros(n)\n",
      "    for i in range(n):\n",
      "        u[i]= 0.75\n",
      "    \n",
      "    return u"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 398
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def q_true(u, x, time):\n",
      "    qtrue = q_exact(x-time*u)\n",
      "    return qtrue"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 399
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Step 2: Solve problem for ONE time step on coarse grid, for time step dt and 2*dt\n",
      "# u is the coefficient function in the form of a vector where u[i] = u(x) evaluated at x[i]\n",
      "# Similarly, q0 is a vector with the initial condition evaluated at cell mid-points x[i]\n",
      "# THE SOLUTION INCLUDES TWO GHOST CELLS ON EACH SIDE\n",
      "\n",
      "def color_onestep(u,q0,dt,dx):\n",
      "    \n",
      "    mx = len(q0)\n",
      "    qn = hstack([0,0, q0, 0,0])   # add two ghost cells on each end\n",
      "                                  # qn[2] through qn[mx+1]\n",
      "    \n",
      "    u = hstack([0,0,u,0,0])\n",
      "    u[1] = 2*u[2] - u[3]\n",
      "    u[0] = 2*u[1] - u[2]\n",
      "    u[mx+2] = 2*u[mx+1] - u[mx]\n",
      "    u[mx+3] = 2*u[mx+2] - u[mx+1]\n",
      "    \n",
      "    # Extrapolate boundary conditions (linear extrapolation)\n",
      "    qn[1] = 2*qn[2] - qn[3]\n",
      "    qn[0] = 2*qn[1] - qn[2]\n",
      "    qn[mx+2] = 2*qn[mx+1] - qn[mx]\n",
      "    qn[mx+3] = 2*qn[mx+2] - qn[mx+1]\n",
      "    \n",
      "    qnp = zeros(qn.shape)   # initialize array\n",
      "    flux = zeros(qn.shape)\n",
      "    \n",
      "    waves = zeros(mx+4)\n",
      "    for i in range(1,mx+4):\n",
      "        waves[i] = qn[i] - qn[i-1]\n",
      "        \n",
      "    theta = zeros(mx+4)    # smoothness indicator for limiter function\n",
      "    for i in range(1,mx+3):\n",
      "        if u[i] >= 0:\n",
      "            I = i-1\n",
      "        else: \n",
      "            I = i+1\n",
      "        if waves[i] <= 1e-6:\n",
      "            if waves[I] <= 1e-6:\n",
      "                theta[i] = 1\n",
      "        else:\n",
      "            theta[i] = waves[I]/waves[i]\n",
      "    \n",
      "    # implement minmod limiter by defining a vector phi_theta\n",
      "    phi_theta = zeros(mx+4)\n",
      "    \n",
      "    for i in range(1,mx+3):\n",
      "        # minmod\n",
      "        #phi_theta[i] = max(0.,min(1., theta[i]))\n",
      "        # no limiter\n",
      "        phi_theta[i] = 1\n",
      "        # superbee\n",
      "        #phi_theta[i] = max(0., min(1,2*theta[i]), min([2.,theta[i]]))\n",
      "        # van Leer\n",
      "        #phi_theta[i] = (theta[i]+ abs(theta[i]))/(1+abs(theta[i]))\n",
      "        \n",
      "    for i in range(1,mx+3):\n",
      "        flux[i] = 0.5*abs(u[i])*(1-(dt/dx)*abs(u[i]))*(phi_theta[i]*waves[i])\n",
      "    \n",
      "    for i in range(mx+3):\n",
      "        if u[i] > 0:\n",
      "            qnp[i] = qn[i] - (dt/dx) * u[i]*(waves[i])\n",
      "        else:\n",
      "            qnp[i] = qn[i] - (dt/dx) * u[i]*(waves[i+1])\n",
      "            \n",
      "        qnp[i] = qnp[i]-(dt/dx)*(flux[i+1]-flux[i])\n",
      "        \n",
      "    # Extrapolate boundary conditions (linear extrapolation)\n",
      "    qnp[mx+3] = 2*qnp[mx+2] - qnp[mx+1]\n",
      "    \n",
      "    return qnp"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 400
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Step 2: Solve problem for ONE time step on coarse grid, for time step dt and 2*dt\n",
      "# u is the coefficient function in the form of a vector where u[i] = u(x) evaluated at x[i]\n",
      "# Similarly, q0 is a vector with the initial condition evaluated at cell mid-points x[i]\n",
      "# THE SOLUTION INCLUDES TWO GHOST CELLS ON EACH SIDE\n",
      "\n",
      "def color_error(u,q0,dt,dx):\n",
      "    \n",
      "    mx = len(q0)\n",
      "    qn = hstack([0,0,0, q0, 0,0,0])   # add three ghost cells on each end\n",
      "                                      # qn[2] through qn[mx+1]\n",
      "        \n",
      "    #print max(abs(u))*dt/dx, 'cfl'\n",
      "    qnp = zeros(qn.shape)   # initialize array\n",
      "    flux = zeros(qn.shape)\n",
      "    \n",
      "    waves = zeros(mx+4)\n",
      "    for i in range(0,mx+4):\n",
      "        waves[i] = qn[i+2] - qn[i]\n",
      "        \n",
      "    for i in range(1,mx+1):\n",
      "        flux[i] = 0.5*abs(u[i-1])*(1-(dt/dx)*abs(u[i-1]))*waves[i-1]\n",
      "    \n",
      "    for i in range(3,mx+3):\n",
      "        if u[i-3] > 0:\n",
      "            qnp[i] = qn[i] - (dt/dx) * u[i-3]*(waves[i])\n",
      "        else:\n",
      "            qnp[i] = qn[i] - (dt/dx) * u[i-3]*(waves[i+1])\n",
      "            \n",
      "        qnp[i] = qnp[i]-(dt/dx)*(flux[i+1]-flux[i])\n",
      "        \n",
      "    # Extrapolate boundary conditions (linear extrapolation)\n",
      "    qnp[2] = 2*qnp[3] - qnp[4]\n",
      "    qnp[1] = 2*qnp[2] - qnp[3]\n",
      "    qnp[mx+3] = 2*qnp[mx+2] - qnp[mx+1]\n",
      "    qnp[mx+4] = 2*qnp[mx+3] - qnp[mx+2]\n",
      "    \n",
      "    return qnp[1:mx+5]"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 401
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def bilin_interp(x_l,x_r,y_b,y_t,ql_b,qr_b,ql_t,qr_t,in_x,in_y):\n",
      "    value = (1/((x_r-x_l)*(y_t-y_b)))*(ql_b*(x_r-in_x)*(y_t-in_y) + qr_b*(in_x-x_l)*(y_t-in_y) + ql_t*(x_r-in_x)*(in_y-y_b) + qr_t*(in_x-x_l)*(in_y-y_b))\n",
      "    return value"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 402
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def amr_color(u_fnc,q0,xc,tfinal,nsteps):\n",
      "    # q0 = initial condition evaluated on coarse grid \n",
      "    # u_fnc = a function name evaluating the coefficient function u at the grid points\n",
      "    # xc = coarse grid, domain [a,b] is divided into mx_c cells, and xc[i] is the mid-point of cell i\n",
      "        \n",
      "    dx = xc[1]- xc[0]             # length of cells in (given) coasrse grid\n",
      "    dt = float(tfinal)/nsteps     # length of time-step\n",
      "    \n",
      "    u = u_fnc(xc)\n",
      "    print max(abs(u))*dt/dx, 'cfl'\n",
      "    # u = velocity (coefficient) function, u[i] is cell-centered value in cell [i],\n",
      "    #     currently evaluated at cell-midpoint\n",
      "    \n",
      "    mx_c = len(xc)              # number of grid cells in coarse grid\n",
      "    \n",
      "    \n",
      "    \n",
      "    qn = hstack([0,0,q0,0,0])\n",
      "    qn[1] = 2*qn[2] - qn[3]\n",
      "    qn[0] = 2*qn[1] - qn[2]\n",
      "    qn[mx_c+2] = 2*qn[mx_c+1] - qn[mx_c]\n",
      "    qn[mx_c+3] = 2*qn[mx_c+2] - qn[mx_c+1]\n",
      "        \n",
      "    # Step 1: initalize coarse grid values q0 and u: Done before calling this function\n",
      "    # The rest of the steps are done over nsteps/K time steps, where K is the the period for regridding\n",
      "    for istep in range(nsteps):\n",
      "        \n",
      "        # Step 2: Solve problem on coarse grid for single time step\n",
      "        #         Done by calling a function 'color_onestep', which advances the solution on a grid\n",
      "        #         with spacing dx, over a single time step of length dt, returns solution with two ghost cell\n",
      "        #         values on each side\n",
      "        \n",
      "        \n",
      "        qnp = color_onestep(u,qn[2:mx_c+2],dt,dx) # with two ghost cell values on each side\n",
      "        \n",
      "    \n",
      "        # Step 3: Estimate error and flag cells based on Richardson extrapolation\n",
      "        #         We only allow for single level refinement.\n",
      "        # This estimation and grid generation procedure is only performed every K time steps. Buffer zone size and K\n",
      "        # are inter-related.\n",
      "        \n",
      "        p = 4\n",
      "        R = 2**1 # refinement factor, chosen as a power of 2\n",
      "        dx_hat = dx/R \n",
      "        dt_hat = dt/R\n",
      "        \n",
      "        # Compare errors between solving problem on coarse grid for 2 time steps of size 2*dt_hat \n",
      "        # and one step of size 4*dt_hat\n",
      "        qnew_D = color_onestep(u,qn[2:mx_c+2],dt,dx) # returned solutions had ghost cells\n",
      "        qnew_D = color_onestep(u,qnew_D[2:mx_c+2],dt,dx)\n",
      "        \n",
      "        qnew_D2h = color_error(u,qn[2:mx_c+2],2*dt,dx)\n",
      "\n",
      "        # Error\n",
      "        error = zeros(mx_c)\n",
      "        s = 2             # order of accuracy\n",
      "        flag_cnt = 0      # flagged cell count\n",
      "        tol = (dx**(s+2))/100\n",
      "\n",
      "        for i in range(0,mx_c):\n",
      "            error[i] = (qnew_D[i+2]-qnew_D2h[i+2])/(2**(s+1)-2)     \n",
      "            if abs(error[i]) >= tol:\n",
      "                flag_cnt = flag_cnt + 1\n",
      "        \n",
      "        # Use AMR only if any cells are flagged for refinement!\n",
      "        if flag_cnt > 0:\n",
      "            \n",
      "            # Cell flagging\n",
      "            flagged = zeros(flag_cnt)  # contains indices of cells flagged for refinement\n",
      "            \n",
      "            flag_cnt = 0\n",
      "            for i in range(0,mx_c):\n",
      "                if abs(error[i]) >= tol:\n",
      "                    flagged[flag_cnt] = i\n",
      "                    flag_cnt = flag_cnt + 1\n",
      "            \n",
      "            # cluster generation\n",
      "            bfs = 2 # buffer-size\n",
      "            gc = 1  # number of ghost cells\n",
      "            cls_tol = 2*(bfs+gc) # threshold for clustering\n",
      "    \n",
      "            temp = zeros(2*flag_cnt) # will be storing start and end points of clusters here\n",
      "            idx = 0\n",
      "            temp[idx] = max(flagged[0]-bfs,0.) # lower end of first cluster\n",
      "            clstr_cnt = 1\n",
      "\n",
      "            for i in range(0,flag_cnt-1):\n",
      "                if (flagged[i+1]-flagged[i]) > cls_tol:\n",
      "                    temp[idx+1] = flagged[i] + bfs  # upper end of (s/2 + 1)th cluster\n",
      "                    clstr_cnt = clstr_cnt + 1\n",
      "                    idx = idx+2\n",
      "                    temp[idx] = flagged[i+1] - bfs\n",
      "\n",
      "            temp[idx+1] = min(flagged[flag_cnt-1]+ bfs, mx_c-1)\n",
      "            clstr_cnt = int(clstr_cnt)\n",
      "\n",
      "            clusters = [0 for i in xrange(2*clstr_cnt)]\n",
      "\n",
      "            for i in range(0,2*clstr_cnt):\n",
      "                clusters[i] = int(temp[i])\n",
      "            #print clusters\n",
      "    \n",
      "            # Generate fine grid\n",
      "            \n",
      "            # Look for the largest sized cluster, to create array qhat of appropriate size\n",
      "            cellnum = [0 for i in xrange(clstr_cnt)]  # Array to store number of cells in each cluster, sans ghost cells.\n",
      "            \n",
      "            for i in range(0,clstr_cnt):\n",
      "                cellnum[i] = int(clusters[(2*i)+1] - clusters[2*i]) +1\n",
      "            maxclsize = max(cellnum)\n",
      "   \n",
      "            qhat = zeros([clstr_cnt, R*maxclsize+4])  # include two ghost cells (on fine grid) on either side of each cluster\n",
      "            \n",
      "            xhat = zeros([clstr_cnt, R*maxclsize+4])\n",
      "        \n",
      "            corr_l = zeros(clstr_cnt)   # store correction terms for fix up step\n",
      "            corr_r = zeros(clstr_cnt)\n",
      "            \n",
      "            # Create arrays for plotting coarse and fine grid values together\n",
      "            xplot = []\n",
      "            xplot.extend(xc[0:clusters[0]])\n",
      "            qplot = []\n",
      "            qplot.extend(qnp[2:clusters[0]+2])\n",
      "            \n",
      "            for iclstr in range(0,clstr_cnt):\n",
      "                \n",
      "                clidx = range(clusters[(2*iclstr)]-1, clusters[(2*iclstr)+1]+2)\n",
      "                # vector of cellnum + 2 coarse cell indices in cluster iclstr\n",
      "                # (could start at -1, if the first cell in the first cluster is 0, which is fine\n",
      "                # since we already have ghost cell values stored in qn and qnp; in fact, qn[1] and qnp[1]\n",
      "                # correspond to function values in cell '-1' on coarse grid. The value for cell j is stored in \n",
      "                # qn[j+2] and qnp[j+2])\n",
      "                \n",
      "                mxf_clstr = R*(cellnum[iclstr]) + 4   # number of fine grid cells in cluster \n",
      "                \n",
      "                # cell midpoints for fine grid\n",
      "                # Left fine grid ghost-cells (may be out of domain)\n",
      "                xhat[iclstr,0] = xc[clidx[1]] - (dx/(2*R))*(3+R)\n",
      "                xhat[iclstr,1] = xc[clidx[1]] - (dx/(2*R))*(1+R)\n",
      "                \n",
      "                # Right ghost cells\n",
      "                xhat[iclstr,mxf_clstr-2] = xc[clidx[cellnum[iclstr]]] + (dx/(2*R))*(R+1)\n",
      "                xhat[iclstr,mxf_clstr-1] = xc[clidx[cellnum[iclstr]]] + (dx/(2*R))*(R+3)\n",
      "                \n",
      "                # Interior fine grid cells\n",
      "                for j in range(0,cellnum[iclstr]):\n",
      "                    for icell in range(0,R):\n",
      "                        xhat[iclstr,R*j + 2 + icell] = xc[clidx[j+1]] + (dx/(2*R))*(2*icell + 1 - R)\n",
      "                        \n",
      "                \n",
      "                # Initializing fine grid q values (including ghost cells)\n",
      "                # Left ghost cells\n",
      "                slope1 = (qn[clidx[0]+2] - qn[clidx[0]+1])/dx\n",
      "                slope2 = (qn[clidx[0]+3] - qn[clidx[0]+2])/dx\n",
      "                if (slope1*slope2)<= 0:\n",
      "                    slope = 0\n",
      "                elif slope1 > 0:\n",
      "                    slope = min(slope1,slope2)\n",
      "                else:\n",
      "                    slope = max(slope1,slope2)\n",
      "                \n",
      "                qhat[iclstr, 0] = qn[clidx[0]+2] + slope*(xc[clidx[1]]-dx - xhat[iclstr, 0])\n",
      "                qhat[iclstr, 1] = qn[clidx[0]+2] + slope*(xc[clidx[1]]-dx - xhat[iclstr, 1])\n",
      "                \n",
      "                \n",
      "                # Right ghost cells\n",
      "                slope1 = (qn[clidx[cellnum[iclstr]]+3] - qn[clidx[cellnum[iclstr]]+2])/dx\n",
      "                slope2 = (qn[clidx[cellnum[iclstr]]+4] - qn[clidx[cellnum[iclstr]]+3])/dx\n",
      "                if (slope1*slope2)<= 0:\n",
      "                    slope = 0\n",
      "                elif slope1 > 0:\n",
      "                    slope = min(slope1,slope2)\n",
      "                else:\n",
      "                    slope = max(slope1,slope2)\n",
      "                \n",
      "                qhat[iclstr, mxf_clstr-2] = qn[clidx[cellnum[iclstr]]+3] + slope*(xc[clidx[cellnum[iclstr]]]+dx - xhat[iclstr, mxf_clstr-2])\n",
      "                qhat[iclstr, mxf_clstr-1] = qn[clidx[cellnum[iclstr]]+3] + slope*(xc[clidx[cellnum[iclstr]]]+dx - xhat[iclstr, mxf_clstr-1])\n",
      "                \n",
      "                # Interior fine grid cells\n",
      "                for j in range(0,cellnum[iclstr]):\n",
      "                    slope1 = (qn[clidx[j+1]+2] - qn[clidx[j+1]+1])/dx\n",
      "                    slope2 = (qn[clidx[j+1]+3] - qn[clidx[j+1]+2])/dx\n",
      "                    if (slope1*slope2)<= 0:\n",
      "                        slope = 0\n",
      "                    elif slope1 > 0:\n",
      "                        slope = min(slope1,slope2)\n",
      "                    else:\n",
      "                        slope = max(slope1,slope2)\n",
      "\n",
      "                    # minmod slope. can use other slopes, limiters, etc\n",
      "                    \n",
      "                    for icell in range(R):\n",
      "                        qhat[iclstr,R*j +2 +icell] = qn[clidx[j+1]+2] + slope*(xc[clidx[j+1]] - xhat[iclstr,R*j +2 +icell] )\n",
      "                    \n",
      "                uhat = u_fnc(xhat[iclstr,0:mxf_clstr])\n",
      "            \n",
      "                for tstep in range(0,R-1):\n",
      "                                    \n",
      "                    # compute correction for this fine grid time step\n",
      "                    # Fix-up in ghost cells to ensure conservation\n",
      "                    # Left ghost cell, indexed by clidx[0]. Look at the fine grid cell closest to the coarse-fine grid\n",
      "                    # interface, given by xhat[iclstr,1]\n",
      "                    \n",
      "                    u_lgc = uhat[1]\n",
      "                    # if u_lgc >= 0: # Right going correction, nothing to do\n",
      "                    if u_lgc < 0: # Left going correction\n",
      "                        corr_l[iclstr] = corr_l[iclstr] + (dt_hat/dx) * u_lgc*(qhat[iclstr,1] - qnp[clidx[0]+2]) \n",
      "            \n",
      "                    # Right ghost cell, indexed by clidx[cellnum[iclstr]+1]. Look at the fine grid cell closest to the \n",
      "                    # coarse-fine grid interface, given by xhat[iclstr,mxf_clstr-2]\n",
      "                    u_rgc = uhat[mxf_clstr-2]\n",
      "            \n",
      "                    if u_rgc >= 0: # Right going correction\n",
      "                        corr_r[iclstr] = corr_r[iclstr] + (dt_hat/dx) * u_rgc*(qnp[clidx[cellnum[iclstr]]+3]- qhat[iclstr,mxf_clstr-2])\n",
      "                    \n",
      "                    # Compute qhat for time t_n + dt_hat\n",
      "                    qhat[iclstr,0:mxf_clstr] = color_onestep(uhat,qhat[iclstr,0:mxf_clstr],dt_hat,dx_hat)[2:mxf_clstr+2]\n",
      "            \n",
      "                    \n",
      "                    # bilin_interp(x_l,x_r,y_b,y_t,ql_b,qr_b,ql_t,qr_t,in_x,in_y)\n",
      "                    # The way this is done depends on R. If R<=2 i.e. p < 2, need to look at two coarse grid ghost cells, \n",
      "                    # otherwise one suffices\n",
      "                    \n",
      "                    if R <= 2:\n",
      "                        qhat[iclstr,0] = bilin_interp(xc[clidx[1]]-2*dx, xc[clidx[1]]-dx, 0, dt, qn[clidx[0]+1],\n",
      "                                                      qn[clidx[0]+2], qnp[clidx[0]+1], qnp[clidx[0]+2],\n",
      "                                                      xhat[iclstr,0], (tstep+1)*dt_hat)\n",
      "                                       \n",
      "                        qhat[iclstr,mxf_clstr-1] = bilin_interp(xc[clidx[cellnum[iclstr]]]+dx,\n",
      "                                                                xc[clidx[cellnum[iclstr]]]+2*dx, 0, dt,\n",
      "                                                                qn[clidx[cellnum[iclstr]]+3], qn[clidx[cellnum[iclstr]]+4],\n",
      "                                                                qnp[clidx[cellnum[iclstr]]+3], qnp[clidx[cellnum[iclstr]]+4] ,\n",
      "                                                                xhat[iclstr,mxf_clstr-1], (tstep+1)*dt_hat)\n",
      "                    else:\n",
      "                        qhat[iclstr,0] = bilin_interp(xc[clidx[1]]-dx, xc[clidx[1]], 0, dt, qn[clidx[0]+2],\n",
      "                                                      qn[clidx[1]+2], qnp[clidx[0]+2], qnp[clidx[1]+2],\n",
      "                                                      xhat[iclstr,0], (tstep+1)*dt_hat)\n",
      "    \n",
      "                        qhat[iclstr,mxf_clstr-1] = bilin_interp(xc[clidx[cellnum[iclstr]]], xc[clidx[cellnum[iclstr]]]+dx,\n",
      "                                                                0, dt, qn[clidx[cellnum[iclstr]]+2],\n",
      "                                                                qn[clidx[cellnum[iclstr]]+3], qnp[clidx[cellnum[iclstr]]+2],\n",
      "                                                                qnp[clidx[cellnum[iclstr]]+3],\n",
      "                                                                xhat[iclstr,mxf_clstr-1], (tstep+1)*dt_hat)\n",
      "    \n",
      "    \n",
      "                    qhat[iclstr,1] = bilin_interp(xc[clidx[1]]-dx, xc[clidx[1]], 0, dt, qn[clidx[0]+2],\n",
      "                                                      qn[clidx[1]+2], qnp[clidx[0]+2], qnp[clidx[1]+2],\n",
      "                                                      xhat[iclstr,1], (tstep+1)*dt_hat)\n",
      "    \n",
      "                    qhat[iclstr,mxf_clstr-2] = bilin_interp(xc[clidx[cellnum[iclstr]]], xc[clidx[cellnum[iclstr]]]+dx,\n",
      "                                                                0, dt, qn[clidx[cellnum[iclstr]]+2],\n",
      "                                                                qn[clidx[cellnum[iclstr]]+3], qnp[clidx[cellnum[iclstr]]+2],\n",
      "                                                                qnp[clidx[cellnum[iclstr]]+3],\n",
      "                                                                xhat[iclstr,mxf_clstr-2], (tstep+1)*dt_hat)\n",
      "                    \n",
      "                    \n",
      "                # Solve over another time step dt_hat to get fine grid values for time t_n+1\n",
      "                qhat[iclstr,0:mxf_clstr] = color_onestep(uhat,qhat[iclstr,0:mxf_clstr],dt_hat,dx_hat)[2:mxf_clstr+2]  \n",
      "                    \n",
      "                # Add corrections to function values in the ghost cells at time t_(n+1)\n",
      "                qnp[clidx[0]+2] = qnp[clidx[0]+2] + corr_l[iclstr]\n",
      "                qnp[clidx[cellnum[iclstr]]+3] = qnp[clidx[cellnum[iclstr]]+3] + corr_r[iclstr]\n",
      "            \n",
      "                # Average to get coarse grid values within fine grid region\n",
      "                for k in range(0,cellnum[iclstr]):\n",
      "                    a = 0\n",
      "                    for icell in range(0,R):\n",
      "                        a = a + qhat[iclstr, R*k +icell + 2]\n",
      "                        \n",
      "                    qnp[clidx[k+1]+2] = a/R\n",
      "                \n",
      "                # Append fine grid values to plotting vectors\n",
      "                xplot.extend(xhat[iclstr, 0:mxf_clstr])\n",
      "                qplot.extend(qhat[iclstr, 0:mxf_clstr])\n",
      "                \n",
      "               \n",
      "        \n",
      "        \n",
      "        # updated values for next time step iteration\n",
      "        qn = qnp\n",
      "        \n",
      "        plt.figure(istep) \n",
      "        plot(xc,qn[3:mx_c+3],'b-')\n",
      "        \n",
      "    \n",
      "          \n",
      "        xplot.extend(xc[clusters[2*clstr_cnt-1]:mx_c])\n",
      "        #print xplot\n",
      "        #print len(xplot)\n",
      "        qplot.extend(qnp[clusters[2*clstr_cnt-1]+2:mx_c+2])\n",
      "        #print qplot\n",
      "        #print len(qplot)\n",
      "        \n",
      "        sort_index = numpy.argsort(xplot)\n",
      "        #xplot = sorted(xplot)\n",
      "        #qplot = [qplot[i] for i in sort_index]\n",
      "        \n",
      "        #plt.figure(istep) \n",
      "        plot(xplot,qplot,'b.')\n",
      "        qreal = q_exact(xc - u*(istep*dt))\n",
      "        plot(xc,qreal,'r-')\n",
      "        ylim([-0.5,1.5])\n",
      "\n",
      "       # plt.show()\n",
      "        \n",
      "        \n",
      "        \n",
      "    \n",
      "    return qn[3:mx_c+3]"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 403
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "#plot(xc,u_fnc(xc))\n",
      "qfinal = amr_color(u_fnc,q0,xc,tfinal,nsteps)\n",
      "plot(xc,qfinal)\n",
      "ylim([-0.5,1.5])"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "0.9375 cfl\n"
       ]
      },
      {
       "ename": "IndexError",
       "evalue": "list index out of range",
       "output_type": "pyerr",
       "traceback": [
        "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m\n\u001b[0;31mIndexError\u001b[0m                                Traceback (most recent call last)",
        "\u001b[0;32m<ipython-input-404-5cb7278143f0>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n\u001b[1;32m      1\u001b[0m \u001b[0;31m#plot(xc,u_fnc(xc))\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 2\u001b[0;31m \u001b[0mqfinal\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mamr_color\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mu_fnc\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mq0\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mxc\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mtfinal\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mnsteps\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m      3\u001b[0m \u001b[0mplot\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mxc\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mqfinal\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      4\u001b[0m \u001b[0mylim\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m-\u001b[0m\u001b[0;36m0.5\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;36m1.5\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
        "\u001b[0;32m<ipython-input-403-5cb1a64cf7bd>\u001b[0m in \u001b[0;36mamr_color\u001b[0;34m(u_fnc, q0, xc, tfinal, nsteps)\u001b[0m\n\u001b[1;32m    274\u001b[0m                 \u001b[0mqplot\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mextend\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mqhat\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0miclstr\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m0\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0mmxf_clstr\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    275\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 276\u001b[0;31m                 \u001b[0mxplot\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mextend\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mxc\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mclusters\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m2\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0miclstr\u001b[0m\u001b[0;34m+\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0mclusters\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m2\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0miclstr\u001b[0m\u001b[0;34m+\u001b[0m\u001b[0;36m2\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    277\u001b[0m                 \u001b[0mqplot\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mextend\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mqnp\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mclusters\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m2\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0miclstr\u001b[0m\u001b[0;34m+\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m+\u001b[0m\u001b[0;36m2\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0mclusters\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m2\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0miclstr\u001b[0m\u001b[0;34m+\u001b[0m\u001b[0;36m2\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m+\u001b[0m\u001b[0;36m2\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    278\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
        "\u001b[0;31mIndexError\u001b[0m: list index out of range"
       ]
      }
     ],
     "prompt_number": 404
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "plot(xc,q0)\n",
      "ylim(-0.5,1.5)\n",
      "plt.figure(2)\n",
      "plot(xc,u_fnc(xc))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def no_amr_color(u_fnc,q0,xc,tfinal,nsteps):\n",
      "    # q0 = initial condition evaluated on coarse grid (cell averages)\n",
      "    # u_fnc = a function name evaluating the coefficient function u at the grid points\n",
      "    # u = velocity (coefficient) function, u[i] is cell-centered value in cell [i],\n",
      "    #     currently evaluated at cell-midpoint\n",
      "    # xc = coarse grid, domain [a,b] is divided into mx cells, and xc[i] is the mid-point of cell i\n",
      "        \n",
      "    dx = xc[1]- xc[0]             # length of cells in (given) coasrse grid\n",
      "    dt = float(tfinal)/nsteps   # length of time-step\n",
      "    u = u_fnc(xc)\n",
      "    \n",
      "    mx_c = len(xc)              # number of grid cells in coarse grid\n",
      "    \n",
      "    # Step 1: initalize coarse grid values q0 and u: Done before calling this function\n",
      "    qn = q0\n",
      "    \n",
      "    for istep in range(nsteps):\n",
      "        \n",
      "    \n",
      "        # Step 2: Solve problem on coarse grid for single time step\n",
      "        #         Done by calling a function 'color_onestep', which advances the solution on a given grid\n",
      "        #         with spacing dx, over a single time step of length dt, returns solution with two ghost cell\n",
      "        #         values on each side\n",
      "        qnp = color_onestep(u,qn,dt,dx)\n",
      "        qn = qnp[2:mx_c+2]\n",
      "        \n",
      "        plt.figure(istep) \n",
      "        plot(xc,qn,'b.-')\n",
      "        qexact = q_true(u, xc, istep*dt)\n",
      "        plot(xc,qexact,'r-')\n",
      "        ylim([-0.5,1.5])\n",
      "        plt.show()\n",
      "        \n",
      "    return qn"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "qnew = no_amr_color(u_fnc,q0,xc,tfinal,nsteps)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": []
    }
   ],
   "metadata": {}
  }
 ]
}