Index: anuga_work/development/ted_rigby/gauge_locations.csv
===================================================================
--- anuga_work/development/ted_rigby/gauge_locations.csv	(revision 7475)
+++ anuga_work/development/ted_rigby/gauge_locations.csv	(revision 7475)
@@ -0,0 +1,34 @@
+name     ,easting    ,northing   ,elevation
+Gauge0_at2m,2.0,0.0,10.00
+Gauge1_at2m,2.0,1.0,9.00
+Gauge2_at2m,2.0,2.0,8.00
+Gauge3_at2m,2.0,3.0,7.00
+Gauge4_at2m,2.0,4.0,6.00
+Gauge5_at2m,2.0,5.0,5.00
+Gauge6_at2m,2.0,6.0,4.00
+Gauge7_at2m,2.0,7.0,3.00
+Gauge8_at2m,2.0,8.0,2.00
+Gauge9_at2m,2.0,9.0,1.00
+Gauge10_at2m,2.0,10.0,0.00
+Gauge0_at500m,500.0,0.0,9.50
+Gauge1_at500m,500.0,1.0,8.50
+Gauge2_at500m,500.0,2.0,7.50
+Gauge3_at500m,500.0,3.0,6.50
+Gauge4_at500m,500.0,4.0,5.50
+Gauge5_at500m,500.0,5.0,4.50
+Gauge6_at500m,500.0,6.0,3.50
+Gauge7_at500m,500.0,7.0,2.50
+Gauge8_at500m,500.0,8.0,1.50
+Gauge9_at500m,500.0,9.0,0.50
+Gauge10_at500m,500.0,10.0,-0.50
+Gauge0_at998m,998.0,0.0,9.00
+Gauge1_at998m,998.0,1.0,8.00
+Gauge2_at998m,998.0,2.0,7.00
+Gauge3_at998m,998.0,3.0,6.00
+Gauge4_at998m,998.0,4.0,5.00
+Gauge5_at998m,998.0,5.0,4.00
+Gauge6_at998m,998.0,6.0,3.00
+Gauge7_at998m,998.0,7.0,2.00
+Gauge8_at998m,998.0,8.0,1.00
+Gauge9_at998m,998.0,9.0,0.00
+Gauge10_at998m,998.0,10.0,-1.0
Index: anuga_work/development/ted_rigby/plane_gauge_locations.csv
===================================================================
--- anuga_work/development/ted_rigby/plane_gauge_locations.csv	(revision 7475)
+++ anuga_work/development/ted_rigby/plane_gauge_locations.csv	(revision 7475)
@@ -0,0 +1,12 @@
+name     ,easting    ,northing   ,elevation
+EDGE_200m,200.0,0.0,0.00
+1_200m,200.0,1.0,0.00
+2_200m,200.0,1.0,0.00
+3_200m,200.0,1.0,0.00
+4_200m,200.0,1.0,0.00
+5_200m,200.0,1.0,0.00
+6_200m,200.0,1.0,0.00
+7_200m,200.0,1.0,0.00
+8_200m,200.0,1.0,0.00
+9_200m,200.0,1.0,0.00
+CL_200m,200.0,1.0,0.00
Index: anuga_work/development/ted_rigby/test_flowhydraulics_in_vee_using_inflow.py
===================================================================
--- anuga_work/development/ted_rigby/test_flowhydraulics_in_vee_using_inflow.py	(revision 7475)
+++ anuga_work/development/ted_rigby/test_flowhydraulics_in_vee_using_inflow.py	(revision 7475)
@@ -0,0 +1,147 @@
+"""test_flow_in_vee_using_inflow.py
+Test the ability of ANUGA to 
+a) maintain the flow rate input into the channel at the upstream end
+b) replicate normal depth half way down the channel after half an hour
+c) replicate a typical velocity profile across the channel after half an hour
+
+The vee channel is 20m wide by 1000m long with side slopes of 1:1 
+The equations for dn are generalised with side slope as 1 V in Side_Slope H
+The test location is in the invert of (y=10m)and half way down the channel (x=500m) 
+This script explores flow behaviour predicted by ANUGA and compares it with known solutions
+
+This test inputs flow using the 'inflow' function .
+
+"""
+
+verbose = True
+import numpy as num
+
+#------------------------------------------------------------------------------
+# Import necessary modules
+#------------------------------------------------------------------------------
+from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular_cross
+from anuga.shallow_water import Domain
+from anuga.shallow_water.shallow_water_domain import Reflective_boundary
+from anuga.shallow_water.shallow_water_domain import Dirichlet_boundary
+from anuga.shallow_water.shallow_water_domain import Transmissive_Momentum_Set_Stage_boundary
+from anuga.shallow_water.shallow_water_domain import Transmissive_boundary
+from anuga.shallow_water.shallow_water_domain import Inflow
+from anuga.shallow_water.data_manager import get_flow_through_cross_section
+from anuga.abstract_2d_finite_volumes.util import sww2csv_gauges, csv2timeseries_graphs
+
+
+#------------------------------------------------------------------------------
+# Setup computational domain
+#------------------------------------------------------------------------------
+
+finaltime = 1800.0    # 1/2 hour
+
+length     = 1000.
+width      = 20.
+side_slope = 1.0
+dx = dy    = 1          # Resolution: of points (elev) grid on both axes
+ref_flow   = 50.0       # Inflow at head of vee
+slope      = 1.0/1000
+mannings_n = 0.150
+
+points, vertices, boundary = rectangular_cross(int(length/dx), int(width/dy),
+                                               len1=length, len2=width)
+
+
+
+domain = Domain(points, vertices, boundary)   
+domain.set_name('test_flowhydraulics_in_vee_using_inflow')  # Output name
+
+
+#------------------------------------------------------------------------------
+# Setup initial conditions
+#------------------------------------------------------------------------------
+
+def topography(x, y):
+    z=-x * slope + abs((y-10)* side_slope)
+    return z
+
+domain.set_quantity('elevation', topography)  # Use function for elevation
+domain.set_quantity('friction', mannings_n)   # Constant friction of conc surface   
+domain.set_quantity('stage',
+                    expression='elevation')   # zero depth initial condition
+
+
+#------------------------------------------------------------------------------
+# Compute normal depth flow characteristics for vee using Mannings equation
+#------------------------------------------------------------------------------
+
+normal_depth  = ( ref_flow*mannings_n* ( (2.0*( (1.0+(1.0/side_slope)**2.0)**0.5) )**0.6667 ) / (side_slope*(slope**0.5))  )**0.375
+area          = (normal_depth**2.0) *side_slope
+wet_perimeter = 2* normal_depth*((1+(side_slope)**2.0)**0.5)
+hyd_radius    = area/wet_perimeter
+flow_velocity = ref_flow/area
+top_width     = 2* side_slope*normal_depth
+
+print 'Mannings solution for normal depth flow in vee is; '
+print '    Channel sideslopes = ', side_slope
+print '    Channel n          = ', mannings_n
+print '    Channel slope      = ', slope
+print '    Reference flow     = ', ref_flow
+print '    Normal depth       = ', normal_depth
+print '    Flow area          = ', area
+print '    Wet perimeter      = ', wet_perimeter
+print '    Hydraulic radius   = ', hyd_radius
+print '    Avg velocity       = ', flow_velocity
+print '    Flow width         = ', top_width
+
+#------------------------------------------------------------------------------
+# Calculate inflow and append to domain forcing terms
+#------------------------------------------------------------------------------
+
+# Fixed Flowrate onto Area 
+fixed_inflow = Inflow(domain,
+                      center=(10.0, 10.0),
+                      radius=5.00,
+                      rate=ref_flow)   
+domain.forcing_terms.append(fixed_inflow)
+
+#------------------------------------------------------------------------------
+# Setup boundary conditions
+#------------------------------------------------------------------------------
+
+
+Br = Reflective_boundary(domain)                       # Solid reflective wall on edge of sides
+# set outflow depth to mannings normal depth for flow to avoid drawdown impact upstream
+def outflow_depth(t):
+    return (-slope*length) + normal_depth 
+Bo = Transmissive_Momentum_Set_Stage_boundary(domain=domain,function=outflow_depth)
+
+domain.set_boundary({'left': Br, 'right': Bo, 'top': Br, 'bottom': Br})
+
+#------------------------------------------------------------------------------
+# Evolve system through time
+#------------------------------------------------------------------------------
+
+for t in domain.evolve(yieldstep=100.0, finaltime=finaltime):
+    if verbose :
+        print domain.timestepping_statistics()
+
+
+#------------------------------------------------------------------------------
+# Compute flow across flowlines -v- ref_flow 
+#------------------------------------------------------------------------------
+
+# Check flow across flowline at 500m (where  it should have normalised)
+q=domain.get_flow_through_cross_section([[500.0,0.0],[500.0,width]])
+print '90 degree flowline half way down channel: ANUGA Q = %f, Ref_flow = %f' %(q, ref_flow)
+
+# Check flow across flowline at 998m (wimmediately US of outflow bdry)
+q=domain.get_flow_through_cross_section([[998.0,0.0],[998.0,width]])
+print '90 degree flowline just US of outflow bdry: ANUGA Q = %f, Ref_flow = %f' %(q, ref_flow)
+
+#------------------------------------------------------------------------------
+# Record temporal flow characteristics just ds inflow bdry
+#------------------------------------------------------------------------------
+
+# Store temporal pattern of gauges at head/middle/bottom of channel in csv files 
+sww2csv_gauges(sww_file  = 'test_flowhydraulics_in_vee_using_inflow.sww',
+               gauge_file= 'gauge_locations.csv',
+               out_name  = 'inflow_',
+               quantities=['stage','speed','elevation'],
+               verbose=True)
Index: anuga_work/development/ted_rigby/test_flowhydraulics_in_vee_using_momentum.py
===================================================================
--- anuga_work/development/ted_rigby/test_flowhydraulics_in_vee_using_momentum.py	(revision 7475)
+++ anuga_work/development/ted_rigby/test_flowhydraulics_in_vee_using_momentum.py	(revision 7475)
@@ -0,0 +1,155 @@
+"""test_flow_in_vee_using_momentum.py
+
+Test the ability of ANUGA to 
+a) maintain the flow rate input into the channel at the upstream end
+b) replicate normal depth half way down the channel after half an hour
+c) replicate a typical velocity profile across the channel after half an hour
+
+The vee channel is 20m wide by 1000m long with side slopes of 1:1 
+The equations for dn are generalised with side slope as 1 V in Side_Slope H
+The test location is in the invert of (y=10m)and half way down the channel (x=500m) 
+This script explores flow behaviour predicted by ANUGA and compares it with known solutions
+
+This test inputs flow using momentum=flow/width .
+
+"""
+
+verbose = True
+import numpy as num
+
+#------------------------------------------------------------------------------
+# Import necessary modules
+#------------------------------------------------------------------------------
+from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular_cross
+from anuga.shallow_water import Domain
+from anuga.shallow_water.shallow_water_domain import Reflective_boundary
+from anuga.shallow_water.shallow_water_domain import Dirichlet_boundary
+from anuga.shallow_water.shallow_water_domain import Transmissive_Momentum_Set_Stage_boundary
+from anuga.shallow_water.shallow_water_domain import Transmissive_boundary
+from anuga.shallow_water.shallow_water_domain import Inflow
+from anuga.shallow_water.data_manager import get_flow_through_cross_section
+from anuga.abstract_2d_finite_volumes.util import sww2csv_gauges, csv2timeseries_graphs
+
+
+#------------------------------------------------------------------------------
+# Setup computational domain
+#------------------------------------------------------------------------------
+
+finaltime = 1800.0    # 1/2 hour
+
+length     = 1000.
+width      = 20.
+side_slope = 1.0
+dx = dy    = 1          # Resolution: of points (elev) grid on both axes
+ref_flow   = 20.0       # Inflow at head of vee
+slope      = 1.0/1000
+mannings_n = 0.150
+
+points, vertices, boundary = rectangular_cross(int(length/dx), int(width/dy),
+                                               len1=length, len2=width)
+
+
+
+domain = Domain(points, vertices, boundary)   
+domain.set_name('test_flowhydraulics_in_vee_using_momentum') # Output name
+
+
+#------------------------------------------------------------------------------
+# Setup initial conditions
+#------------------------------------------------------------------------------
+
+def topography(x, y):
+    z=-x * slope + abs((y-10)* side_slope)
+    return z
+
+domain.set_quantity('elevation', topography)  # Use function for elevation
+domain.set_quantity('friction', mannings_n)   # Constant friction of conc surface   
+domain.set_quantity('stage',
+                    expression='elevation')   # zero depth initial condition
+
+
+
+#------------------------------------------------------------------------------
+# Compute normal depth flow characteristics for vee using Mannings equation
+#------------------------------------------------------------------------------
+
+normal_depth  = ( ref_flow*mannings_n* ( (2.0*( (1.0+(1.0/side_slope)**2.0)**0.5) )**0.6667 ) / (side_slope*(slope**0.5))  )**0.375
+area          = (normal_depth**2.0) *side_slope
+wet_perimeter = 2* normal_depth*((1+(side_slope)**2.0)**0.5)
+hyd_radius    = area/wet_perimeter
+flow_velocity = ref_flow/area
+top_width     = 2* side_slope*normal_depth
+
+print 'Mannings solution for normal depth flow in vee is; '
+print '    Channel sideslopes = ', side_slope
+print '    Channel n          = ', mannings_n
+print '    Channel slope      = ', slope
+print '    Reference flow     = ', ref_flow
+print '    Normal depth       = ', normal_depth
+print '    Flow area          = ', area
+print '    Wet perimeter      = ', wet_perimeter
+print '    Hydraulic radius   = ', hyd_radius
+print '    Avg velocity       = ', flow_velocity
+print '    Flow width         = ', top_width
+
+
+#------------------------------------------------------------------------------
+# Calculate xmomentum and depth to apply to the inflow boundary
+#------------------------------------------------------------------------------
+
+# NOTE: momentum is applied 'uniformly' to the 'full' inflow boundary by ANUGA
+# this is all very unrealistic as only part of the bdry is 'wet'
+# but should reasonably replicate ref_flow in ANUGA nevertheless
+inflow_xmomentum = ref_flow/width  
+inflow_depth     = normal_depth
+
+#------------------------------------------------------------------------------
+# Setup boundary conditions
+#------------------------------------------------------------------------------
+
+
+Br = Reflective_boundary(domain)       # Solid reflective wall on edge of sides
+Bi = Dirichlet_boundary([inflow_depth,inflow_xmomentum,0.0])    # inflow at top
+# set outflow depth to mannings normal depth for flow to avoid drawdown impact upstream
+def outflow_depth(t):
+    return (-slope*length) + normal_depth 
+Bo = Transmissive_Momentum_Set_Stage_boundary(domain=domain,function=outflow_depth)
+
+
+domain.set_boundary({'left': Bi, 'right': Bo, 'top': Br, 'bottom': Br})
+
+#------------------------------------------------------------------------------
+# Evolve system through time
+#------------------------------------------------------------------------------
+
+for t in domain.evolve(yieldstep=100.0, finaltime=finaltime):
+    if verbose :
+        print domain.timestepping_statistics()
+
+
+#------------------------------------------------------------------------------
+# Compute flow across flowlines -v- ref_flow 
+#------------------------------------------------------------------------------
+
+# Check flow across flowline at 2m (immediately DS of inflow bdry)
+q=domain.get_flow_through_cross_section([[2.0,0.0],[2.0,width]])
+print '90 degree flowline just ds of inflow bdry: ANUGA Q = %f, Ref_flow = %f' %(q, ref_flow)
+
+# Check flow across flowline at 500m (where  it should have normalised)
+q=domain.get_flow_through_cross_section([[500.0,0.0],[500.0,width]])
+print '90 degree flowline half way down channel: ANUGA Q = %f, Ref_flow = %f' %(q, ref_flow)
+
+# Check flow across flowline at 998m (wimmediately US of outflow bdry)
+q=domain.get_flow_through_cross_section([[998.0,0.0],[998.0,width]])
+print '90 degree flowline just US of outflow bdry: ANUGA Q = %f, Ref_flow = %f' %(q, ref_flow)
+
+#------------------------------------------------------------------------------
+# Record temporal flow characteristics just ds inflow bdry
+#------------------------------------------------------------------------------
+
+# Store temporal pattern of gauges at head/middle/bottom of channel in csv files 
+sww2csv_gauges(sww_file  = 'test_flowhydraulics_in_vee_using_momentum.sww',
+               gauge_file= 'gauge_locations.csv',
+               out_name  = 'momentum_',
+               quantities=['stage','speed','elevation'],
+               verbose=True)
Index: anuga_work/development/ted_rigby/test_flowhydraulics_on_plane.py
===================================================================
--- anuga_work/development/ted_rigby/test_flowhydraulics_on_plane.py	(revision 7475)
+++ anuga_work/development/ted_rigby/test_flowhydraulics_on_plane.py	(revision 7475)
@@ -0,0 +1,134 @@
+
+
+"""test_flow_on_plane.py
+Test the ability of ANUGA to maintain flow and match normal depth on an inclined plane
+using constant inflow thru the upstream bdry at normal_depth and equivalent momentum
+The plane is 20m wide by 400m long dipping at various slopes with a range of roughnesses 
+A recording gauges is located  200m downstream
+"""
+
+verbose = True
+import numpy as num
+
+#------------------------------------------------------------------------------
+# Import necessary modules
+#------------------------------------------------------------------------------
+from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular_cross
+from anuga.shallow_water import Domain
+from anuga.shallow_water.shallow_water_domain import Reflective_boundary
+from anuga.shallow_water.shallow_water_domain import Dirichlet_boundary
+from anuga.shallow_water.shallow_water_domain import Transmissive_Momentum_Set_Stage_boundary
+from anuga.shallow_water.shallow_water_domain import Inflow
+from anuga.shallow_water.data_manager import get_flow_through_cross_section
+from anuga.abstract_2d_finite_volumes.util import sww2csv_gauges, csv2timeseries_graphs
+
+
+#------------------------------------------------------------------------------
+# Setup computational domain
+#------------------------------------------------------------------------------
+
+finaltime = 1800.0    # 1/2 hour
+
+length   = 400.
+width    = 20.
+dx = dy  = 2          # Resolution: of grid on both axes
+ref_flow = 20.0       # Inflow at head of plane
+
+points, vertices, boundary = rectangular_cross(int(length/dx), int(width/dy),
+                                               len1=length, len2=width)
+
+for mannings_n in [0.150, 0.070, 0.03, 0.015]:
+    # Loop over a range of roughnesses  
+    
+    for slope in [1.0/1000, 1.0/300, 1.0/100]:
+        # Loop over a range of bedslopes  
+        
+
+        domain = Domain(points, vertices, boundary)   
+        domain.set_name('test_flowhydraulics_on_plane')              # Output name
+
+
+        #------------------------------------------------------------------------------
+        # Setup initial conditions
+        #------------------------------------------------------------------------------
+
+        def topography(x, y):
+            z=-x * slope
+            return z
+
+        domain.set_quantity('elevation', topography)  # Use function for elevation
+        domain.set_quantity('friction', mannings_n)   # Constant friction of conc surface   
+        domain.set_quantity('stage',
+                            expression='elevation')   # Dry initial condition
+
+
+        #------------------------------------------------------------------------------
+        # Seup Inflow and compute flow characteristics
+        #------------------------------------------------------------------------------
+      
+        # Compute normal depth, Froude and momentum on plane using Mannings equation
+        normal_depth  =(ref_flow*mannings_n/(slope**0.5*width))**0.6
+        flow_velocity = ref_flow/(width*normal_depth)
+        xmomentum     = flow_velocity*normal_depth
+        froude_no     = flow_velocity/(9.8*normal_depth)**0.5
+        if verbose:
+            print
+            print 'Slope:', slope, 'Mannings n:', mannings_n, 'Velocity:',flow_velocity,'Froude:',froude_no
+
+
+        #------------------------------------------------------------------------------
+        # Setup boundary conditions
+        #------------------------------------------------------------------------------
+
+        Br = Reflective_boundary(domain)                       # Solid reflective wall
+        Bo = Dirichlet_boundary([-(slope*length)+normal_depth, xmomentum, 0])  # Outflow stage at normal depth
+        Bi = Dirichlet_boundary([normal_depth, xmomentum, 0])  # Inflow stage at normal depth
+       
+        domain.set_boundary({'left': Bi, 'right': Bo, 'top': Br, 'bottom': Br})
+
+        #------------------------------------------------------------------------------
+        # Evolve system through time
+        #------------------------------------------------------------------------------
+        
+        for t in domain.evolve(yieldstep=100.0, finaltime=finaltime):
+            if verbose :
+                print domain.timestepping_statistics()
+
+
+        #------------------------------------------------------------------------------
+        # Compute depth and check final depth against normal depth
+        #------------------------------------------------------------------------------
+
+        
+        # Store temporal pattern at 200m in a csv file to check for steady state free of waves
+        sww2csv_gauges('test_flowhydraulics_on_plane.sww','plane_gauge_locations.csv',quantities=['stage','elevation'],verbose=True)
+
+        # Compute domain depth in middle of plane at 200m
+        x=200.0 
+        y=10.00
+        w = domain.get_quantity('stage').get_values(interpolation_points=[[x, y]])[0]
+        z = domain.get_quantity('elevation').get_values(interpolation_points=[[x, y]])[0]
+        domain_depth = w-z
+        
+        msg = 'Predicted depth of flow at 200m was %f, should have been %f' % (domain_depth, normal_depth)
+        if verbose:
+            print 'Depth at 200m: ANUGA = %f, Mannings = %f' % (domain_depth, normal_depth)
+        assert num.allclose(domain_depth,normal_depth, rtol=5.0e-2), msg
+
+        #------------------------------------------------------------------------------
+        # Compute flow thru flowlines at 200m - confirm equals inflow
+        #------------------------------------------------------------------------------
+
+        # Square on flowline at 200m
+        q=domain.get_flow_through_cross_section([[200.0,0.0],[200.0,20.0]])
+        msg = 'Predicted square on flow was %f, should have been %f' % (q, ref_flow)
+        if verbose:
+            print '90 degree flowline: ANUGA = %f, Ref = %f' % (q, ref_flow)
+        assert num.allclose(q, ref_flow, rtol=1.0e-2), msg         
+        
+        # 45 degree flowline at 200m
+        q=domain.get_flow_through_cross_section([[200.0,0.0],[220.0,20.0]])
+        msg = 'Predicted 45 flow was %f, should have been %f' % (q, ref_flow)
+        if verbose:
+            print '45 degree flowline: ANUGA = %f, Ref = %f' % (q, ref_flow)
+        assert num.allclose(q, ref_flow, rtol=1.0e-2), msg         
