Index: development/pyvolution-1d/pt.py
===================================================================
--- development/pyvolution-1d/pt.py	(revision 3294)
+++ development/pyvolution-1d/pt.py	(revision 3294)
@@ -0,0 +1,194 @@
+from math import sqrt, pi
+from shallow_water_1d import *
+from Numeric import allclose, array, zeros, ones, Float, take
+from config import g, epsilon
+
+def test_sqrt():
+    for i in range (80000):
+        a = sqrt(4356)
+        b = sqrt(2031)
+
+def test_flux1():
+    #Use data from previous version of pyvolution
+    #normal = array([1.,0])
+    #ql = array([-0.2, 2, 3])
+    #qr = array([-0.2, 2, 3])
+    ql = array([-0.2, 2])
+    qr = array([-0.2, 2])
+    zl = zr = -0.5
+    #flux, max_speed = flux_function(normal, ql, qr, zl, zr)
+    flux, max_speed = flux_function(1.0, ql, qr, zl, zr)
+    print 'flux', flux
+
+def test_compute_fluxes0():
+    #Do a full triangle and check that fluxes cancel out for
+    #the constant stage case
+    
+    print 'check min time step in compute fluxes is ok, John'
+    
+    #a = [0.0, 0.0]
+    #b = [0.0, 2.0]
+    #c = [2.0,0.0]
+    #d = [0.0, 4.0]
+    #e = [2.0, 2.0]
+    #f = [4.0,0.0]
+    a=0.0
+    b=2.0
+    c=4.0
+    d=6.0
+    e=8.0
+    f=10.0
+    
+    points = [a, b, c, d, e, f]
+    #bac, bce, ecf, dbe
+    #vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
+    
+    domain = Domain(points)
+    #domain.set_quantity('stage', [[2,2,2], [2,2,2],
+    #                              [2,2,2], [2,2,2]])
+    domain.set_quantity('stage', [[2,2], [2,2], [2,2], [2,2], [2,2]])
+    domain.check_integrity()
+    
+    #assert allclose(domain.neighbours, [[-1,1,-1], [2,3,0], [-1,-1,1],[1,-1,-1]])
+    assert allclose(domain.neighbours, [[-1,1], [0,2], [1,3],[2,4], [3,-1]])
+    #assert allclose(domain.neighbour_edges, [[-1,2,-1], [2,0,1], [-1,-1,0],[1,-1,-1]])
+    #assert allclose(domain.neighbour_edges, [[-1,0], [1,0], [1,0], [1,0], [1,-1]])
+    assert allclose(domain.neighbour_vertices, [[-1,0], [1,0], [1,0], [1,0], [1,-1]])
+    
+    zl=zr=0. #Assume flat bed
+    
+    #Flux across right edge of volume 1
+    #normal = domain.get_normal(1,0)
+    #ql = domain.get_conserved_quantities(vol_id=1, edge=0)
+    #qr = domain.get_conserved_quantities(vol_id=2, edge=2)
+    #flux0, max_speed = flux_function(normal, ql, qr, zl, zr)
+    
+    #Flux across right edge of element 1
+    ql = domain.get_conserved_quantities(vol_id=1, vertex=1)
+    qr = domain.get_conserved_quantities(vol_id=2, vertex=0)
+    #print 'qr and ql 1'
+    #print qr
+    #print ql
+    flux0, max_speed = flux_function(1.0, ql, qr, zl, zr)
+    
+    #Check that flux seen from other triangles is inverse
+    tmp = qr; qr=ql; ql=tmp
+    #print 'qr and ql 2'
+    #print qr
+    #print ql
+    #normal = domain.get_normal(2,2)
+    #flux, max_speed = flux_function(normal, ql, qr, zl, zr)
+    flux, max_speed = flux_function(-1.0, ql, qr, zl, zr)
+    #print 'fluxes'
+    #print flux0
+    #print flux
+    assert allclose(flux + flux0, 0.)
+    
+    #Flux across upper edge of volume 1
+    #normal = domain.get_normal(1,1)
+    #ql = domain.get_conserved_quantities(vol_id=1, edge=1)
+    #qr = domain.get_conserved_quantities(vol_id=3, edge=0)
+    #flux1, max_speed = flux_function(normal, ql, qr, zl, zr)
+    
+    #Flux across left edge of element 1
+    #ql = domain.get_conserved_quantities(vol_id=1, edge=0)
+    #qr = domain.get_conserved_quantities(vol_id=0, edge=1)
+    ql = domain.get_conserved_quantities(vol_id=1, vertex=0)
+    qr = domain.get_conserved_quantities(vol_id=0, vertex=1)
+    flux1, max_speed = flux_function(-1.0, ql, qr, zl, zr)
+    
+    #Check that flux seen from other triangles is inverse
+    tmp = qr; qr=ql; ql=tmp
+    #normal = domain.get_normal(3,0)
+    #flux, max_speed = flux_function(normal, ql, qr, zl, zr)
+    flux, max_speed = flux_function(1.0, ql, qr, zl, zr)
+    assert allclose(flux + flux1, 0.)
+    
+    #Flux across lower left hypotenuse of volume 1
+    #normal = domain.get_normal(1,2)
+    #ql = domain.get_conserved_quantities(vol_id=1, edge=2)
+    #qr = domain.get_conserved_quantities(vol_id=0, edge=1)
+    #flux2, max_speed = flux_function(normal, ql, qr, zl, zr)
+    
+    #Check that flux seen from other triangles is inverse
+    #tmp = qr; qr=ql; ql=tmp
+    #normal = domain.get_normal(0,1)
+    #flux, max_speed = flux_function(normal, ql, qr, zl, zr)
+    #assert allclose(flux + flux2, 0.)
+    
+    
+    #Scale by edgelengths, add up anc check that total flux is zero
+    #e0 = domain.edgelengths[1, 0]
+    #e1 = domain.edgelengths[1, 1]
+    #e2 = domain.edgelengths[1, 2]
+    
+    #assert allclose(e0*flux0+e1*flux1+e2*flux2, 0.)
+    print 'flux0',flux0
+    print 'flux1',flux1
+    assert allclose(flux0+flux1, 0.)
+    
+    #Now check that compute_flux yields zeros as well
+    domain.compute_fluxes()
+    
+    #for name in ['stage', 'xmomentum', 'ymomentum']:
+    for name in ['stage', 'xmomentum']:
+        #print name, domain.quantities[name].explicit_update
+        assert allclose(domain.quantities[name].explicit_update[1], 0)
+        
+def test_1d_solution_I():
+    print "TEST 1D-SOLUTION I"
+    
+    L = 2000.0     # Length of channel (m)
+    N = 100        # Number of compuational cells
+    cell_len = L/N # Origin = 0.0
+    
+    points = zeros(N+1,Float)
+    for i in range(N+1):
+        points[i] = i*cell_len
+        
+    domain = Domain(points)
+    domain.order = 2
+    def stage(x):
+        for i in range(len(x)):
+            if x[i]<=1000.0:
+                x[i] = 10.0
+            else:
+                x[i] = 5.0
+        return x
+
+    domain.set_quantity('stage', stage)
+    #L = domain.quantities['stage'].vertex_values
+    #print "Initial Stage"
+    #print L
+
+    domain.set_boundary({'exterior': Reflective_boundary(domain)})
+    
+    import time
+    t0 = time.time()
+    yieldstep = 1.0
+    finaltime = 50.0
+    for t in domain.evolve(yieldstep = yieldstep, finaltime = finaltime):
+        #xmom = domain.quantities['xmomentum']
+        #xmom = xmom.centroid_values
+        #stage = domain.quantities['stage']
+        #stage = stage.centroid_values
+        #print 'stage', stage
+        #print 'xmom', xmom
+        #for i in range N
+        #u = xmom/stage
+        pass
+    
+    print 'That took %.2f seconds' %(time.time()-t0)
+    
+    #L = domain.quantities['stage'].vertex_values
+    #print "Final Stage"
+    #print L
+    
+    C = domain.quantities['stage'].vertex_values
+    #print C
+    
+    f = file('test_solution_Ix.out', 'w')
+    for i in range(N):
+        f.write(str(C[i,1]))
+        f.write("\n")
+    f.close
