# Set up controller and controller parameters
    #===========================================================================
    claw = pyclaw.Controller()
    claw.tfinal = 2.5
    claw.solution = pyclaw.Solution(state,domain)
    claw.solver = solver
    claw.outdir = outdir
    claw.num_output_times = 10

    #===========================================================================
    # Solve the problem
    #===========================================================================
    status = claw.run()

    #===========================================================================
    # Plot results
    #===========================================================================
    if iplot:     plot.interactive_plot(outdir=outdir,format=claw.output_format)
    if htmlplot:  plot.html_plot(outdir=outdir,format=claw.output_format)


if __name__=="__main__":
    from pyclaw.util import run_app_from_main
    output = run_app_from_main(shallow2D)
    print 'Error: ', output

            a2=1.0 #a2values[ii]
            KB = Z
            rhoB = Z
            state.aux_global['KB'] = KB
            state.aux_global['rhoB'] = rhoB
            state.aux_global['trdone'] = False
            state.aux=setaux(xc,rhoB,KB,rhoA,KA,alpha,bc_lower,xupper=xupper)
            grid.x.bc_lower=2
            grid.x.bc_upper=2
            state.t = 0.0
            qinit(state,ic=2,a2=a2)
            init_solution = Solution(state)
            claw.solution = init_solution
            claw.solution.t = 0.0

            claw.tfinal = tfinal
            claw.outdir = './_output_Z'+str(Z)+'_'+str(cellsperlayer)
            status = claw.run()

    else:
        # Solve
        status = claw.run()

    if htmlplot:  pyclaw.plot.html_plot(outdir=outdir)
    if iplot:     pyclaw.plot.interactive_plot(outdir=outdir)


if __name__=="__main__":
    from pyclaw.util import run_app_from_main
    output = run_app_from_main(stegoton)

    #qinit(state,mx,my)


    #===========================================================================
    # Set up controller and controller parameters
    #===========================================================================
    claw = pyclaw.Controller()
    claw.keep_copy = False
    claw.output_style = 1
    claw.num_output_times = 10
    claw.tfinal = 10
    claw.solution = pyclaw.Solution(state,domain)
    claw.solver = solver
    claw.outdir = outdir

    #===========================================================================
    # Solve the problem
    #===========================================================================
    status = claw.run()

    #===========================================================================
    # Plot results
    #===========================================================================
    if htmlplot:  pyclaw.plot.html_plot(outdir=outdir)
    if iplot:     pyclaw.plot.interactive_plot(outdir=outdir)


if __name__=="__main__":
    from pyclaw.util import run_app_from_main
    output = run_app_from_main(shallow_4_Rossby_Haurwitz)

    claw.keep_copy = False
    claw.outstyle = 1
    claw.nout = 10
    claw.tfinal = 1.0
    claw.solution = pyclaw.Solution(state)
    claw.solver = solver
    claw.outdir = outdir

    #===========================================================================
    # Solve the problem
    #===========================================================================
    status = claw.run()

    #===========================================================================
    # Plot results
    #===========================================================================
    if htmlplot:  pyclaw.plot.html_plot(outdir=outdir)
    if iplot:     pyclaw.plot.interactive_plot(outdir=outdir)



if __name__=="__main__":
    from pyclaw.util import run_app_from_main
    output = run_app_from_main(advection_annulus)
    print 'Error: ',output

    state = pyclaw.State(grid,meqn)

    state.aux_global['gamma']= gamma
    state.aux_global['gamma1']= gamma1

    state.q[0,:] = 1.
    state.q[1,:] = 0.
    x =state.grid.x.center
    state.q[2,:] = ( (x<0.1)*1.e3 + (0.1<=x)*(x<0.9)*1.e-2 + (0.9<=x)*1.e2 ) / gamma1

    solver.limiters = 4

    claw = pyclaw.Controller()
    claw.tfinal = 0.038
    claw.solution = pyclaw.Solution(state)
    claw.solver = solver
    claw.nout = 10
    claw.outdir = outdir

    # Solve
    status = claw.run()

    if htmlplot:  pyclaw.plot.html_plot(outdir=outdir)
    if iplot:     pyclaw.plot.interactive_plot(outdir=outdir)

    return claw.solution.q

