Index: uga_validation/analytical solutions/Analytical solution_oblique_shock_01.py
===================================================================
--- /anuga_validation/analytical solutions/Analytical solution_oblique_shock_01.py	(revision 7623)
+++ 	(revision )
@@ -1,194 +1,0 @@
-"""
-Example of shallow water wave equation
-consisting of an asymetrical converging channel.
-
-Copyright 2004
-Christopher Zoppou, Stephen Roberts, Ole Nielsen, Duncan Gray
-Geoscience Australia, ANU
-
-Specific methods pertaining to the 2D shallow water equation
-are imported from shallow_water
-for use with the generic finite volume framework
-
-Conserved quantities are h, uh and vh stored as elements 0, 1 and 2 in the
-numerical vector named conserved_quantities.
-
-"""
-
-######################
-# Module imports
-#
-
-#Were these used?
-#import visualise2_chris as visualise
-#import Image, ImageGrab
-
-from anuga.shallow_water import Domain
-from anuga.shallow_water import Transmissive_boundary, Reflective_boundary,\
-     Dirichlet_boundary , Transmissive_Momentum_Set_Stage_boundary
-
-from math import sqrt, cos, sin, pi
-from mesh_factory import oblique_cross
-
-
-######################
-# Domain
-#
-leny = 30.
-lenx = 40.
-f = 2
-n = 50*f
-m = 80*f
-theta = 25
-h_bc = 1.5
-p_bc = 5.0
-
-points, elements, boundary = oblique_cross(m, n, lenx, leny, theta = theta)
-domain = Domain(points, elements, boundary)
-
-
-print 'Number of Elements = ',domain.number_of_elements
-# Order of solver
-domain.default_order=2
-domain.beta_w = 0.8
-
-# Store output
-domain.store=True
-
-# Provide file name for storing output
-domain.set_name('oblique_cross_%g_%g_%g_%g_%g' % (n,m,theta,h_bc,p_bc))
-print domain.get_name()
-
-# Output format
-domain.format="sww" #NET.CDF binary format
-                    # "dat" for ASCII
-
-
-# Visualization smoothing
-domain.smooth=True
-domain.visualise=False
-
-#######################
-#Bed-slope and friction
-class ConstantFunctionT:
-
-    def __init__(self,value):
-        self.value = value
-
-    def __call__(self,t):
-        return self.value
-
-
-shock_hb = ConstantFunctionT(h_bc)
-shock_hh = ConstantFunctionT(h_bc)
-shock_pb = ConstantFunctionT(p_bc)
-shock_pp = ConstantFunctionT(p_bc)
-
-
-def x_slope(x, y):
-    return 0*x
-
-
-
-def shock_h(x,y):
-    n = x.shape[0]
-    w = 0*x
-    for i in range(n):
-        if x[i]<0:
-            w[i] = shock_hh(0.0)
-        else:
-            w[i] = 1.0
-    return w
-
-
-def shock_p(x,y):
-    n = x.shape[0]
-    w = 0*x
-    for i in range(n):
-        if x[i]<0:
-            w[i] = shock_pp(0.0)
-        else:
-            w[i] = 0.0
-    return w
-
-
-
-
-
-domain.set_quantity('elevation', x_slope)
-domain.set_quantity('friction', 0.0)
-
-######################
-# Boundary conditions
-#
-R = Reflective_boundary(domain)
-T = Transmissive_boundary(domain)
-D = Dirichlet_boundary([shock_hb(0.0) , shock_pb(0.0), 0.0])
-Bts = Transmissive_Momentum_Set_Stage_boundary(domain, shock_hb)
-
-domain.set_boundary({'left': D, 'right': T, 'top': R, 'bottom': R})
-
-######################
-#Initial condition
-
-domain.set_quantity('stage', shock_h )
-domain.set_quantity('xmomentum',shock_p )
-
-class SetValueWhere:
-
-    def __init__(self,quantity,value=0,xrange=0.0):
-        self.quantity =  quantity
-        self.value = value
-        self.xrange = xrange
-
-    def __call__(self,x,y):
-        n = x.shape[0]
-        w = self.quantity.centroid_values
-        for i in range(n):
-            if x[i]<self.xrange:
-                w[i] = self.value
-        return w
-
-
-Stage = domain.quantities['stage']
-Xmom =  domain.quantities['xmomentum']
-Ymom =  domain.quantities['ymomentum']
-
-#for id, face in domain.boundary:
-#    print '(',id,',',face,') ',domain.boundary[(id,face)]
-
-######################
-#Evolution
-import time
-t0 = time.time()
-for t in domain.evolve(yieldstep = 0.1, finaltime = 20.0):
-    domain.write_time()
-    id = 3399
-    print Stage.get_values(location='centroids',indices=[id])
-    print Xmom.get_values(location='centroids',indices=[id])
-    print Ymom.get_values(location='centroids',indices=[id])
-    vstage =  Stage.get_values(location='centroids',indices=[id])
-    vxmom  =  Xmom.get_values(location='centroids',indices=[id])
-    id = 12719
-    print Stage.get_values(location='centroids',indices=[id])
-    print Xmom.get_values(location='centroids',indices=[id])
-    print Ymom.get_values(location='centroids',indices=[id])
-    tclean = 2.0
-#    if t > tclean-0.09 and t < tclean+0.01:
-#        print 'Cleaning up Initial profile'
-#        setstage = SetValueWhere(quantity=Stage,value=vstage[0],xrange=10)
-#        setxmom  = SetValueWhere(quantity=Xmom,value=vxmom[0],xrange=10)
-#        Stage.set_values(setstage,location='centroids')
-#        Xmom.set_values(setxmom,location='centroids')
-#        Dnew = Dirichlet_boundary([vstage[0] , vxmom[0], 0.0])
-#        domain.set_boundary({'left': Dnew, 'right': T, 'top': R, 'bottom': R})
-
-print 'That took %.2f seconds' %(time.time()-t0)
-
-#FIXME: Compute average water depth on either side of shock and compare
-#to expected values. And also Froude numbers.
-
-
-#print "saving file?"
-#im = ImageGrab.grab()
-#im.save("ccube.eps")
Index: uga_validation/analytical solutions/Analytical solution_wave_runup.py
===================================================================
--- /anuga_validation/analytical solutions/Analytical solution_wave_runup.py	(revision 7623)
+++ 	(revision )
@@ -1,177 +1,0 @@
-"""Example of shallow water wave equation analytical solution of the
-one-dimensional Thacker and Greenspan wave run-up treated as a two-dimensional solution.
-
-   Copyright 2004
-   Christopher Zoppou, Stephen Roberts, Ole Nielsen, Duncan Gray
-   Geoscience Australia
-
-Specific methods pertaining to the 2D shallow water equation
-are imported from shallow_water
-for use with the generic finite volume framework
-
-Conserved quantities are h, uh and vh stored as elements 0, 1 and 2 in the
-numerical vector named conserved_quantities.
-"""
-
-######################
-# Module imports
-import sys
-from os import sep
-sys.path.append('..'+sep+'pyvolution')
-
-from shallow_water import Domain, Reflective_boundary, Time_boundary
-from math import sqrt, cos, sin, pi
-from mesh_factory import strang_mesh
-
-#Convenience functions
-def imag(a):
-    return a.imag
-
-def real(a):
-    return a.real
-
-
-######################
-# Domain
-# Strang_domain will search through the file and test to see if there are
-# two or three entries. Two entries are for points and three for triangles.
-
-#points, elements = strang_mesh('Run-up.pt')
-points, elements = strang_mesh('strang_7389.pt')
-domain = Domain(points, elements)
-
-domain.default_order = 2
-domain.smooth = True
-
-#Set a default tagging
-epsilon = 1.0e-12
-
-for id, face in domain.boundary:
-    domain.boundary[(id,face)] = 'external'
-    if domain.get_vertex_coordinate(id,(face+1)%3)[0] < -200.0+1.0e-10:
-        domain.boundary[(id,face)] = 'left'
-
-
-
-
-# Provide file name for storing output
-domain.store = True
-domain.format = 'sww'
-domain.set_name('run-up_second_order_again')
-
-print "Number of triangles = ", len(domain)
-
-#Reduction operation for get_vertex_values
-from anuga.pyvolution.util import mean
-domain.reduction = mean
-#domain.reduction = min  #Looks better near steep slopes
-
-#Define the boundary condition
-def stage_setup(x,t_star1):
-    vh = 0
-    alpha = 0.1
-    eta = 0.1
-    a = 1.5*sqrt(1.+0.9*eta)
-    l_0 = 200.
-    ii = complex(0,1)
-    g = 9.81
-    v_0 = sqrt(g*l_0*alpha)
-    v1 = 0.
-
-    sigma_max = 100.
-    sigma_min = -100.
-    for j in range (1,50):
-        sigma0 = (sigma_max+sigma_min)/2.
-        lambda_prime = 2./a*(t_star1/sqrt(l_0/alpha/g)+v1)
-        sigma_prime = sigma0/a
-        const = (1.-ii*lambda_prime)**2+sigma_prime**2
-
-        v1 = 8.*eta/a*imag(1./const**(3./2.) \
-            -3./4.*(1.-ii*lambda_prime)/const**(5./2.))
-
-        x1 = -v1**2/2.-a**2*sigma_prime**2/16.+eta \
-            *real(1.-2.*(5./4.-ii*lambda_prime) \
-            /const**(3./2.)+3./2.*(1.-ii*lambda_prime)**2 \
-            /const**(5./2.))
-
-        neta1 = x1 + a*a*sigma_prime**2/16.
-
-        v_star1 = v1*v_0
-        x_star1 = x1*l_0
-        neta_star1 = neta1*alpha*l_0
-        stage = neta_star1
-        z = stage - x_star1*alpha
-        uh = z*v_star1
-
-        if x_star1-x > 0:
-            sigma_max = sigma0
-        else:
-            sigma_min = sigma0
-
-        if abs(abs(sigma0)-100.) < 10:
-
-#     solution does not converge because bed is dry
-            stage = 0.
-            uh = 0.
-            z = 0.
-
-    return [stage, uh, vh]
-
-def boundary_stage(t):
-    x = -200
-    return stage_setup(x,t)
-
-
-######################
-#Initial condition
-print 'Initial condition'
-t_star1 = 0.0
-slope = -0.1
-
-#Set bed-elevation and friction(None)
-def x_slope(x,y):
-    n = x.shape[0]
-    z = 0*x
-    for i in range(n):
-        z[i] = -slope*x[i]
-    return z
-
-domain.set_quantity('elevation', x_slope)
-
-#Set the water depth
-print 'Initial water depth'
-def stage(x,y):
-    z = x_slope(x,y)
-    n = x.shape[0]
-    w = 0*x
-    for i in range(n):
-        w[i], uh, vh = stage_setup(x[i],t_star1)
-        h = w[i] - z[i]
-        if h < 0:
-            h = 0
-            w[i] = z[i]
-    return w
-
-domain.set_quantity('stage', stage)
-
-
-#####################
-#Set up boundary conditions
-Br = Reflective_boundary(domain)
-Bw = Time_boundary(domain, boundary_stage)
-domain.set_boundary({'left': Bw, 'external': Br})
-
-
-
-#domain.visualise = True
-
-######################
-#Evolution
-import time
-t0 = time.time()
-for t in domain.evolve(yieldstep = 1., finaltime = 100):
-    domain.write_time()
-    print boundary_stage(domain.time)
-
-print 'That took %.2f seconds' %(time.time()-t0)
-
Index: /anuga_validation/analytical solutions/Analytical_solution_Yoon_cross_mesh.py
===================================================================
--- /anuga_validation/analytical solutions/Analytical_solution_Yoon_cross_mesh.py	(revision 7623)
+++ /anuga_validation/analytical solutions/Analytical_solution_Yoon_cross_mesh.py	(revision 7624)
@@ -1,5 +1,5 @@
 """Example of shallow water wave equation analytical solution
 consists of a parabolic profile in a parabolic basin. Analytical
-solutiuon to this problem was derived by Carrier and Greenspan
+solutiuon to this problem was derived by Thacker
 and used by Yoon and Chou.
 
@@ -14,20 +14,21 @@
 #from os import sep
 #sys.path.append('..'+sep+'pyvolution')
-from anuga.pyvolution.shallow_water import Domain, Transmissive_boundary, Reflective_boundary,\
+from anuga.shallow_water_balanced.swb_domain import Domain, Transmissive_boundary, Reflective_boundary,\
      Dirichlet_boundary
+
+#from anuga.interface import Domain, Transmissive_boundary, Reflective_boundary,\
+#     Dirichlet_boundary
 from math import sqrt, cos, sin, pi
-from anuga.pyvolution.mesh_factory import rectangular_cross
-from anuga.utilities.polygon import inside_polygon
-from Numeric import asarray
-from anuga.pyvolution.least_squares import Interpolation
+from anuga.interface import rectangular_cross
+from anuga.utilities.polygon import inside_polygon, is_inside_triangle
+from numpy import asarray
+
 
 #-------------------------------
 # Domain
-n = 100
-m = 100
-delta_x = 80.0
-delta_y = 80.0
-lenx = delta_x*n
-leny = delta_y*m
+n = 200
+m = 200
+lenx = 8000.0
+leny = 8000.0
 origin = (-4000.0, -4000.0)
 
@@ -37,5 +38,11 @@
 #----------------
 # Order of scheme
-domain.default_order = 1
+# Good compromise between
+# limiting and CFL
+#---------------
+domain.set_default_order(2)
+domain.set_timestepping_method(2)
+domain.set_beta(0.7)
+domain.set_CFL(0.6)
 
 domain.smooth = True
@@ -54,5 +61,5 @@
 #------------------------------------------
 # Reduction operation for get_vertex_values
-from anuga.pyvolution.util import mean
+from anuga.utilities.numerical_tools import mean
 domain.reduction = mean #domain.reduction = min  #Looks better near steep slopes
 
@@ -107,16 +114,19 @@
 # Find triangle that contains the point points
 # and print to file
-points = [0.,0.]
+
+
+
+points = (0.0, 0.0)
 for n in range(len(domain.triangles)):
-    t = domain.triangles[n]
-    tri = [domain.coordinates[t[0]],domain.coordinates[t[1]],domain.coordinates[t[2]]]
+    tri = domain.get_vertex_coordinates(n)
 
-    if inside_polygon(points,tri):
-        print 'Point is within triangle with vertices '+'%s'%tri
+    if is_inside_triangle(points,tri):
+        #print 'Point is within triangle with vertices '+'%s'%tri
         n_point = n
 
 print 'n_point = ',n_point
 t = domain.triangles[n_point]
-tri = [domain.coordinates[t[0]],domain.coordinates[t[1]],domain.coordinates[t[2]]]
+print 't = ', t
+tri = domain.get_vertex_coordinates(n)
 
 filename=domain.get_name()
@@ -137,12 +147,13 @@
     domain.write_time()
 
-    tri_array = asarray(tri)
-    t_array = asarray([[0,1,2]])
-    interp = Interpolation(tri_array,t_array,[points])
+    #tri_array = asarray(tri)
+    #t_array = asarray([[0,1,2]])
+    #interp = Interpolation(tri_array,t_array,[points])
 
 
-    stage     = Stage.get_values(location='centroids',indices=[n_point])[0]
-    xmomentum = Xmomentum.get_values(location='centroids',indices=[n_point])[0]
-    ymomentum = Ymomentum.get_values(location='centroids',indices=[n_point])[0]
+    stage     = Stage.get_values(location='centroids',indices=[n_point])
+    xmomentum = Xmomentum.get_values(location='centroids',indices=[n_point])
+    ymomentum = Ymomentum.get_values(location='centroids',indices=[n_point])
+    #print '%10.6f   %10.6f  %10.6f   %10.6f\n'%(t,stage,xmomentum,ymomentum)
     file.write( '%10.6f   %10.6f  %10.6f   %10.6f\n'%(t,stage,xmomentum,ymomentum) )
 
@@ -161,5 +172,5 @@
 xlabel('time(s)')
 ylabel('stage (m)')
-legend(('Observed', 'Modelled'), shadow=True, loc='upper left')
+#legend(('Observed', 'Modelled'), shadow=True, loc='upper left')
 #savefig(name, dpi = 300)
 
Index: /anuga_validation/analytical solutions/Analytical_solution_oblique_shock_01.py
===================================================================
--- /anuga_validation/analytical solutions/Analytical_solution_oblique_shock_01.py	(revision 7624)
+++ /anuga_validation/analytical solutions/Analytical_solution_oblique_shock_01.py	(revision 7624)
@@ -0,0 +1,194 @@
+"""
+Example of shallow water wave equation
+consisting of an asymetrical converging channel.
+
+Copyright 2004
+Christopher Zoppou, Stephen Roberts, Ole Nielsen, Duncan Gray
+Geoscience Australia, ANU
+
+Specific methods pertaining to the 2D shallow water equation
+are imported from shallow_water
+for use with the generic finite volume framework
+
+Conserved quantities are h, uh and vh stored as elements 0, 1 and 2 in the
+numerical vector named conserved_quantities.
+
+"""
+
+######################
+# Module imports
+#
+
+#Were these used?
+#import visualise2_chris as visualise
+#import Image, ImageGrab
+
+from anuga.shallow_water import Domain
+from anuga.shallow_water import Transmissive_boundary, Reflective_boundary,\
+     Dirichlet_boundary , Transmissive_Momentum_Set_Stage_boundary
+
+from math import sqrt, cos, sin, pi
+from mesh_factory import oblique_cross
+
+
+######################
+# Domain
+#
+leny = 30.
+lenx = 40.
+f = 2
+n = 50*f
+m = 80*f
+theta = 25
+h_bc = 1.5
+p_bc = 5.0
+
+points, elements, boundary = oblique_cross(m, n, lenx, leny, theta = theta)
+domain = Domain(points, elements, boundary)
+
+
+print 'Number of Elements = ',domain.number_of_elements
+# Order of solver
+domain.default_order=2
+domain.beta_w = 0.8
+
+# Store output
+domain.store=True
+
+# Provide file name for storing output
+domain.set_name('oblique_cross_%g_%g_%g_%g_%g' % (n,m,theta,h_bc,p_bc))
+print domain.get_name()
+
+# Output format
+domain.format="sww" #NET.CDF binary format
+                    # "dat" for ASCII
+
+
+# Visualization smoothing
+domain.smooth=True
+domain.visualise=False
+
+#######################
+#Bed-slope and friction
+class ConstantFunctionT:
+
+    def __init__(self,value):
+        self.value = value
+
+    def __call__(self,t):
+        return self.value
+
+
+shock_hb = ConstantFunctionT(h_bc)
+shock_hh = ConstantFunctionT(h_bc)
+shock_pb = ConstantFunctionT(p_bc)
+shock_pp = ConstantFunctionT(p_bc)
+
+
+def x_slope(x, y):
+    return 0*x
+
+
+
+def shock_h(x,y):
+    n = x.shape[0]
+    w = 0*x
+    for i in range(n):
+        if x[i]<0:
+            w[i] = shock_hh(0.0)
+        else:
+            w[i] = 1.0
+    return w
+
+
+def shock_p(x,y):
+    n = x.shape[0]
+    w = 0*x
+    for i in range(n):
+        if x[i]<0:
+            w[i] = shock_pp(0.0)
+        else:
+            w[i] = 0.0
+    return w
+
+
+
+
+
+domain.set_quantity('elevation', x_slope)
+domain.set_quantity('friction', 0.0)
+
+######################
+# Boundary conditions
+#
+R = Reflective_boundary(domain)
+T = Transmissive_boundary(domain)
+D = Dirichlet_boundary([shock_hb(0.0) , shock_pb(0.0), 0.0])
+Bts = Transmissive_Momentum_Set_Stage_boundary(domain, shock_hb)
+
+domain.set_boundary({'left': D, 'right': T, 'top': R, 'bottom': R})
+
+######################
+#Initial condition
+
+domain.set_quantity('stage', shock_h )
+domain.set_quantity('xmomentum',shock_p )
+
+class SetValueWhere:
+
+    def __init__(self,quantity,value=0,xrange=0.0):
+        self.quantity =  quantity
+        self.value = value
+        self.xrange = xrange
+
+    def __call__(self,x,y):
+        n = x.shape[0]
+        w = self.quantity.centroid_values
+        for i in range(n):
+            if x[i]<self.xrange:
+                w[i] = self.value
+        return w
+
+
+Stage = domain.quantities['stage']
+Xmom =  domain.quantities['xmomentum']
+Ymom =  domain.quantities['ymomentum']
+
+#for id, face in domain.boundary:
+#    print '(',id,',',face,') ',domain.boundary[(id,face)]
+
+######################
+#Evolution
+import time
+t0 = time.time()
+for t in domain.evolve(yieldstep = 0.1, finaltime = 20.0):
+    domain.write_time()
+    id = 3399
+    print Stage.get_values(location='centroids',indices=[id])
+    print Xmom.get_values(location='centroids',indices=[id])
+    print Ymom.get_values(location='centroids',indices=[id])
+    vstage =  Stage.get_values(location='centroids',indices=[id])
+    vxmom  =  Xmom.get_values(location='centroids',indices=[id])
+    id = 12719
+    print Stage.get_values(location='centroids',indices=[id])
+    print Xmom.get_values(location='centroids',indices=[id])
+    print Ymom.get_values(location='centroids',indices=[id])
+    tclean = 2.0
+#    if t > tclean-0.09 and t < tclean+0.01:
+#        print 'Cleaning up Initial profile'
+#        setstage = SetValueWhere(quantity=Stage,value=vstage[0],xrange=10)
+#        setxmom  = SetValueWhere(quantity=Xmom,value=vxmom[0],xrange=10)
+#        Stage.set_values(setstage,location='centroids')
+#        Xmom.set_values(setxmom,location='centroids')
+#        Dnew = Dirichlet_boundary([vstage[0] , vxmom[0], 0.0])
+#        domain.set_boundary({'left': Dnew, 'right': T, 'top': R, 'bottom': R})
+
+print 'That took %.2f seconds' %(time.time()-t0)
+
+#FIXME: Compute average water depth on either side of shock and compare
+#to expected values. And also Froude numbers.
+
+
+#print "saving file?"
+#im = ImageGrab.grab()
+#im.save("ccube.eps")
Index: /anuga_validation/analytical solutions/Analytical_solution_wave_runup.py
===================================================================
--- /anuga_validation/analytical solutions/Analytical_solution_wave_runup.py	(revision 7624)
+++ /anuga_validation/analytical solutions/Analytical_solution_wave_runup.py	(revision 7624)
@@ -0,0 +1,177 @@
+"""Example of shallow water wave equation analytical solution of the
+one-dimensional Thacker and Greenspan wave run-up treated as a two-dimensional solution.
+
+   Copyright 2004
+   Christopher Zoppou, Stephen Roberts, Ole Nielsen, Duncan Gray
+   Geoscience Australia
+
+Specific methods pertaining to the 2D shallow water equation
+are imported from shallow_water
+for use with the generic finite volume framework
+
+Conserved quantities are h, uh and vh stored as elements 0, 1 and 2 in the
+numerical vector named conserved_quantities.
+"""
+
+######################
+# Module imports
+import sys
+from os import sep
+sys.path.append('..'+sep+'pyvolution')
+
+from shallow_water import Domain, Reflective_boundary, Time_boundary
+from math import sqrt, cos, sin, pi
+from mesh_factory import strang_mesh
+
+#Convenience functions
+def imag(a):
+    return a.imag
+
+def real(a):
+    return a.real
+
+
+######################
+# Domain
+# Strang_domain will search through the file and test to see if there are
+# two or three entries. Two entries are for points and three for triangles.
+
+#points, elements = strang_mesh('Run-up.pt')
+points, elements = strang_mesh('strang_7389.pt')
+domain = Domain(points, elements)
+
+domain.default_order = 2
+domain.smooth = True
+
+#Set a default tagging
+epsilon = 1.0e-12
+
+for id, face in domain.boundary:
+    domain.boundary[(id,face)] = 'external'
+    if domain.get_vertex_coordinate(id,(face+1)%3)[0] < -200.0+1.0e-10:
+        domain.boundary[(id,face)] = 'left'
+
+
+
+
+# Provide file name for storing output
+domain.store = True
+domain.format = 'sww'
+domain.set_name('run-up_second_order_again')
+
+print "Number of triangles = ", len(domain)
+
+#Reduction operation for get_vertex_values
+from anuga.pyvolution.util import mean
+domain.reduction = mean
+#domain.reduction = min  #Looks better near steep slopes
+
+#Define the boundary condition
+def stage_setup(x,t_star1):
+    vh = 0
+    alpha = 0.1
+    eta = 0.1
+    a = 1.5*sqrt(1.+0.9*eta)
+    l_0 = 200.
+    ii = complex(0,1)
+    g = 9.81
+    v_0 = sqrt(g*l_0*alpha)
+    v1 = 0.
+
+    sigma_max = 100.
+    sigma_min = -100.
+    for j in range (1,50):
+        sigma0 = (sigma_max+sigma_min)/2.
+        lambda_prime = 2./a*(t_star1/sqrt(l_0/alpha/g)+v1)
+        sigma_prime = sigma0/a
+        const = (1.-ii*lambda_prime)**2+sigma_prime**2
+
+        v1 = 8.*eta/a*imag(1./const**(3./2.) \
+            -3./4.*(1.-ii*lambda_prime)/const**(5./2.))
+
+        x1 = -v1**2/2.-a**2*sigma_prime**2/16.+eta \
+            *real(1.-2.*(5./4.-ii*lambda_prime) \
+            /const**(3./2.)+3./2.*(1.-ii*lambda_prime)**2 \
+            /const**(5./2.))
+
+        neta1 = x1 + a*a*sigma_prime**2/16.
+
+        v_star1 = v1*v_0
+        x_star1 = x1*l_0
+        neta_star1 = neta1*alpha*l_0
+        stage = neta_star1
+        z = stage - x_star1*alpha
+        uh = z*v_star1
+
+        if x_star1-x > 0:
+            sigma_max = sigma0
+        else:
+            sigma_min = sigma0
+
+        if abs(abs(sigma0)-100.) < 10:
+
+#     solution does not converge because bed is dry
+            stage = 0.
+            uh = 0.
+            z = 0.
+
+    return [stage, uh, vh]
+
+def boundary_stage(t):
+    x = -200
+    return stage_setup(x,t)
+
+
+######################
+#Initial condition
+print 'Initial condition'
+t_star1 = 0.0
+slope = -0.1
+
+#Set bed-elevation and friction(None)
+def x_slope(x,y):
+    n = x.shape[0]
+    z = 0*x
+    for i in range(n):
+        z[i] = -slope*x[i]
+    return z
+
+domain.set_quantity('elevation', x_slope)
+
+#Set the water depth
+print 'Initial water depth'
+def stage(x,y):
+    z = x_slope(x,y)
+    n = x.shape[0]
+    w = 0*x
+    for i in range(n):
+        w[i], uh, vh = stage_setup(x[i],t_star1)
+        h = w[i] - z[i]
+        if h < 0:
+            h = 0
+            w[i] = z[i]
+    return w
+
+domain.set_quantity('stage', stage)
+
+
+#####################
+#Set up boundary conditions
+Br = Reflective_boundary(domain)
+Bw = Time_boundary(domain, boundary_stage)
+domain.set_boundary({'left': Bw, 'external': Br})
+
+
+
+#domain.visualise = True
+
+######################
+#Evolution
+import time
+t0 = time.time()
+for t in domain.evolve(yieldstep = 1., finaltime = 100):
+    domain.write_time()
+    print boundary_stage(domain.time)
+
+print 'That took %.2f seconds' %(time.time()-t0)
+
