Index: trunk/anuga_work/development/2010-projects/anuga_1d/sqpipe/culvert.py
===================================================================
--- trunk/anuga_work/development/2010-projects/anuga_1d/sqpipe/culvert.py	(revision 8862)
+++ trunk/anuga_work/development/2010-projects/anuga_1d/sqpipe/culvert.py	(revision 8862)
@@ -0,0 +1,200 @@
+import os
+from math import sqrt, pi, sin, cos
+import numpy
+import time
+
+from anuga_1d.config import g, epsilon
+from anuga_1d.base.generic_mesh import uniform_mesh
+import anuga_1d.sqpipe.sqpipe_domain as dom
+
+#===============================================================================
+# setup problem
+#===============================================================================
+
+L_x = 50.0         ## length of channel
+
+def elevation(x):
+    z = -x*0.1
+    return z
+
+def stage(x):
+
+    w = elevation(x) + 1.0
+
+    return w
+
+
+def height(x):
+
+    h = stage(x) - elevation(x)
+
+    return h
+
+def width(x):
+    return numpy.ones_like(x)
+
+def top(x):
+
+    t = numpy.ones_like(x)
+
+    t = numpy.where(x<-40, 20, t)
+    t = numpy.where(x>40, 20, t)
+    return t
+
+def area(x):
+    return height(x)*width(x)
+
+
+#===============================================================================
+
+def get_domain():
+    N = 100
+    print "Evaluating domain with %d cells" %N
+
+    points, boundary = uniform_mesh(N, x_0 = -L_x, x_1 = L_x)
+
+    domain = dom.Domain(points, boundary, bulk_modulus = 500.0)
+
+    domain.set_spatial_order(2)
+    domain.set_timestepping_method('euler')
+    domain.set_CFL(0.5)
+    domain.set_limiter("vanleer")
+
+    domain.set_beta(1.0)
+
+    domain.set_quantity('area', area)
+    domain.set_quantity('stage', stage)
+    domain.set_quantity('elevation',elevation)
+    domain.set_quantity('width',width)
+    domain.set_quantity('top',top)
+
+    Br = dom.Reflective_boundary(domain)
+
+    domain.set_boundary({'left': Br, 'right' : Br})
+
+    return domain
+
+def animate_domain(domain, yieldstep, finaltime):
+    import pylab
+
+    pylab.ion()
+
+    x, z, w, h, v, t, s, m, M = get_quantities(domain)
+
+    zplot, wplot, ztplot, hplot, tplot, vplot, splot, Mplot, mplot = make_plots(x, z, w, h, v, t, s, m, M)
+
+    for t in domain.evolve(yieldstep = yieldstep, finaltime = finaltime):
+
+        domain.write_time()
+        x, z, w, h, v, t, s, m, M = get_quantities(domain)
+        
+        zplot.set_ydata(z)
+        ztplot.set_ydata(z+t)
+        wplot.set_ydata(w)
+        hplot.set_ydata(h)
+        tplot.set_ydata(t)
+        vplot.set_ydata(v)
+        splot.set_ydata(s)
+        mplot.set_ydata(m)
+        Mplot.set_ydata(M)
+
+
+        pylab.ion()
+        pylab.draw()
+
+
+
+
+def plot_domain(domain):
+    import pylab
+
+    x, z, w, h, v, t, s, m, M = get_quantities(domain)
+
+    pylab.ioff()
+    pylab.hold(False)
+
+    make_plots(x, z, w, h, v, t, s, m, M)
+    
+    pylab.show()
+
+def write_domain(domain, outfile):
+    x = domain.get_centroids()
+    z = domain.get_quantity('elevation', 'centroids')
+    w = domain.get_quantity('stage', 'centroids')
+
+    f = open(outfile, 'w')
+    for i in range(len(x)):
+        f.write("%s %s %s\n" % (x[i], z[i], w[i]))
+    f.close()
+
+def get_quantities(domain):
+    x = domain.get_centroids()
+    z = domain.get_quantity('elevation', 'centroids')
+    w = domain.get_quantity('stage', 'centroids')
+    h = domain.get_quantity('height', 'centroids')
+    v = domain.get_quantity('velocity', 'centroids')
+    t = domain.get_quantity('top', 'centroids')
+    s = domain.state
+    m = domain.get_mass()
+    M = m.sum() * numpy.ones_like(x)
+
+    return x, z, w, h, v, t, s, m, M
+
+def make_plots(x, z, w, h, v, t, s, m, M):
+    import pylab
+
+    #fig = pylab.gcf()
+    #fig.set_size_inches(12,12, forward=True)
+
+    
+
+
+    plot1 = pylab.subplot(321)
+    zplot, = pylab.plot(x, z)
+    wplot, = pylab.plot(x, w)
+    ztplot, = pylab.plot(x, z+t)
+    plot1.set_xlim([-60,60])
+    plot1.set_ylim([-10,10])
+    pylab.xlabel('Position')
+    pylab.ylabel('Stage')
+
+    plot2 = pylab.subplot(322)
+    hplot, = pylab.plot(x, h)
+    tplot, = pylab.plot(x, t)
+    plot2.set_xlim([-60,60])
+    plot2.set_ylim([-1,5])
+    pylab.xlabel('Position')
+    pylab.ylabel('Height')
+
+    plot3 = pylab.subplot(323)
+    vplot, = pylab.plot(x, v)
+    plot3.set_xlim([-60,60])
+    plot3.set_ylim([-6,6])
+    pylab.xlabel('Position')
+    pylab.ylabel('Velocity')
+
+    plot4 = pylab.subplot(324)
+    splot, = pylab.plot(x, s)
+    plot4.set_xlim([-60,60])
+    plot4.set_ylim([-1,2])
+    pylab.xlabel('Position')
+    pylab.ylabel('State')
+
+    plot5 = pylab.subplot(325)
+    mplot, = pylab.plot(x, m)
+    plot5.set_xlim([-60,60])
+    plot5.set_ylim([-1,10])
+    pylab.xlabel('Position')
+    pylab.ylabel('Mass')
+
+    plot6 = pylab.subplot(326)
+    Mplot, = pylab.plot(x, M)
+    plot6.set_xlim([-60,60])
+    plot6.set_ylim([-1,450])
+    pylab.xlabel('Position')
+    pylab.ylabel('Total Mass')
+
+
+
+
+    return zplot, wplot, ztplot, hplot, tplot, vplot, splot, Mplot, mplot
Index: trunk/anuga_work/development/2010-projects/anuga_1d/sqpipe/parabolic_canal.py
===================================================================
--- trunk/anuga_work/development/2010-projects/anuga_1d/sqpipe/parabolic_canal.py	(revision 8861)
+++ trunk/anuga_work/development/2010-projects/anuga_1d/sqpipe/parabolic_canal.py	(revision 8862)
@@ -6,4 +6,5 @@
 from anuga_1d.config import g, epsilon
 from anuga_1d.base.generic_mesh import uniform_mesh