if __name__=="__main__":
    from pyclaw.util import run_app_from_main
    output = run_app_from_main(wcblast)

    meqn = 5
    maux=8
    state = pyclaw.State(grid,meqn,maux)

    state.aux_global['gamma']= gamma
    state.aux_global['gamma1']= gamma1

    qinit(state)
    auxinit(state)

    solver.user_bc_lower=shockbc

    claw = pyclaw.Controller()
    claw.tfinal = 0.75
    claw.solution = pyclaw.Solution(state)
    claw.solver = solver
    claw.nout = 10
    claw.outdir = outdir

    # Solve
    status = claw.run()

    if htmlplot:  pyclaw.plot.html_plot(outdir=outdir)
    if iplot:     pyclaw.plot.interactive_plot(outdir=outdir)

    return claw.solution.q

if __name__=="__main__":
    from pyclaw.util import run_app_from_main
    output = run_app_from_main(shockbubble)

    # ================
    qinit(state) # This function is defined above


    #===========================================================================
    # Set up controller and controller parameters
    #===========================================================================
    claw = pyclaw.Controller()
    claw.tfinal = 2.0
    claw.solution = pyclaw.Solution(state,domain)
    claw.solver = solver
    claw.num_output_times = 10
    print claw.num_output_times
    claw.outdir = outdir

    #===========================================================================
    # Solve the problem
    #===========================================================================
    status = claw.run()

    #===========================================================================
    # Plot results
    #===========================================================================
    if htmlplot:  pyclaw.plot.html_plot(outdir=outdir)
    if iplot:     pyclaw.plot.interactive_plot(outdir=outdir)


if __name__=="__main__":
    from pyclaw.util import run_app_from_main
    output = run_app_from_main(advection2D)

        ts.setTimeStep(solver.dt)
        ts.setDuration(claw.tfinal, max_steps=1000)
        ts.setFromOptions()
        q0 = x[:].copy()
        ts.solve(x)
        qfinal = x[:].copy()
        dx = grid.d[0]
    else:
        status = claw.run()
        #This test is set up so that the waves pass through the domain
        #exactly once, and the final solution should be equal to the
        #initial condition.  Here we output the 1-norm of their difference.
        q0=claw.frames[0].grid.gqVec.getArray().reshape([-1])
        qfinal=claw.frames[claw.nout].grid.gqVec.getArray().reshape([-1])
        dx=claw.frames[0].grid.d[0]

    if htmlplot:  plot.html_plot(outdir=outdir,format=output_format)
    if iplot:     plot.interactive_plot(outdir=outdir,format=output_format)
    if petscPlot: plot.plotPetsc(output_object)

    print 'Max error:', np.max(qfinal - q0)

    return dx*np.sum(np.abs(qfinal-q0))


if __name__=="__main__":
    from pyclaw.util import run_app_from_main
    error = run_app_from_main(acoustics)
    print 'Error: ',error

    #===========================================================================
    x = pyclaw.Dimension('x',0.0,1.0,500)
    grid = pyclaw.Grid(x)
    meqn = 1
    state = pyclaw.State(grid,meqn)

    xc=grid.x.center
    state.q[0,:] = np.sin(np.pi*2*xc) + 0.50
    state.aux_global['efix']=True

    #===========================================================================
    # Setup controller and controller parameters. Then solve the problem
    #===========================================================================
    claw = pyclaw.Controller()
    claw.tfinal =0.5
    claw.solution = pyclaw.Solution(state)
    claw.solver = solver
    claw.outdir = outdir

    status = claw.run()

    if htmlplot:  pyclaw.plot.html_plot(outdir=outdir)
    if iplot:     pyclaw.plot.interactive_plot(outdir=outdir)


