Index: /anuga_core/source/anuga/culvert_flows/test_culvert_class.py
===================================================================
--- /anuga_core/source/anuga/culvert_flows/test_culvert_class.py	(revision 7495)
+++ /anuga_core/source/anuga/culvert_flows/test_culvert_class.py	(revision 7496)
@@ -22,4 +22,105 @@
 
 import numpy as num
+
+
+# Helper functions
+def run_culvert_flow_problem(depth):
+    """Run flow with culvert given depth
+    """
+
+    length = 40.
+    width = 5.
+
+    dx = dy = 1           # Resolution: Length of subdivisions on both axes
+
+    points, vertices, boundary = rectangular_cross(int(length/dx),
+                                                   int(width/dy),
+                                                   len1=length, 
+                                                   len2=width)
+    domain = Domain(points, vertices, boundary)   
+    domain.set_name('Test_culvert_shallow') # Output name
+    domain.set_default_order(2)
+
+
+    #----------------------------------------------------------------------
+    # Setup initial conditions
+    #----------------------------------------------------------------------
+
+    def topography(x, y):
+        """Set up a weir
+        
+        A culvert will connect either side
+        """
+        # General Slope of Topography
+        z=-x/1000
+        
+        N = len(x)
+        for i in range(N):
+
+           # Sloping Embankment Across Channel
+            if 5.0 < x[i] < 10.1:
+                # Cut Out Segment for Culvert face                
+                if  1.0+(x[i]-5.0)/5.0 <  y[i]  < 4.0 - (x[i]-5.0)/5.0: 
+                   z[i]=z[i]
+                else:
+                   z[i] +=  0.5*(x[i] -5.0)    # Sloping Segment  U/S Face
+            if 10.0 < x[i] < 12.1:
+               z[i] +=  2.5                    # Flat Crest of Embankment
+            if 12.0 < x[i] < 14.5:
+                # Cut Out Segment for Culvert face                
+                if  2.0-(x[i]-12.0)/2.5 <  y[i]  < 3.0 + (x[i]-12.0)/2.5:
+                   z[i]=z[i]
+                else:
+                   z[i] +=  2.5-1.0*(x[i] -12.0) # Sloping D/S Face
+            		   
+    		
+        return z
+
+
+    domain.set_quantity('elevation', topography) 
+    domain.set_quantity('friction', 0.01)         # Constant friction 
+    domain.set_quantity('stage',
+                        expression='elevation + %f' % depth) # Shallow initial condition
+
+    # Boyd culvert
+    culvert = Culvert_flow(domain,
+                           label='Culvert No. 1',
+                           description='This culvert is a test unit 1.2m Wide by 0.75m High',   
+                           end_point0=[9.0, 2.5], 
+                           end_point1=[13.0, 2.5],
+                           width=1.20, height=0.75,
+                           culvert_routine=boyd_generalised_culvert_model,
+                           number_of_barrels=1,
+                           update_interval=2,
+                           verbose=False)
+    
+
+    domain.forcing_terms.append(culvert)
+    
+
+    #-----------------------------------------------------------------------
+    # Setup boundary conditions
+    #-----------------------------------------------------------------------
+
+    # Inflow based on Flow Depth and Approaching Momentum
+
+    Br = Reflective_boundary(domain)              # Solid reflective wall
+    domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
+
+
+    
+    #-----------------------------------------------------------------------
+    # Evolve system through time
+    #-----------------------------------------------------------------------
+
+    #print 'depth', depth
+    ref_volume = domain.get_quantity('stage').get_integral()
+    for t in domain.evolve(yieldstep = 0.1, finaltime = 25):
+        new_volume = domain.get_quantity('stage').get_integral()
+        
+        msg = ('Total volume has changed: Is %.8f m^3 should have been %.8f m^3'
+               % (new_volume, ref_volume))
+        assert num.allclose(new_volume, ref_volume), msg        
+        
 
 
@@ -100,4 +201,5 @@
                                end_point0=[9.0, 2.5], 
                                end_point1=[13.0, 2.5],
+                               width=1.00,
                                use_velocity_head=True,
                                verbose=False)