+import anuga_1d.sqpipe.sqpipe_domain as dom
 
 #===============================================================================
@@ -66,7 +67,8 @@
 #===============================================================================
 
-def get_domain(dom):
+def get_domain():
     N = 100
     print "Evaluating domain with %d cells" %N
+    
     domain = dom.Domain(*uniform_mesh(N, x_0 = -2.0*L_x, x_1 = 2.0*L_x), bulk_modulus = 100.0)
 
@@ -99,4 +101,6 @@
 
     for t in domain.evolve(yieldstep = yieldstep, finaltime = finaltime):
+
+        domain.write_time()
         x, z, w, h, v, t, s, m, M = get_quantities(domain)
         
Index: trunk/anuga_work/development/2010-projects/anuga_1d/sqpipe/run_culvert.py
===================================================================
--- trunk/anuga_work/development/2010-projects/anuga_1d/sqpipe/run_culvert.py	(revision 8862)
+++ trunk/anuga_work/development/2010-projects/anuga_1d/sqpipe/run_culvert.py	(revision 8862)
@@ -0,0 +1,14 @@
+
+import time
+
+
+import culvert as model
+
+domain = model.get_domain()
+
+finaltime = 100.0
+yieldstep = 0.1
+
+model.animate_domain(domain, yieldstep, finaltime)
+print "finished"
+model.plot_domain(domain)
Index: trunk/anuga_work/development/2010-projects/anuga_1d/sqpipe/run_parabolic_canal.py
===================================================================
--- trunk/anuga_work/development/2010-projects/anuga_1d/sqpipe/run_parabolic_canal.py	(revision 8862)
+++ trunk/anuga_work/development/2010-projects/anuga_1d/sqpipe/run_parabolic_canal.py	(revision 8862)
@@ -0,0 +1,14 @@
+
+import time
+
+
+import parabolic_canal as model
+
+domain = model.get_domain()
+finaltime = 1000.0
+yieldstep = 10.0
+
+
+model.animate_domain(domain, yieldstep, finaltime)
+print "finished"
+model.plot_domain(domain)
Index: trunk/anuga_work/development/2010-projects/anuga_1d/sqpipe/sqpipe.h
===================================================================
--- trunk/anuga_work/development/2010-projects/anuga_1d/sqpipe/sqpipe.h	(revision 8861)
+++ trunk/anuga_work/development/2010-projects/anuga_1d/sqpipe/sqpipe.h	(revision 8862)
@@ -133,6 +133,8 @@
   struct sqpipe_params *sp = (struct sqpipe_params *) p;
 