if __name__=="__main__":
    from pyclaw.util import run_app_from_main
    output = run_app_from_main(burgers)
    print 'Error: ',output

    solver.cfl_desired = 0.8

    x = pyclaw.Dimension('x',0.0,1.0,mx)
    grid = pyclaw.Grid(x)
    num_eqn = 1
    state = pyclaw.State(grid,num_eqn)
    state.problem_data['u']=1.

    xc=grid.x.center
    if IC=='gauss_square':
        state.q[0,:] = np.exp(-beta * (xc-x0)**2) + (xc>0.6)*(xc<0.8)
    elif IC=='wavepacket':
        state.q[0,:] = np.exp(-beta * (xc-x0)**2) * np.sin(80.*xc)
    else:
        raise Exception('Unrecognized initial condition specification.')

    claw = pyclaw.Controller()
    claw.solution = pyclaw.Solution(state)
    claw.solver = solver
    claw.outdir = outdir

    claw.tfinal =10.0
    status = claw.run()

    if htmlplot:  pyclaw.plot.html_plot(outdir=outdir)
    if iplot:     pyclaw.plot.interactive_plot(outdir=outdir)

if __name__=="__main__":
    from pyclaw.util import run_app_from_main
    output = run_app_from_main(fig_61_62_63)

    solver.limiters = pyclaw.limiters.tvd.MC
    solver.bc_lower[0] = 2
    solver.bc_upper[0] = 2
    solver.aux_bc_lower[0] = 2
    solver.aux_bc_upper[0] = 2

    xlower=0.0; xupper=1.0; mx=100
    x    = pyclaw.Dimension('x',xlower,xupper,mx)
    grid = pyclaw.Grid(x)
    maux=1
    meqn = 1
    state = pyclaw.State(grid,meqn,maux)

    qinit(state)
    auxinit(state)

    claw = pyclaw.Controller()
    claw.outdir = outdir
    claw.solution = pyclaw.Solution(state)
    claw.solver = solver

    claw.tfinal = 1.0
    status = claw.run()

    if htmlplot:  pyclaw.plot.html_plot(outdir=outdir)
    if iplot:     pyclaw.plot.interactive_plot(outdir=outdir)

if __name__=="__main__":
    from pyclaw.util import run_app_from_main
    output = run_app_from_main(vc_advection)

    state.aux_global['cc']=np.sqrt(bulk/rho)

    #========================================================================
    # Set the initial condition
    #========================================================================
    xc=grid.x.center
    beta=100; gamma=0; x0=0.75
    state.q[0,:] = 0.5*np.exp(-80 * xc**2) + 0.5*(np.abs(xc+0.2)<0.1)
    state.q[1,:] = 0.
    
    #========================================================================
    # Set up the controller object
    #========================================================================
    claw = pyclaw.Controller()
    claw.solution = pyclaw.Solution(state)
    claw.solver = solver
    claw.outdir = outdir
    claw.tfinal = 3.0
    claw.nout   = 30

    # Solve
    status = claw.run()

    # Plot results
    if htmlplot:  pyclaw.plot.html_plot(outdir=outdir)
    if iplot:     pyclaw.plot.interactive_plot(outdir=outdir)

if __name__=="__main__":
    from pyclaw.util import run_app_from_main
    output = run_app_from_main(fig_31_38)

    x = pyclaw.Dimension('x',-2.0,2.0,mx)
    y = pyclaw.Dimension('y',-2.0,2.0,my)
    domain = pyclaw.Domain([x,y])
    num_eqn = 1
    state = pyclaw.State(domain,num_eqn)

    qinit(state)

    solver.dimensional_split = 1
    solver.cfl_max = 1.0
    solver.cfl_desired = 0.9
    solver.num_waves = 2
    solver.limiters = pyclaw.limiters.tvd.minmod

    claw = pyclaw.Controller()
    claw.tfinal = 1.0
    claw.solution = pyclaw.Solution(state,domain)
    claw.solver = solver
    claw.num_output_times = 10

    # Solve
    status = claw.run()

    if htmlplot:  pyclaw.plot.html_plot(outdir=outdir)
    if iplot:     pyclaw.plot.interactive_plot(outdir=outdir)