@@ -214,4 +316,5 @@
                                end_point0=[9.0, 2.5], 
                                end_point1=[13.0, 2.5],
+                               height=0.75,
                                verbose=False)
 
@@ -243,5 +346,22 @@
 
     
-    def test_that_culvert_rating_limits_flow_in_shallow_inlet_condition(self):
+
+        
+    def test_that_culvert_flows_conserves_volume(self):
+        """test_that_culvert_flows_conserves_volume
+
+        Test that culvert on a sloping dry bed limits flows when very little water
+        is present at inlet.
+
+        Uses helper function: run_culvert_flow_problem(depth):
+
+        """
+
+        # Try this for a range of depths
+        for depth in [0.1, 0.2, 0.5, 1.0]:
+            run_culvert_flow_problem(depth)
+
+
+    def OBSOLETE_XXXtest_that_culvert_rating_limits_flow_in_shallow_inlet_condition(self):
         """test_that_culvert_rating_limits_flow_in_shallow_inlet_condition
         
@@ -251,4 +371,6 @@
         This one is using the rating curve variant
         """
+
+        
 
         path = get_pathname_from_package('anuga.culvert_flows')    
@@ -307,13 +429,25 @@
         domain.set_quantity('stage',
                             expression='elevation + 0.1') # Shallow initial condition
-                            
-
-        filename = os.path.join(path, 'example_rating_curve.csv')
+
+        # Boyd culvert
         culvert = Culvert_flow(domain,
-                               culvert_description_filename=filename,
+                               label='Culvert No. 1',
+                               description='This culvert is a test unit 1.2m Wide by 0.75m High',   
                                end_point0=[9.0, 2.5], 
                                end_point1=[13.0, 2.5],
-                               trigger_depth=0.01,
+                               width=1.20, height=0.75,
+                               culvert_routine=boyd_generalised_culvert_model,
+                               number_of_barrels=1,
+                               update_interval=2,
                                verbose=False)
+        
+        # Rating curve
+        #filename = os.path.join(path, 'example_rating_curve.csv')
+        #culvert = Culvert_flow(domain,
+        #                       culvert_description_filename=filename,
+        #                       end_point0=[9.0, 2.5], 
+        #                       end_point1=[13.0, 2.5],
+        #                       trigger_depth=0.01,
+        #                       verbose=False)
 
         domain.forcing_terms.append(culvert)
@@ -335,9 +469,9 @@
         #-----------------------------------------------------------------------
 
+        print 'depth', 0.1
         ref_volume = domain.get_quantity('stage').get_integral()
         for t in domain.evolve(yieldstep = 0.1, finaltime = 25):
             new_volume = domain.get_quantity('stage').get_integral()
 
-
             msg = ('Total volume has changed: Is %.8f m^3 should have been %.8f m^3'
                    % (new_volume, ref_volume))
@@ -346,8 +480,11 @@
         
         return
-        # Now try this for a range of other depths
-        for depth in [1.0, 0.5, 0.2, 0.1, 0.05]:
+        # Now try this again for a depth of 10 cm and for a range of other depths
+        for depth in [0.1, 0.2, 0.5, 1.0]:
+            print 'depth', depth
             domain.set_time(0.0)
             
+            domain.set_quantity('elevation', topography) 
+            domain.set_quantity('friction', 0.01)         # Constant friction 
             domain.set_quantity('stage',
                                 expression='elevation + %f' % depth)
@@ -358,5 +495,5 @@
                 new_volume = domain.get_quantity('stage').get_integral()
             
-                msg = 'Total volume has changed: Is %.2f m^3 should have been %.2f m^3'\
+                msg = 'Total volume has changed: Is %.8f m^3 should have been %.8f m^3'\
                     % (new_volume, ref_volume)
 
@@ -439,9 +576,9 @@
                                end_point0=[9.0, 2.5], 
                                end_point1=[13.0, 2.5],
-                               width=1.20,height=0.75,
+                               width=1.20, height=0.75,
                                culvert_routine=boyd_generalised_culvert_model,
                                number_of_barrels=1,
                                update_interval=2,