-  quantityflux[0] = (normal * pq->u) * quantity_get_conserved(q, 0, normal);
-  quantityflux[1] = normal * (pq->u * quantity_get_conserved(q, 1, normal) + 0.5 * sp->g * pq->h *pq->h);
+  //quantityflux[0] = (normal * pq->u) * quantity_get_conserved(q, 0, normal);
+  //quantityflux[1] = normal * (pq->u * quantity_get_conserved(q, 1, normal) + 0.5 * sp->g * pq->h *pq->h);
+  quantityflux[0] = (normal * pq->u) * pq->a;
+  quantityflux[1] = normal * (pq->u * pq->d + 0.5 * sp->g * pq->h *pq->h);
 
   return quantityflux;
Index: trunk/anuga_work/development/2010-projects/anuga_1d/sqpipe/sqpipe_domain.py
===================================================================
--- trunk/anuga_work/development/2010-projects/anuga_1d/sqpipe/sqpipe_domain.py	(revision 8861)
+++ trunk/anuga_work/development/2010-projects/anuga_1d/sqpipe/sqpipe_domain.py	(revision 8862)
@@ -152,5 +152,6 @@
         # Import flux method and call it
         from anuga_1d.sqpipe.sqpipe_comp_flux import compute_fluxes_ext
- 	self.flux_timestep = compute_fluxes_ext(timestep,self)        
+
+        self.flux_timestep = compute_fluxes_ext(timestep,self)
 
     def update_state(self):
Index: trunk/anuga_work/development/2010-projects/anuga_1d/sqpipe/test_sqpipe_domain.py
===================================================================
--- trunk/anuga_work/development/2010-projects/anuga_1d/sqpipe/test_sqpipe_domain.py	(revision 8861)
+++ 	(revision )
@@ -1,22 +1,0 @@
-import anuga_1d.sqpipe.parabolic_canal
-import time
-import anuga_1d.sqpipe.sqpipe_domain as dom
-
-
-domain = anuga_1d.sqpipe.parabolic_canal.get_domain(dom)
-finaltime = 1000.0
-yieldstep = 10.0
-#m0 = (domain.get_mass()).sum()
-#for t in domain.evolve(yieldstep = yieldstep, finaltime = finaltime):
-#	m = (domain.get_mass()).sum()
-#	if abs(m - m0) > 0.5:
-#		print m0 - m
-#		
-#	m0 = m
-	
-#outfile="out"
-#anuga_1d.sqpipe.parabolic_canal.write_domain(domain, outfile)
-
-anuga_1d.sqpipe.parabolic_canal.animate_domain(domain, yieldstep, finaltime)
-print "finished"
-anuga_1d.sqpipe.parabolic_canal.plot_domain(domain)