if __name__=="__main__":
    from pyclaw.util import run_app_from_main
    output = run_app_from_main(kpp)

        state.aux_global = {}
        state.aux_global['turnZero_half_2D'] = turnZero_half_2D
        state.aux_global['t_turnZero'] = t_turnZero
        state.mp = 1

        state.aux = setaux(grid.x.center,grid.y.center)
        #Initial condition
        qinit(state,A,x0,y0,varx,vary)

        claw.solution = pyclaw.Solution(state)
        claw.nout = nout

    #claw.p_function = p_function
    claw.compute_F = compute_F
    grid.add_gauges([[0.25,0.25],[0.75,0.25],[0.25,0.75],[0.75,0.75]])
    solver.compute_gauge_values = gauge_pfunction
    claw.write_aux_init = False

    #Solve
    status = claw.run()
    
    #strain=claw.frames[claw.nout].state.gqVec.getArray().reshape([grid.ng[0],grid.ng[1],meqn])[:,:,0]
    #return strain

    if iplot:    pyclaw.plot.interactive_plot()
    if htmlplot: pyclaw.plot.html_plot()

if __name__=="__main__":
    from pyclaw.util import run_app_from_main
    output = run_app_from_main(psystem2D)

    x0 = 0.75
    state.q[0,:] = np.exp(-beta * (xc-x0)**2) * np.cos(gamma * (xc - x0))

    #===========================================================================
    # Set up controller and controller parameters
    #===========================================================================
    claw = pyclaw.Controller()
    claw.solution = pyclaw.Solution(state)
    claw.solver = solver
    claw.outdir = outdir

    claw.tfinal = 1.0

    #===========================================================================
    # Solve the problem
    #===========================================================================
    status = claw.run()

    #===========================================================================
    # Plot results
    #===========================================================================
    #if htmlplot:  pyclaw.plot.html_plot(outdir=outdir)
    #if iplot:     pyclaw.plot.interactive_plot(outdir=outdir)

if __name__=="__main__":

    from pyclaw.util import run_app_from_main
    output = run_app_from_main(advection_implicitLW)
    print 'Error: ',output

    state.problem_data['cc']=sqrt(bulk/rho)
 

    #========================================================================
    # Set the initial condition
    #========================================================================
    xc=domain.grid.x.centers
    beta=100; gamma=0; x0=0.75
    state.q[0,:] = exp(-beta * (xc-x0)**2) * cos(gamma * (xc - x0))
    state.q[1,:] = 0.
    
    #========================================================================
    # Set up the controller object
    #========================================================================
    claw = pyclaw.Controller()
    claw.solution = pyclaw.Solution(state,domain)
    claw.solver = solver
    claw.outdir = outdir
    claw.tfinal = 1.0

    # Solve
    status = claw.run()

    # Plot results
    if htmlplot:  pyclaw.plot.html_plot(outdir=outdir)
    if iplot:     pyclaw.plot.interactive_plot(outdir=outdir)

if __name__=="__main__":
    from pyclaw.util import run_app_from_main
    output = run_app_from_main(acoustics)

    hv = np.zeros((mx,my))
    hw = np.zeros((mx,my))

    # Read solution
    for j in range(my):
        tmp = np.fromfile(f,dtype='float',sep=" ",count=4*mx)
        tmp = tmp.reshape((mx,4))
        h[:,j] = tmp[:,0]
        hu[:,j] = tmp[:,1]
        hv[:,j] = tmp[:,2]
        hw[:,j] = tmp[:,3]

    
    # Plot solution in the computational domain
    # =========================================

    # Fluid height
    plt.figure()
    CS = plt.contour(xc,yc,h)
    plt.title('Fluid height (computational domain)')
    plt.xlabel('xc')
    plt.ylabel('yc')
    plt.clabel(CS, inline=1, fontsize=10)
    plt.show()


if __name__=="__main__":
    from pyclaw.util import run_app_from_main
    output = run_app_from_main(contourLineSphere)
    print 'Error: ',output