-                               verbose=True)
+                               verbose=False)
         
         domain.forcing_terms.append(culvert)
@@ -543,5 +680,5 @@
                                culvert_routine=boyd_generalised_culvert_model,        
                                number_of_barrels=1,
-                               verbose=True)
+                               verbose=False)
                                
         
@@ -551,4 +688,113 @@
         culvert(domain)
     
+
+    
+
+    def test_momentum_jet(self):
+        """test_momentum_jet
+        
+        Test that culvert_class can accept keyword use_momentum_jet
+        This does not yet imply that the values have been tested. FIXME
+        """
+
+
+        length = 40.
+        width = 5.
+
+        dx = dy = 1           # Resolution: Length of subdivisions on both axes
+
+        points, vertices, boundary = rectangular_cross(int(length/dx),
+                                                       int(width/dy),
+                                                       len1=length, 
+                                                       len2=width)
+        domain = Domain(points, vertices, boundary)   
+        domain.set_name('Test_culvert_shallow') # Output name
+        domain.set_default_order(2)
+
+
+        #----------------------------------------------------------------------
+        # Setup initial conditions
+        #----------------------------------------------------------------------
+
+        def topography(x, y):
+            """Set up a weir
+            
+            A culvert will connect either side
+            """
+            # General Slope of Topography
+            z=-x/1000
+            
+            N = len(x)
+            for i in range(N):
+
+               # Sloping Embankment Across Channel
+                if 5.0 < x[i] < 10.1:
+                    # Cut Out Segment for Culvert face                
+                    if  1.0+(x[i]-5.0)/5.0 <  y[i]  < 4.0 - (x[i]-5.0)/5.0: 
+                       z[i]=z[i]
+                    else:
+                       z[i] +=  0.5*(x[i] -5.0)    # Sloping Segment  U/S Face
+                if 10.0 < x[i] < 12.1:
+                   z[i] +=  2.5                    # Flat Crest of Embankment
+                if 12.0 < x[i] < 14.5:
+                    # Cut Out Segment for Culvert face                
+                    if  2.0-(x[i]-12.0)/2.5 <  y[i]  < 3.0 + (x[i]-12.0)/2.5:
+                       z[i]=z[i]
+                    else:
+                       z[i] +=  2.5-1.0*(x[i] -12.0) # Sloping D/S Face
+                		   
+        		
+            return z
+
+
+        domain.set_quantity('elevation', topography) 
+        domain.set_quantity('friction', 0.01)         # Constant friction 
+        domain.set_quantity('stage',
+                            expression='elevation + 1.0') # Shallow initial condition
+
+        # Boyd culvert
+        culvert = Culvert_flow(domain,
+                               label='Culvert No. 1',
+                               description='This culvert is a test unit 1.2m Wide by 0.75m High',   
+                               end_point0=[9.0, 2.5], 
+                               end_point1=[13.0, 2.5],
+                               width=1.20, height=0.75,
+                               culvert_routine=boyd_generalised_culvert_model,
+                               number_of_barrels=1,
+                               use_momentum_jet=True,
+                               update_interval=2,
+                               verbose=False)
+    
+
+        domain.forcing_terms.append(culvert)
+
+        
+        # Call
+        culvert(domain)
+    
+
+        #-----------------------------------------------------------------------
+        # Setup boundary conditions
+        #-----------------------------------------------------------------------
+
+
+        Br = Reflective_boundary(domain)              # Solid reflective wall
+        domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
+
+        #-----------------------------------------------------------------------
+        # Evolve system through time
+        #-----------------------------------------------------------------------
+
+        ref_volume = domain.get_quantity('stage').get_integral()
+        for t in domain.evolve(yieldstep = 0.1, finaltime = 25):
+            new_volume = domain.get_quantity('stage').get_integral()
+            
+            msg = ('Total volume has changed: Is %.8f m^3 should have been %.8f m^3'
+                   % (new_volume, ref_volume))
+            assert num.allclose(new_volume, ref_volume), msg        
+        
+
+        
+
                
 #-------------------------------------------------------------
