Changeset 9740 for trunk/anuga_core/anuga/structures

Timestamp:

Jan 23, 2017, 4:44:46 PM (8 years ago)

Author:

steve

Message:

Update to svn

Location:

trunk/anuga_core/anuga/structures

Files:

• boyd_box_operator.py (modified) (5 diffs)
• boyd_pipe_operator.py (modified) (13 diffs)
• tests/test_boyd_box_operator.py (modified) (1 diff)
• tests/test_boyd_pipe_operator.py (modified) (11 diffs)
• tests/test_weir_orifice_trapezoid_operator.py (modified) (20 diffs)
• weir_orifice_trapezoid_operator.py (modified) (10 diffs)

trunk/anuga_core/anuga/structures/boyd_box_operator.py

@@ -28 +28 @@
                  width,                                    
                  height=None,                        
-                 blockage=0.0,  # added by DPM 24/7/2016 
+                 blockage=0.0,
                  z1=0.0,
                  z2=0.0,                 
@@ -55 +55 @@
                                           width=width,
                                           height=height,
-                                          blockage=blockage,  # added by DPM 24/7/2016
+                                          blockage=blockage,
                                           diameter=None,                                         
                                           apron=apron,
@@ -83 +83 @@
         self.culvert_width = self.get_culvert_width()
         self.culvert_height = self.get_culvert_height()
-        self.culvert_blockage = self.get_culvert_blockage()#added DPM 24/7/2016
+        self.culvert_blockage = self.get_culvert_blockage()
 
         #FIXME SR: Why is this hard coded!
@@ -196 +196 @@
                 print 'width ',self.culvert_width
                 print 'depth ',self.culvert_height
-                print 'culvert blockage ',self.culvert_blockage #added DPM 24/7/2016
+                print 'culvert blockage ',self.culvert_blockage
                 print 'flow_width ',self.culvert_width
                 print 'length ' ,self.culvert_length
@@ -275 +275 @@
 
     local_debug = False
-           
-    Q_inlet_unsubmerged = 0.544*anuga.g**0.5*(1-blockage)*width*driving_energy**1.50 # Flow based on Inlet Ctrl Inlet Unsubmerged
-    Q_inlet_submerged = 0.702*anuga.g**0.5*(1-blockage)**2*width*depth**0.89*driving_energy**0.61  # Flow based on Inlet Ctrl Inlet Submerged
-
+    
+    if blockage >= 1.0:
+        Q = barrel_velocity = outlet_culvert_depth = 0.0
+        flow_area = 0.00001
+        case = '100 blocked culvert'
+        return Q, barrel_velocity, outlet_culvert_depth, flow_area, case
+    else:                      
+        Q_inlet_unsubmerged = 0.544*anuga.g**0.5*(1-blockage)*width*driving_energy**1.50 # Flow based on Inlet Ctrl Inlet Unsubmerged
+        Q_inlet_submerged = 0.702*anuga.g**0.5*(1-blockage)**2*width*depth**0.89*driving_energy**0.61  # Flow based on Inlet Ctrl Inlet Submerged
+
+    #print 'blockage ', blockage
     # FIXME(Ole): Are these functions really for inlet control?
     if Q_inlet_unsubmerged < Q_inlet_submerged:

trunk/anuga_core/anuga/structures/boyd_pipe_operator.py

@@ -26 +26 @@
                  losses,
                  diameter=None,
+                 blockage=0.0,
                  z1=0.0,
                  z2=0.0,
@@ -52 +53 @@
                                           height=None,
                                           diameter=diameter,
+                                          blockage=blockage,
                                           apron=apron,
                                           manning=manning,
@@ -77 +79 @@
         self.culvert_length = self.get_culvert_length()
         self.culvert_diameter = self.get_culvert_diameter()
+        self.culvert_blockage = self.get_culvert_blockage()
 
         #print self.culvert_diameter
@@ -136 +139 @@
                               boyd_pipe_function(depth               =self.inflow.get_enquiry_depth(),
                                                  diameter            =self.culvert_diameter,
+                                                 blockage            =self.culvert_blockage,
                                                  length              =self.culvert_length,
                                                  driving_energy      =self.driving_energy,
@@ -143 +147 @@
                                                  manning             =self.manning)
         else:
-             Q = barrel_velocity = outlet_culvert_depth = 0.0
-             case = 'Inlet dry'
+            Q = barrel_velocity = outlet_culvert_depth = 0.0
+            case = 'Inlet dry'
 
 
@@ -160 +164 @@
 # define separately so that can be imported in parallel code.
 #=============================================================================
-def boyd_pipe_function(depth, diameter, length, driving_energy, delta_total_energy, outlet_enquiry_depth, sum_loss, manning):
+def boyd_pipe_function(depth, diameter, blockage, length, driving_energy, delta_total_energy, outlet_enquiry_depth, sum_loss, manning):
 
 
@@ -175 +179 @@
     """
 
+
+    # Note this errors if blockage is set to 1.0 (ie 100% blockaage) and i have no idea how to fix it   
+    if blockage >= 1.0:
+        Q = barrel_velocity = outlet_culvert_depth = 0.0
+        flow_area = 0.00001
+        case = '100 blocked culvert'
+        return Q, barrel_velocity, outlet_culvert_depth, flow_area, case
+    if blockage > 0.9:
+        bf = 3.333-3.333*blockage     
+    else:
+        bf = 1.0-0.4012316798*blockage-0.3768350138*(blockage**2) 
+    #print 'blockage ', blockage
+    #print '      bf ', bf
+
     # Calculate flows for inlet control for circular pipe
-    Q_inlet_unsubmerged = 0.421*anuga.g**0.5*diameter**0.87*driving_energy**1.63 # Inlet Ctrl Inlet Unsubmerged
-    Q_inlet_submerged = 0.530*anuga.g**0.5*diameter**1.87*driving_energy**0.63   # Inlet Ctrl Inlet Submerged
+    Q_inlet_unsubmerged = 0.421*anuga.g**0.5*((bf*diameter)**0.87)*driving_energy**1.63 # Inlet Ctrl Inlet Unsubmerged
+    Q_inlet_submerged = 0.530*anuga.g**0.5*((bf*diameter)**1.87)*driving_energy**0.63   # Inlet Ctrl Inlet Submerged
     # Note for to SUBMERGED TO OCCUR operator.inflow.get_average_specific_energy() should be > 1.2 x diameter.... Should Check !!!
 
@@ -185 +203 @@
     # THE LOWEST Value will Control Calcs From here
     # Calculate Critical Depth Based on the Adopted Flow as an Estimate
-    dcrit1 = diameter/1.26*(Q/anuga.g**0.5*diameter**2.5)**(1/3.75)
-    dcrit2 = diameter/0.95*(Q/anuga.g**0.5*diameter**2.5)**(1/1.95)
+    dcrit1 = (bf*diameter)/1.26*(Q/anuga.g**0.5*((bf*diameter)**2.5))**(1/3.75)
+    dcrit2 = (bf*diameter)/0.95*(Q/anuga.g**0.5*(bf*diameter)**2.5)**(1/1.95)
     # From Boyd Paper ESTIMATE of Dcrit has 2 criteria as
-    if dcrit1/diameter  > 0.85:
+    if dcrit1/(bf*diameter)  > 0.85:
         outlet_culvert_depth = dcrit2
     else:
@@ -194 +212 @@
     #outlet_culvert_depth = min(outlet_culvert_depth, diameter)
     # Now determine Hydraulic Radius Parameters Area & Wetted Perimeter
-    if outlet_culvert_depth >= diameter:
-        outlet_culvert_depth = diameter  # Once again the pipe is flowing full not partfull
-        flow_area = (diameter/2)**2 * math.pi  # Cross sectional area of flow in the culvert
-        perimeter = diameter * math.pi
-        flow_width= diameter
+    if outlet_culvert_depth >= (bf*diameter):
+        outlet_culvert_depth = (bf*diameter)  # Once again the pipe is flowing full not partfull
+        flow_area = (bf*diameter/2)**2 * math.pi  # Cross sectional area of flow in the culvert
+        perimeter = bf * diameter * math.pi
+        flow_width= bf * diameter
         case = 'Inlet CTRL Outlet submerged Circular PIPE FULL'
         if local_debug:
@@ -205 +223 @@
     else:
         #alpha = anuga.acos(1 - outlet_culvert_depth/diameter)    # Where did this Come From ????/
-        alpha = anuga.acos(1-2*outlet_culvert_depth/diameter)*2
+        alpha = anuga.acos(1-2*outlet_culvert_depth/(bf*diameter))*2
         #flow_area = diameter**2 * (alpha - sin(alpha)*cos(alpha))        # Pipe is Running Partly Full at the INLET   WHRE did this Come From ?????
-        flow_area = diameter**2/8*(alpha - math.sin(alpha))   # Equation from  GIECK 5th Ed. Pg. B3
-        flow_width= diameter*math.sin(alpha/2.0)
-        perimeter = alpha*diameter/2.0
+        flow_area = (bf*diameter)**2/8*(alpha - math.sin(alpha))   # Equation from  GIECK 5th Ed. Pg. B3
+        flow_width= bf*diameter*math.sin(alpha/2.0)
+        perimeter = alpha*bf*diameter/2.0
         case = 'INLET CTRL Culvert is open channel flow we will for now assume critical depth'
         if local_debug:
@@ -221 +239 @@
 
         # Determine the depth at the outlet relative to the depth of flow in the Culvert
-        if outlet_enquiry_depth > diameter:       # Outlet is submerged Assume the end of the Pipe is flowing FULL
-            outlet_culvert_depth=diameter
-            flow_area = (diameter/2)**2 * math.pi  # Cross sectional area of flow in the culvert
-            perimeter = diameter * math.pi
-            flow_width= diameter
+        if outlet_enquiry_depth > bf*diameter:       # Outlet is submerged Assume the end of the Pipe is flowing FULL
+            outlet_culvert_depth=bf*diameter
+            flow_area = (bf*diameter/2)**2 * math.pi  # Cross sectional area of flow in the culvert
+            perimeter = bf*diameter * math.pi
+            flow_width= bf*diameter
             case = 'Outlet submerged'
             if local_debug:
@@ -231 +249 @@
         else:   # Culvert running PART FULL for PART OF ITS LENGTH   Here really should use the Culvert Slope to calculate Actual Culvert Depth & Velocity
             # IF  operator.outflow.get_average_depth() < diameter
-            dcrit1 = diameter/1.26*(Q/anuga.g**0.5*diameter**2.5)**(1/3.75)
-            dcrit2 = diameter/0.95*(Q/anuga.g**0.5*diameter**2.5)**(1/1.95)
-            if dcrit1/diameter >0.85:
+            dcrit1 = (bf*diameter)/1.26*(Q/anuga.g**0.5*((bf*diameter)**2.5))**(1/3.75)
+            dcrit2 = (bf*diameter)/0.95*(Q/anuga.g**0.5*((bf*diameter)**2.5))**(1/1.95)
+            if dcrit1/(bf*diameter) >0.85:
                 outlet_culvert_depth= dcrit2
             else:
                 outlet_culvert_depth = dcrit1
-            if outlet_culvert_depth > diameter:
-                outlet_culvert_depth = diameter  # Once again the pipe is flowing full not partfull
-                flow_area = (diameter/2)**2 * math.pi  # Cross sectional area of flow in the culvert
-                perimeter = diameter * math.pi
-                flow_width= diameter
+            if outlet_culvert_depth > bf*diameter:
+                outlet_culvert_depth = bf*diameter  # Once again the pipe is flowing full not partfull
+                flow_area = (bf*diameter/2)**2 * math.pi  # Cross sectional area of flow in the culvert
+                perimeter = bf*diameter * math.pi
+                flow_width= bf*diameter
                 case = 'Outlet unsubmerged PIPE FULL'
                 if local_debug:
                     anuga.log.critical('Outlet unsubmerged PIPE FULL')
             else:
-                alpha = anuga.acos(1-2*outlet_culvert_depth/diameter)*2
-                flow_area = diameter**2/8*(alpha - math.sin(alpha))   # Equation from  GIECK 5th Ed. Pg. B3
-                flow_width= diameter*math.sin(alpha/2.0)
-                perimeter = alpha*diameter/2.0
+                alpha = anuga.acos(1-2*outlet_culvert_depth/(bf*diameter))*2
+                flow_area = (bf*diameter)**2/8*(alpha - math.sin(alpha))   # Equation from  GIECK 5th Ed. Pg. B3
+                flow_width= bf*diameter*math.sin(alpha/2.0)
+                perimeter = alpha*bf*diameter/2.0
                 case = 'Outlet is open channel flow we will for now assume critical depth'
                 if local_debug:
@@ -301 +319 @@
     return Q, barrel_velocity, outlet_culvert_depth, flow_area, case
 
-
-
-

trunk/anuga_core/anuga/structures/tests/test_boyd_box_operator.py

@@ -2656 +2656 @@
         assert numpy.allclose(v, v_expected, rtol=2.0e-2) #outflow velocity
         assert numpy.allclose(d, d_expected, rtol=2.0e-2) #depth at outlet used to calc v 
-        
-       
+
+
+    def test_boyd_13(self):
+        """test_boyd_13
+        
+        This tests the Boyd routine with data obtained from ??? by Petar Milevski    
+        """
+        # FIXME(Ole): This test fails (20 Feb 2009)
+
+        g=9.81
+
+
+        inlet_depth=0.150
+        outlet_depth=0.15
+        inlet_velocity=1.00
+        outlet_velocity=0.5
+        
+        culvert_length=10.0
+        culvert_width=3.6
+        culvert_height=1.20
+        culvert_blockage = 0.50
+                
+        culvert_type='box'
+        manning=0.013
+        sum_loss=1.5
+
+        inlet_specific_energy=inlet_depth + 0.5*inlet_velocity**2/g 
+        culvert_slope=1  # % Downward
+        z_in = 10.0
+        z_out = z_in-culvert_length*culvert_slope/100
+        E_in = z_in+inlet_depth + 0.5*inlet_velocity**2/g
+        E_out = z_out+outlet_depth + 0.5*outlet_velocity**2/g
+        delta_total_energy = E_in-E_out
+        inlet_specific_energy=inlet_depth + 0.5*inlet_velocity**2/g 
+
+
+        Q_expected = 0.28
+        v_expected = 1.15
+        d_expected = 0.13
+
+        if verbose:
+            print 50*'='
+            print 'width ',culvert_width
+            print 'depth ',culvert_height
+            print 'blockage',culvert_blockage
+            print 'flow_width ',culvert_width
+            print 'length ' ,culvert_length
+            print 'driving_energy ',inlet_specific_energy
+            print 'delta_total_energy ',delta_total_energy
+            print 'outlet_enquiry_depth ',outlet_depth
+            print 'sum_loss ',sum_loss
+            print 'manning ',manning
+            
+        Q, v, d, flow_area, case= boyd_box_function(culvert_width, 
+                                                    culvert_height,
+                                                    culvert_blockage, 
+                                                    culvert_width, 
+                                                    culvert_length, 
+                                                    inlet_specific_energy, 
+                                                    delta_total_energy, 
+                                                    outlet_depth, 
+                                                    sum_loss,
+                                                    manning)
+        
+#         Q, v, d = boyd_generalised_culvert_model(inlet_depth,
+#                                                  outlet_depth,
+#                                                  inlet_velocity,
+#                                                  outlet_velocity,
+#                                                  inlet_specific_energy, 
+#                                                  delta_total_energy, 
+#                                                  g,
+#                                                  culvert_length,
+#                                                  culvert_width,
+#                                                  culvert_height,
+#                                                  culvert_type,
+#                                                  manning,
+#                                                  sum_loss)
+        if verbose:
+            print ('%s,%.2f,%.2f,%.2f' %('ANUGAcalcsTEST01 Q-v-d',Q,v,d))
+            print('%s,%.2f,%.2f,%.2f' %('Spreadsheet_Boydcalcs', Q_expected, v_expected, d_expected))
+
+        assert numpy.allclose(Q, Q_expected, rtol=1.0e-1) #inflow
+        assert numpy.allclose(v, v_expected, rtol=1.0e-1) #outflow velocity
+        assert numpy.allclose(d, d_expected, rtol=1.0e-1) #depth at outlet used to calc v 
+
+    def test_boyd_13_operator(self):
+        """test_boyd_non_skew
+        
+        This tests the Boyd routine with data obtained from culvertw application 1.1 by IceMindserer  BD Parkinson, 
+        calculation code by MJ Boyd 
+        """
+
+        g=9.81
+
+
+        inlet_depth=0.150
+        outlet_depth=0.15
+        inlet_velocity=1.00
+        outlet_velocity=0.5
+        
+        culvert_length=10.0
+        culvert_width=3.6
+        culvert_height=1.20
+        
+        culvert_type='box'
+        manning=0.013
+        sum_loss=1.5
+
+        inlet_specific_energy=inlet_depth + 0.5*inlet_velocity**2/g 
+        culvert_slope=1  # % Downward
+        z_in = 10.0
+        z_out = z_in-culvert_length*culvert_slope/100
+        E_in = z_in+inlet_depth + 0.5*inlet_velocity**2/g
+        E_out = z_out+outlet_depth + 0.5*outlet_velocity**2/g
+        delta_total_energy = E_in-E_out
+        inlet_specific_energy=inlet_depth + 0.5*inlet_velocity**2/g 
+
+
+        Q_expected = 0.5526
+        v_expected = 1.146
+        d_expected = 0.1339
+        
+      
+        
+        elevation_0 = z_in
+        elevation_1 = z_out
+        
+        stage_0 = elevation_0 + inlet_depth
+        stage_1 = elevation_1 + outlet_depth
+ 
+
+        domain_length = 200.0
+        domain_width = 200.0
+
+        #culvert_length = 20.0
+        #culvert_width = 3.66
+        #culvert_height = 3.66
+        culvert_losses = {'inlet':0.5, 'outlet':1.0, 'bend':0.0, 'grate':0.0, 'pier': 0.0, 'other': 0.0}
+        culvert_mannings = 0.013
+        
+        culvert_apron = 0.0
+        enquiry_gap = 5.0
+
+
+
+
+        domain = self._create_domain(d_length=domain_length,
+                                     d_width=domain_width,
+                                     dx = 5.0,
+                                     dy = 5.0,
+                                     elevation_0 = elevation_0,
+                                     elevation_1 = elevation_1,
+                                     stage_0 = stage_0,
+                                     stage_1 = stage_1,
+                                     xvelocity_0 = inlet_velocity,
+                                     xvelocity_1 = outlet_velocity)
+ 
+
+        #print 'Defining Structures'
+        
+        ep0 = numpy.array([domain_length/2-culvert_length/2, 100.0])
+        ep1 = numpy.array([domain_length/2+culvert_length/2, 100.0])
+        
+        
+        culvert = Boyd_box_operator(domain,
+                                    losses=culvert_losses,
+                                    width=culvert_width,
+                                    end_points=[ep0, ep1],
+                                    height=culvert_height,
+                                    apron=culvert_apron,
+                                    enquiry_gap=enquiry_gap,
+                                    use_momentum_jet=False,
+                                    use_velocity_head=True,
+                                    manning=culvert_mannings,
+                                    logging=False,
+                                    label='3.6x1.2RCBC',
+                                    verbose=False)
+
+        #culvert.determine_inflow_outflow()
+        
+        ( Q, v, d ) = culvert.discharge_routine()
+        
+        if verbose:
+            print 'test_boyd_operator_1'
+            print 'expected ',Q_expected,v_expected, d_expected
+            print 'calc ',Q,v,d
+        
+
+        assert numpy.allclose(Q, Q_expected, rtol=2.0e-2) #inflow
+        assert numpy.allclose(v, v_expected, rtol=2.0e-2) #outflow velocity
+        assert numpy.allclose(d, d_expected, rtol=2.0e-2) #depth at outlet used to calc v 
+
+
+    def test_boyd_14(self):
+        """test_boyd_14
+        
+        This tests the Boyd routine with data obtained from ??? by Petar Milevski    
+        """
+        # FIXME(Ole): This test fails (20 Feb 2009)
+
+        g=9.81
+        culvert_slope=1  # Downward
+
+        inlet_depth=1.50
+        inlet_velocity= 4.0
+        outlet_depth=0.8
+        outlet_velocity=4.0
+        
+        culvert_length=10.0
+        culvert_width=3.60
+        culvert_height=1.20
+        culvert_blockage = 0.50
+        
+        culvert_type='box'
+        manning=0.013
+        sum_loss=1.5
+
+        inlet_specific_energy=inlet_depth + 0.5*inlet_velocity**2/g 
+        z_in = 10.0
+        z_out = 10.0-culvert_length*culvert_slope/100
+        E_in = z_in+inlet_depth + 0.5*inlet_velocity**2/g
+        E_out = z_out+outlet_depth + 0.5*outlet_velocity**2/g
+        delta_total_energy = E_in-E_out
+
+
+        Q_expected = 3.88
+        v_expected = 2.77
+        d_expected = 0.78
+        
+        if verbose:
+            print 50*'='
+            print 'width ',culvert_width
+            print 'depth ',culvert_height
+            print 'blockage',culvert_blockage
+            print 'flow_width ',culvert_width
+            print 'length ' ,culvert_length
+            print 'driving_energy ',inlet_specific_energy
+            print 'delta_total_energy ',delta_total_energy
+            print 'outlet_enquiry_depth ',outlet_depth
+            print 'sum_loss ',sum_loss
+            print 'manning ',manning
+            
+        Q, v, d, flow_area, case= boyd_box_function(culvert_width, 
+                                                    culvert_height,
+                                                    culvert_blockage, 
+                                                    culvert_width, 
+                                                    culvert_length, 
+                                                    inlet_specific_energy, 
+                                                    delta_total_energy, 
+                                                    outlet_depth, 
+                                                    sum_loss,
+                                                    manning)
+
+#         Q, v, d = boyd_generalised_culvert_model(inlet_depth,
+#                                                  outlet_depth,
+#                                                  inlet_velocity,
+#                                                  outlet_velocity,
+#                                                  inlet_specific_energy, 
+#                                                  delta_total_energy, 
+#                                                  g,
+#                                                  culvert_length,
+#                                                  culvert_width,
+#                                                  culvert_height,
+#                                                  culvert_type,
+#                                                  manning,
+#                                                  sum_loss)        
+        if verbose:
+            print ('%s,%.3f'%('SPEC_E = ',inlet_specific_energy))
+            print ('%s,%.3f'%('Delta E = ',delta_total_energy))
+        
+            print ('%s,%.3f,%.3f,%.3f' %('ANUGAcalcsTEST06 Q-v-d',Q,v,d))
+            print('%s,%.2f,%.2f,%.2f' %('Spreadsheet_Boydcalcs', Q_expected, v_expected, d_expected))
+
+        assert numpy.allclose(Q, Q_expected, rtol=1.0e-1) #inflow
+        assert numpy.allclose(v, v_expected, rtol=1.0e-1) #outflow velocity
+        assert numpy.allclose(d, d_expected, rtol=1.0e-1) #depth at outlet used to calc v 
+
+    def test_boyd_14_operator(self):
+        """test_boyd_non_skew
+        
+        This tests the Boyd routine with data obtained from culvertw application 1.1 by IceMindserer  BD Parkinson, 
+        calculation code by MJ Boyd 
+        """
+
+        g=9.81
+        culvert_slope=1  # Downward
+
+        inlet_depth=1.50
+        inlet_velocity= 4.0
+        outlet_depth=0.8
+        outlet_velocity=4.0
+        culvert_length=10.0
+        culvert_width=3.60
+        culvert_height=1.20
+        
+        culvert_type='box'
+        manning=0.013
+        sum_loss=1.5
+
+        inlet_specific_energy=inlet_depth + 0.5*inlet_velocity**2/g 
+        z_in = 10.0
+        z_out = 10.0-culvert_length*culvert_slope/100
+        E_in = z_in+inlet_depth + 0.5*inlet_velocity**2/g
+        E_out = z_out+outlet_depth + 0.5*outlet_velocity**2/g
+        delta_total_energy = E_in-E_out
+
+
+        Q_expected = 13.546
+        v_expected = 3.136
+        d_expected = 1.20
+        
+        elevation_0 = z_in
+        elevation_1 = z_out
+        
+        stage_0 = elevation_0 + inlet_depth
+        stage_1 = elevation_1 + outlet_depth
+ 
+
+        domain_length = 200.0
+        domain_width = 200.0
+
+        #culvert_length = 20.0
+        #culvert_width = 3.66
+        #culvert_height = 3.66
+        culvert_losses = {'inlet':0.5, 'outlet':1.0, 'bend':0.0, 'grate':0.0, 'pier': 0.0, 'other': 0.0}
+        culvert_mannings = 0.013
+        
+        culvert_apron = 0.0
+        enquiry_gap = 5.0
+
+
+
+
+        domain = self._create_domain(d_length=domain_length,
+                                     d_width=domain_width,
+                                     dx = 5.0,
+                                     dy = 5.0,
+                                     elevation_0 = elevation_0,
+                                     elevation_1 = elevation_1,
+                                     stage_0 = stage_0,
+                                     stage_1 = stage_1,
+                                     xvelocity_0 = inlet_velocity,
+                                     xvelocity_1 = outlet_velocity)
+ 
+
+        #print 'Defining Structures'
+        
+        ep0 = numpy.array([domain_length/2-culvert_length/2, 100.0])
+        ep1 = numpy.array([domain_length/2+culvert_length/2, 100.0])
+        
+        
+        culvert = Boyd_box_operator(domain,
+                                    losses=culvert_losses,
+                                    width=culvert_width,
+                                    end_points=[ep0, ep1],
+                                    height=culvert_height,
+                                    apron=culvert_apron,
+                                    enquiry_gap=enquiry_gap,
+                                    use_momentum_jet=False,
+                                    use_velocity_head=True,
+                                    manning=culvert_mannings,
+                                    logging=False,
+                                    label='3.6x1.2RCBC',
+                                    verbose=False)
+
+        #culvert.determine_inflow_outflow()
+        
+        ( Q, v, d ) = culvert.discharge_routine()
+        
+        if verbose:
+            print 'test_boyd_operator_1'
+            print 'expected ',Q_expected,v_expected, d_expected
+            print 'calc ',Q,v,d
+        
+
+        assert numpy.allclose(Q, Q_expected, rtol=2.0e-2) #inflow
+        assert numpy.allclose(v, v_expected, rtol=2.0e-2) #outflow velocity
+        assert numpy.allclose(d, d_expected, rtol=2.0e-2) #depth at outlet used to calc v 
+        
+
+    def test_boyd_14_operator_invert(self):
+        """test_boyd_non_skew
+        
+        This tests the Boyd routine with data obtained from culvertw application 1.1 by IceMindserer  BD Parkinson, 
+        calculation code by MJ Boyd 
+        """
+
+        #verbose = True
+        g=9.81
+        culvert_slope=1  # Downward
+
+        inlet_depth=1.50
+        inlet_velocity= 4.0
+        outlet_depth=0.8
+        outlet_velocity=4.0
+        culvert_length=10.0
+        culvert_width=3.60
+        culvert_height=1.20
+        
+        culvert_type='box'
+        manning=0.013
+        sum_loss=1.5
+
+        inlet_specific_energy=inlet_depth + 0.5*inlet_velocity**2/g 
+        z_in = 10.0
+        z_out = 10.0-culvert_length*culvert_slope/100
+        E_in = z_in+inlet_depth + 0.5*inlet_velocity**2/g
+        E_out = z_out+outlet_depth + 0.5*outlet_velocity**2/g
+        delta_total_energy = E_in-E_out
+
+        invert0 = 20.0
+        invert1 = 30.0
+
+
+        #Q_expected = 13.546
+        #v_expected = 3.136
+        #d_expected = 1.20
+        
+        Q_expected = 0.0
+        v_expected = 0.0
+        d_expected = 0.0
+        
+        elevation_0 = z_in
+        elevation_1 = z_out
+        
+        stage_0 = elevation_0 + inlet_depth
+        stage_1 = elevation_1 + outlet_depth
+ 
+
+        domain_length = 200.0
+        domain_width = 200.0
+
+        #culvert_length = 20.0
+        #culvert_width = 3.66
+        #culvert_height = 3.66
+        culvert_losses = {'inlet':0.5, 'outlet':1.0, 'bend':0.0, 'grate':0.0, 'pier': 0.0, 'other': 0.0}
+        culvert_mannings = 0.013
+        
+        culvert_apron = 0.0
+        enquiry_gap = 5.0
+
+
+
+
+        domain = self._create_domain(d_length=domain_length,
+                                     d_width=domain_width,
+                                     dx = 5.0,
+                                     dy = 5.0,
+                                     elevation_0 = elevation_0,
+                                     elevation_1 = elevation_1,
+                                     stage_0 = stage_0,
+                                     stage_1 = stage_1,
+                                     xvelocity_0 = inlet_velocity,
+                                     xvelocity_1 = outlet_velocity)
+ 
+
+        #print 'Defining Structures'
+        
+        ep0 = numpy.array([domain_length/2-culvert_length/2, 100.0])
+        ep1 = numpy.array([domain_length/2+culvert_length/2, 100.0])
+        
+        
+        culvert = Boyd_box_operator(domain,
+                                    losses=culvert_losses,
+                                    width=culvert_width,
+                                    end_points=[ep0, ep1],
+                                    height=culvert_height,
+                                    apron=culvert_apron,
+                                    #invert_elevations=[elevation_0,elevation_1],
+                                    invert_elevations=[invert0,invert1],
+                                    enquiry_gap=enquiry_gap,
+                                    use_momentum_jet=False,
+                                    use_velocity_head=True,
+                                    manning=culvert_mannings,
+                                    logging=False,
+                                    label='3.6x1.2RCBC',
+                                    verbose=verbose)
+
+        #culvert.determine_inflow_outflow()
+        
+        ( Q, v, d ) = culvert.discharge_routine()
+        
+        if verbose:
+            print 'test_boyd_operator_12_invert'
+            print 'expected ',Q_expected,v_expected, d_expected
+            print 'calc ',Q,v,d
+        
+
+        assert numpy.allclose(Q, Q_expected, rtol=2.0e-2) #inflow
+        assert numpy.allclose(v, v_expected, rtol=2.0e-2) #outflow velocity
+        assert numpy.allclose(d, d_expected, rtol=2.0e-2) #depth at outlet used to calc v           
 # =========================================================================
 if __name__ == "__main__":

trunk/anuga_core/anuga/structures/tests/test_boyd_pipe_operator.py

@@ -122 +122 @@
         culvert_width = 1.2
         ##culvert_height = 3.66
+        culvert_blockage = 0.0
+        
         culvert_losses = {'inlet':0.5, 'outlet':1.0, 'bend':0.0, 'grate':0.0, 'pier': 0.0, 'other': 0.0}
         culvert_mannings = 0.013
@@ -153 +155 @@
                                     losses=culvert_losses,
                                     diameter=culvert_width,
+                                    blockage=culvert_blockage,
                                     end_points=[ep0, ep1],
                                     #height=culvert_height,
@@ -196 +199 @@
         culvert_width = 1.2
         ##culvert_height = 3.66
+        culvert_blockage = 0.0
+                
         culvert_losses = {'inlet':0.5, 'outlet':1.0, 'bend':0.0, 'grate':0.0, 'pier': 0.0, 'other': 0.0}
         culvert_mannings = 0.013
@@ -227 +232 @@
                                     losses=culvert_losses,
                                     diameter=culvert_width,
+                                    blockage=culvert_blockage,
                                     end_points=[ep0, ep1],
                                     #height=culvert_height,
@@ -270 +276 @@
         culvert_width = 1.2
         ##culvert_height = 3.66
+        culvert_blockage=0.0
+        
         culvert_losses = {'inlet':0.5, 'outlet':1.0, 'bend':0.0, 'grate':0.0, 'pier': 0.0, 'other': 0.0}
         culvert_mannings = 0.013
@@ -301 +309 @@
                                     losses=culvert_losses,
                                     diameter=culvert_width,
+                                    blockage=culvert_blockage,
                                     end_points=[ep0, ep1],
                                     #height=culvert_height,
@@ -345 +354 @@
         culvert_width = 1.2
         ##culvert_height = 3.66
+        culvert_blockage = 0.0
+        
         culvert_losses = {'inlet':0.5, 'outlet':1.0, 'bend':0.0, 'grate':0.0, 'pier': 0.0, 'other': 0.0}
         culvert_mannings = 0.013
@@ -376 +387 @@
                                     losses=culvert_losses,
                                     diameter=culvert_width,
+                                    blockage=culvert_blockage,
                                     end_points=[ep0, ep1],
                                     #height=culvert_height,
@@ -420 +432 @@
         culvert_width = 1.2
         ##culvert_height = 3.66
+        culvert_blockage = 0.0
+        
         culvert_losses = {'inlet':0.5, 'outlet':1.0, 'bend':0.0, 'grate':0.0, 'pier': 0.0, 'other': 0.0}
         culvert_mannings = 0.013
@@ -451 +465 @@
                                     losses=culvert_losses,
                                     diameter=culvert_width,
+                                    blockage=culvert_blockage,
                                     end_points=[ep0, ep1],
                                     #height=culvert_height,
@@ -475 +490 @@
         assert numpy.allclose(Q, expected_Q, rtol=1.0e-2) #inflow
         assert numpy.allclose(v, expected_v, rtol=1.0e-2) #outflow velocity
+        assert numpy.allclose(d, expected_d, rtol=1.0e-2) #depth at outlet used to calc v        
+        
+    def test_boyd_non_skew6(self):
+        """test_boyd_non_skew
+        
+        This tests the Boyd routine with data obtained from culvertw application 1.1 by IceMindserer  BD Parkinson, 
+        calculation code by MJ Boyd
+        This tests the blockage code 
+        """
+
+        stage_0 = 15.0 #change
+        stage_1 = 14.0 #change
+        elevation_0 = 11.0
+        elevation_1 = 10.0
+
+        domain_length = 200.0
+        domain_width = 200.0
+
+        culvert_length = 20.0
+        culvert_width = 1.2
+        ##culvert_height = 3.66
+        culvert_blockage = 0.50
+        
+        culvert_losses = {'inlet':0.5, 'outlet':1.0, 'bend':0.0, 'grate':0.0, 'pier': 0.0, 'other': 0.0}
+        culvert_mannings = 0.013
+        
+        culvert_apron = 0.0
+        enquiry_gap = 5.0
+
+        
+        expected_Q = 1.75
+        expected_v = 3.11
+        expected_d = 0.85
+        
+
+        domain = self._create_domain(d_length=domain_length,
+                                     d_width=domain_width,
+                                     dx = 5.0,
+                                     dy = 5.0,
+                                     elevation_0 = elevation_0,
+                                     elevation_1 = elevation_1,
+                                     stage_0 = stage_0,
+                                     stage_1 = stage_1)
+ 
+
+        #print 'Defining Structures'
+        
+        ep0 = numpy.array([domain_length/2-culvert_length/2, 100.0])
+        ep1 = numpy.array([domain_length/2+culvert_length/2, 100.0])
+        
+        
+        culvert = Boyd_pipe_operator(domain,
+                                    losses=culvert_losses,
+                                    diameter=culvert_width,
+                                    blockage=culvert_blockage,
+                                    end_points=[ep0, ep1],
+                                    #height=culvert_height,
+                                    apron=culvert_apron,
+                                    enquiry_gap=enquiry_gap,
+                                    use_momentum_jet=False,
+                                    use_velocity_head=False,
+                                    manning=culvert_mannings,
+                                    logging=False,
+                                    label='1.2pipe',
+                                    verbose=False)
+
+        #culvert.determine_inflow_outflow()
+        
+        ( Q, v, d ) = culvert.discharge_routine()
+        
+        if verbose:
+            print 'test_boyd_non_skew6'
+            print 'Q: ', Q, 'expected_Q: ', expected_Q
+            print 'v: ', v, 'expected_v: ', expected_v
+            print 'd: ', d, 'expected_d: ', expected_d
+
+
+        assert numpy.allclose(Q, expected_Q, rtol=1.0e-2) #inflow
+        assert numpy.allclose(v, expected_v, rtol=1.0e-2) #outflow velocity
         assert numpy.allclose(d, expected_d, rtol=1.0e-2) #depth at outlet used to calc v  
+
+    def test_boyd_non_skew7(self):
+        """test_boyd_non_skew
+        
+        This tests the Boyd routine with data obtained from culvertw application 1.1 by IceMindserer  BD Parkinson, 
+        calculation code by MJ Boyd
+        This tests the blockage code 
+        """
+
+        stage_0 = 15.0 #change
+        stage_1 = 14.0 #change
+        elevation_0 = 11.0
+        elevation_1 = 10.0
+
+        domain_length = 200.0
+        domain_width = 200.0
+
+        culvert_length = 20.0
+        culvert_width = 1.2
+        ##culvert_height = 3.66
+        culvert_blockage = 1.0
+        
+        culvert_losses = {'inlet':0.5, 'outlet':1.0, 'bend':0.0, 'grate':0.0, 'pier': 0.0, 'other': 0.0}
+        culvert_mannings = 0.013
+        
+        culvert_apron = 0.0
+        enquiry_gap = 5.0
+
+        
+        expected_Q = 0.0
+        expected_v = 0.0
+        expected_d = 0.0
+        
+
+        domain = self._create_domain(d_length=domain_length,
+                                     d_width=domain_width,
+                                     dx = 5.0,
+                                     dy = 5.0,
+                                     elevation_0 = elevation_0,
+                                     elevation_1 = elevation_1,
+                                     stage_0 = stage_0,
+                                     stage_1 = stage_1)
+ 
+
+        #print 'Defining Structures'
+        
+        ep0 = numpy.array([domain_length/2-culvert_length/2, 100.0])
+        ep1 = numpy.array([domain_length/2+culvert_length/2, 100.0])
+        
+        
+        culvert = Boyd_pipe_operator(domain,
+                                    losses=culvert_losses,
+                                    diameter=culvert_width,
+                                    blockage=culvert_blockage,
+                                    end_points=[ep0, ep1],
+                                    #height=culvert_height,
+                                    apron=culvert_apron,
+                                    enquiry_gap=enquiry_gap,
+                                    use_momentum_jet=False,
+                                    use_velocity_head=False,
+                                    manning=culvert_mannings,
+                                    logging=False,
+                                    label='1.2pipe',
+                                    verbose=False)
+
+        #culvert.determine_inflow_outflow()
+        
+        ( Q, v, d ) = culvert.discharge_routine()
+        
+        if verbose:
+            print 'test_boyd_non_skew7'
+            print 'Q: ', Q, 'expected_Q: ', expected_Q
+            print 'v: ', v, 'expected_v: ', expected_v
+            print 'd: ', d, 'expected_d: ', expected_d
+
+
+        assert numpy.allclose(Q, expected_Q, rtol=1.0e-2, atol=1.0e-5) #inflow
+        assert numpy.allclose(v, expected_v, rtol=1.0e-2, atol=1.0e-5) #outflow velocity
+        assert numpy.allclose(d, expected_d, rtol=1.0e-2, atol=1.0e-5) #depth at outlet used to calc v  
+
 # =========================================================================
 if __name__ == "__main__":

trunk/anuga_core/anuga/structures/tests/test_weir_orifice_trapezoid_operator.py

@@ -16 +16 @@
 
 
-verbose = False
+verbose = True
 
 
@@ -206 +206 @@
         culvert_z1=2
         culvert_z2=2
+        culvert_blockage = 0.0
         
         culvert_type='trapezoid'
@@ -228 +229 @@
             print 'width ',culvert_width
             print 'depth ',culvert_height
+            print 'blockage ',culvert_blockage
             print 'flow_width ',culvert_width
             print 'length ' ,culvert_length
@@ -236 +238 @@
             print 'manning ',manning
         
-        Q, v, d, flow_area, case= weir_orifice_trapezoid_function(culvert_width, 
-                                                    culvert_height, 
-                                                    culvert_width,
-                                                    culvert_slope,
-                                                    culvert_z1,
-                                                    culvert_z2, 
-                                                    culvert_length, 
-                                                    inlet_specific_energy, 
-                                                    delta_total_energy, 
-                                                    outlet_depth, 
-                                                    sum_loss,
-                                                    manning)
-        
-
-
+        Q, v, d, flow_area, case= weir_orifice_trapezoid_function(
+                                                    width      = culvert_width, 
+                                                    depth      = culvert_height,
+                                                    blockage   = culvert_blockage, 
+                                                    z1         = culvert_z1,
+                                                    z2         = culvert_z2, 
+                                                    flow_width = culvert_width,
+                                                    length     = culvert_length, 
+                                                    culvert_slope      = culvert_slope,
+                                                    driving_energy     = inlet_specific_energy, 
+                                                    delta_total_energy = delta_total_energy, 
+                                                    outlet_enquiry_depth = outlet_depth, 
+                                                    sum_loss   = sum_loss,
+                                                    manning    = manning)
+        
 
 
@@ -284 +286 @@
         culvert_z1=2
         culvert_z2=2
+        culvert_blockage = 0.0
         
         culvert_type='trapezoid'
@@ -303 +306 @@
         if verbose:
             print 50*'='
-            print 'UNITTEST ',inspect.stack()[0][3] 
+            print 'UNITTEST ',inspect.stack()[0][3]
             print 'width ',culvert_width
             print 'depth ',culvert_height
+            print 'blockage ',culvert_blockage
             print 'flow_width ',culvert_width
             print 'length ' ,culvert_length
@@ -314 +318 @@
             print 'manning ',manning
         
-        Q, v, d, flow_area, case= weir_orifice_trapezoid_function(culvert_width, 
-                                                    culvert_height, 
-                                                    culvert_width,
-                                                    culvert_slope,
-                                                    culvert_z1,
-                                                    culvert_z2, 
-                                                    culvert_length, 
-                                                    inlet_specific_energy, 
-                                                    delta_total_energy, 
-                                                    outlet_depth, 
-                                                    sum_loss,
-                                                    manning)
+        Q, v, d, flow_area, case= weir_orifice_trapezoid_function(
+                                                    width      = culvert_width, 
+                                                    depth      = culvert_height,
+                                                    blockage   = culvert_blockage, 
+                                                    z1         = culvert_z1,
+                                                    z2         = culvert_z2, 
+                                                    flow_width = culvert_width,
+                                                    length     = culvert_length, 
+                                                    culvert_slope      = culvert_slope,
+                                                    driving_energy     = inlet_specific_energy, 
+                                                    delta_total_energy = delta_total_energy, 
+                                                    outlet_enquiry_depth = outlet_depth, 
+                                                    sum_loss   = sum_loss,
+                                                    manning    = manning)
 
 
@@ -360 +366 @@
         culvert_z1=2
         culvert_z2=2
+        culvert_blockage = 0.0
         
         culvert_type='trapezoid'
@@ -382 +389 @@
             print 'width ',culvert_width
             print 'depth ',culvert_height
+            print 'blockage ',culvert_blockage
             print 'flow_width ',culvert_width
             print 'length ' ,culvert_length
@@ -390 +398 @@
             print 'manning ',manning
         
-        Q, v, d, flow_area, case= weir_orifice_trapezoid_function(culvert_width, 
-                                                    culvert_height, 
-                                                    culvert_width,
-                                                    culvert_slope,
-                                                    culvert_z1,
-                                                    culvert_z2, 
-                                                    culvert_length, 
-                                                    inlet_specific_energy, 
-                                                    delta_total_energy, 
-                                                    outlet_depth, 
-                                                    sum_loss,
-                                                    manning)
+        Q, v, d, flow_area, case= weir_orifice_trapezoid_function(
+                                                    width      = culvert_width, 
+                                                    depth      = culvert_height,
+                                                    blockage   = culvert_blockage, 
+                                                    z1         = culvert_z1,
+                                                    z2         = culvert_z2, 
+                                                    flow_width = culvert_width,
+                                                    length     = culvert_length, 
+                                                    culvert_slope      = culvert_slope,
+                                                    driving_energy     = inlet_specific_energy, 
+                                                    delta_total_energy = delta_total_energy, 
+                                                    outlet_enquiry_depth = outlet_depth, 
+                                                    sum_loss   = sum_loss,
+                                                    manning    = manning)
 
 
@@ -436 +446 @@
         culvert_z1=2
         culvert_z2=2
+        culvert_blockage = 0.0
         
         culvert_type='trapezoid'
@@ -458 +469 @@
             print 'width ',culvert_width
             print 'depth ',culvert_height
+            print 'blockage ',culvert_blockage
             print 'flow_width ',culvert_width
             print 'length ' ,culvert_length
@@ -466 +478 @@
             print 'manning ',manning
         
-        Q, v, d, flow_area, case= weir_orifice_trapezoid_function(culvert_width, 
-                                                    culvert_height, 
-                                                    culvert_width,
-                                                    culvert_slope,
-                                                    culvert_z1,
-                                                    culvert_z2, 
-                                                    culvert_length, 
-                                                    inlet_specific_energy, 
-                                                    delta_total_energy, 
-                                                    outlet_depth, 
-                                                    sum_loss,
-                                                    manning)
+        Q, v, d, flow_area, case= weir_orifice_trapezoid_function(
+                                                    width      = culvert_width, 
+                                                    depth      = culvert_height,
+                                                    blockage   = culvert_blockage, 
+                                                    z1         = culvert_z1,
+                                                    z2         = culvert_z2, 
+                                                    flow_width = culvert_width,
+                                                    length     = culvert_length, 
+                                                    culvert_slope      = culvert_slope,
+                                                    driving_energy     = inlet_specific_energy, 
+                                                    delta_total_energy = delta_total_energy, 
+                                                    outlet_enquiry_depth = outlet_depth, 
+                                                    sum_loss   = sum_loss,
+                                                    manning    = manning)
 
 
@@ -513 +527 @@
         culvert_z1=2
         culvert_z2=2
+        culvert_blockage = 0.0
         
         culvert_type='trapezoid'
@@ -535 +550 @@
             print 'width ',culvert_width
             print 'depth ',culvert_height
+            print 'blockage ',culvert_blockage
             print 'flow_width ',culvert_width
             print 'length ' ,culvert_length
@@ -543 +559 @@
             print 'manning ',manning
         
-        Q, v, d, flow_area, case= weir_orifice_trapezoid_function(culvert_width, 
-                                                    culvert_height, 
-                                                    culvert_width,
-                                                    culvert_slope,
-                                                    culvert_z1,
-                                                    culvert_z2, 
-                                                    culvert_length, 
-                                                    inlet_specific_energy, 
-                                                    delta_total_energy, 
-                                                    outlet_depth, 
-                                                    sum_loss,
-                                                    manning)
+        Q, v, d, flow_area, case= weir_orifice_trapezoid_function(
+                                                    width      = culvert_width, 
+                                                    depth      = culvert_height,
+                                                    blockage   = culvert_blockage, 
+                                                    z1         = culvert_z1,
+                                                    z2         = culvert_z2, 
+                                                    flow_width = culvert_width,
+                                                    length     = culvert_length, 
+                                                    culvert_slope      = culvert_slope,
+                                                    driving_energy     = inlet_specific_energy, 
+                                                    delta_total_energy = delta_total_energy, 
+                                                    outlet_enquiry_depth = outlet_depth, 
+                                                    sum_loss   = sum_loss,
+                                                    manning    = manning)
 
 
@@ -589 +607 @@
         culvert_z1=2
         culvert_z2=2
+        culvert_blockage = 0.0
         
         culvert_type='trapezoid'
@@ -611 +630 @@
             print 'width ',culvert_width
             print 'depth ',culvert_height
+            print 'blockage ',culvert_blockage
             print 'flow_width ',culvert_width
             print 'length ' ,culvert_length
@@ -619 +639 @@
             print 'manning ',manning
         
-        Q, v, d, flow_area, case= weir_orifice_trapezoid_function(culvert_width, 
-                                                    culvert_height, 
-                                                    culvert_width,
-                                                    culvert_slope,
-                                                    culvert_z1,
-                                                    culvert_z2, 
-                                                    culvert_length, 
-                                                    inlet_specific_energy, 
-                                                    delta_total_energy, 
-                                                    outlet_depth, 
-                                                    sum_loss,
-                                                    manning)
+        Q, v, d, flow_area, case= weir_orifice_trapezoid_function(
+                                                    width      = culvert_width, 
+                                                    depth      = culvert_height,
+                                                    blockage   = culvert_blockage, 
+                                                    z1         = culvert_z1,
+                                                    z2         = culvert_z2, 
+                                                    flow_width = culvert_width,
+                                                    length     = culvert_length, 
+                                                    culvert_slope      = culvert_slope,
+                                                    driving_energy     = inlet_specific_energy, 
+                                                    delta_total_energy = delta_total_energy, 
+                                                    outlet_enquiry_depth = outlet_depth, 
+                                                    sum_loss   = sum_loss,
+                                                    manning    = manning)
 
 
@@ -644 +666 @@
         assert numpy.allclose(v, v_expected, rtol=1.0e-1) #outflow velocity
         assert numpy.allclose(d, d_expected, rtol=1.0e-1) #depth at outlet used to calc v 
+
 # =========================================================================
 if __name__ == "__main__":

trunk/anuga_core/anuga/structures/weir_orifice_trapezoid_operator.py

@@ -28 +28 @@
                  z1,
                  z2,
+                 blockage=0.0,  #added by DPM 5/10/2016
                  height=None,
                  end_points=None,
@@ -52 +53 @@
                                           width=width,
                                           height=height,
+                                          blockage=blockage, #added by DPM 5/10/2016
                                           #culvert_slope=culvert_slope,
                                           z1=z1,
@@ -82 +84 @@
         self.culvert_width = self.get_culvert_width()
         self.culvert_height = self.get_culvert_height()
+        self.culvert_blockage = self.get_culvert_blockage()
         self.culvert_slope = self.culvert_slope()
         self.culvert_z1 = self.get_culvert_z1()
@@ -161 +164 @@
                 print 'width ',self.culvert_width
                 print 'depth ',self.culvert_height
+                print 'culvert blockage ',self.culvert_blockage #added by DPM 5/10/2016
                 print 'left bank slope ',self.culvert_z1
                 print 'right bank slope ',self.culvert_z2                
@@ -178 +182 @@
                               weir_orifice_trapezoid_function(width =self.culvert_width,
                                                 depth               =self.culvert_height,
+                                                blockage            =self.culvert_blockage, #added by DPM 5/10/2016
                                                 z1                  =self.culvert_z1,
                                                 z2                  =self.culvert_z2, 
@@ -189 +194 @@
                                                 manning             =self.manning)
         else:
-             Q = barrel_velocity = outlet_culvert_depth = 0.0
-             case = 'Inlet dry'
+            Q = barrel_velocity = outlet_culvert_depth = 0.0
+            case = 'Inlet dry'
 
 
@@ -211 +216 @@
 def weir_orifice_trapezoid_function(width, 
                         depth,
+                        blockage, #added by DPM 5/10/2016
                         z1,
                         z2, 
@@ -227 +233 @@
 
     local_debug = False
-    
-    Q_inlet_unsubmerged = 1.7*((2*width+depth*(z1+z2))/2)*driving_energy**1.50 # Flow based on Inlet Ctrl Inlet Unsubmerged (weir flow equation)
-    Q_inlet_submerged = 0.8*anuga.g**0.5*(0.5*depth*(2*width+depth*(z1+z2)))*driving_energy**0.5  # Flow based on Inlet Ctrl Inlet Submerged (orifice equation)
-
-    # FIXME(Ole): Are these functions really for inlet control?
-    if Q_inlet_unsubmerged < Q_inlet_submerged:
-        Q = Q_inlet_unsubmerged
+
+    if blockage >= 1.0:
+        Q = barrel_velocity = outlet_culvert_depth = 0.0
+        flow_area = 0.00001
+        case = '100 blocked culvert'
+        return Q, barrel_velocity, outlet_culvert_depth, flow_area, case
+    else: 
+        b_depth = (1-blockage)*depth
+        b_width = (1-blockage)*width
+          
+        Q_inlet_unsubmerged = 1.7*((2*b_width+b_depth*(z1+z2))/2)*driving_energy**1.50 # Flow based on Inlet Ctrl Inlet Unsubmerged (weir flow equation)
+        Q_inlet_submerged = 0.8*anuga.g**0.5*(0.5*b_depth*(2*b_width+b_depth*(z1+z2)))*driving_energy**0.5  # Flow based on Inlet Ctrl Inlet Submerged (orifice equation)
+
+
+
+        # FIXME(Ole): Are these functions really for inlet control?
+        if Q_inlet_unsubmerged < Q_inlet_submerged:
+            Q = Q_inlet_unsubmerged
+            # Critical Depths Calculation
+            dcrit=0.00001
+            ic=0
+            dyc=0.001
+            while abs(dyc)>0.00001:
+                Tc=b_width+(z1+z2)*dcrit
+                Pc=b_width+((z1**2+1)**0.5+(z2**2+1)**0.5)*dcrit
+                Ac=0.5*dcrit*(b_width+Tc)
+                Rc=Ac/Pc
+                fc=Ac**1.5*Tc**-0.5-Q/(9.81**0.5)
+                ffc=Ac**1.5*-0.5*Tc**-1.5*(z1+z2)+Tc**-0.5*1.5*Ac**0.5*Tc
+                dyc=-fc/ffc
+                dcrit=dcrit+dyc
+                ic=ic+1
+            dcrit = dcrit
+            if dcrit > depth:
+                dcrit = depth
+                flow_area = b_width*dcrit+0.5*(z1+z2)*dcrit**2
+                perimeter= 2.0*b_width+(z1+z2)*dcrit + (dcrit**2+(z1*dcrit)**2)**0.5 + (dcrit**2+(z2*dcrit)**2)**0.5
+            else: # dcrit < depth
+                flow_area = b_width*dcrit+0.5*(z1+z2)*dcrit**2
+                perimeter= b_width + (dcrit**2+(z1*dcrit)**2)**0.5 + (dcrit**2+(z2*dcrit)**2)**0.5
+            outlet_culvert_depth = dcrit
+            case = 'Inlet unsubmerged Box Acts as Weir'
+        else: # Inlet Submerged but check internal culvert flow depth
+            Q = Q_inlet_submerged
+            # Critical Depths Calculation
+            dcrit=0.00001
+            ic=0
+            dyc=0.001
+            while abs(dyc)>0.00001:
+                Tc=b_width+(z1+z2)*dcrit
+                Pc=b_width+((z1**2+1)**0.5+(z2**2+1)**0.5)*dcrit
+                Ac=0.5*dcrit*(b_width+Tc)
+                Rc=Ac/Pc
+                fc=Ac**1.5*Tc**-0.5-Q/(9.81**0.5)
+                ffc=Ac**1.5*-0.5*Tc**-1.5*(z1+z2)+Tc**-0.5*1.5*Ac**0.5*Tc
+                dyc=-fc/ffc
+                dcrit=dcrit+dyc
+                ic=ic+1
+            dcrit = dcrit
+            if dcrit > depth:
+                dcrit = depth
+                flow_area = b_width*dcrit+0.5*(z1+z2)*dcrit**2
+                perimeter= 2.0*b_width+(z1+z2)*dcrit + (dcrit**2+(z1*dcrit)**2)**0.5 + (dcrit**2+(z2*dcrit)**2)**0.5
+            else: # dcrit < depth
+                flow_area = b_width*dcrit+0.5*(z1+z2)*dcrit**2
+                perimeter= b_width + (dcrit**2+(z1*dcrit)**2)**0.5 + (dcrit**2+(z2*dcrit)**2)**0.5
+            outlet_culvert_depth = dcrit
+            case = 'Inlet submerged Box Acts as Orifice'
         # Critical Depths Calculation
         dcrit=0.00001
@@ -239 +306 @@
         dyc=0.001
         while abs(dyc)>0.00001:
-            Tc=width+(z1+z2)*dcrit
-            Pc=width+((z1**2+1)**0.5+(z2**2+1)**0.5)*dcrit
-            Ac=0.5*dcrit*(width+Tc)
+            Tc=b_width+(z1+z2)*dcrit
+            Pc=b_width+((z1**2+1)**0.5+(z2**2+1)**0.5)*dcrit
+            Ac=0.5*dcrit*(b_width+Tc)
             Rc=Ac/Pc
             fc=Ac**1.5*Tc**-0.5-Q/(9.81**0.5)
@@ -249 +316 @@
             ic=ic+1
         dcrit = dcrit
-        if dcrit > depth:
-            dcrit = depth
-            flow_area = width*dcrit+0.5*(z1+z2)*dcrit**2
-            perimeter= 2.0*width+(z1+z2)*dcrit + (dcrit**2+(z1*dcrit)**2)**0.5 + (dcrit**2+(z2*dcrit)**2)**0.5
-        else: # dcrit < depth
-            flow_area = width*dcrit+0.5*(z1+z2)*dcrit**2
-            perimeter= width + (dcrit**2+(z1*dcrit)**2)**0.5 + (dcrit**2+(z2*dcrit)**2)**0.5
+        # May not need this .... check if same is done above
         outlet_culvert_depth = dcrit
-        case = 'Inlet unsubmerged Box Acts as Weir'
-    else: # Inlet Submerged but check internal culvert flow depth
-        Q = Q_inlet_submerged
-        # Critical Depths Calculation
-        dcrit=0.00001
-        ic=0
-        dyc=0.001
-        while abs(dyc)>0.00001:
-            Tc=width+(z1+z2)*dcrit
-            Pc=width+((z1**2+1)**0.5+(z2**2+1)**0.5)*dcrit
-            Ac=0.5*dcrit*(width+Tc)
-            Rc=Ac/Pc
-            fc=Ac**1.5*Tc**-0.5-Q/(9.81**0.5)
-            ffc=Ac**1.5*-0.5*Tc**-1.5*(z1+z2)+Tc**-0.5*1.5*Ac**0.5*Tc
-            dyc=-fc/ffc
-            dcrit=dcrit+dyc
-            ic=ic+1
-        dcrit = dcrit
-        if dcrit > depth:
-            dcrit = depth
-            flow_area = width*dcrit+0.5*(z1+z2)*dcrit**2
-            perimeter= 2.0*width+(z1+z2)*dcrit + (dcrit**2+(z1*dcrit)**2)**0.5 + (dcrit**2+(z2*dcrit)**2)**0.5
-        else: # dcrit < depth
-            flow_area = width*dcrit+0.5*(z1+z2)*dcrit**2
-            perimeter= width + (dcrit**2+(z1*dcrit)**2)**0.5 + (dcrit**2+(z2*dcrit)**2)**0.5
-        outlet_culvert_depth = dcrit
-        case = 'Inlet submerged Box Acts as Orifice'
-    # Critical Depths Calculation
-    dcrit=0.00001
-    ic=0
-    dyc=0.001
-    while abs(dyc)>0.00001:
-        Tc=width+(z1+z2)*dcrit
-        Pc=width+((z1**2+1)**0.5+(z2**2+1)**0.5)*dcrit
-        Ac=0.5*dcrit*(width+Tc)
-        Rc=Ac/Pc
-        fc=Ac**1.5*Tc**-0.5-Q/(9.81**0.5)
-        ffc=Ac**1.5*-0.5*Tc**-1.5*(z1+z2)+Tc**-0.5*1.5*Ac**0.5*Tc
-        dyc=-fc/ffc
-        dcrit=dcrit+dyc
-        ic=ic+1
-    dcrit = dcrit
-    # May not need this .... check if same is done above
-    outlet_culvert_depth = dcrit
-    if outlet_culvert_depth > depth:
-        outlet_culvert_depth = depth  # Once again the pipe is flowing full not partfull
-        flow_area = width*depth+0.5*(z1+z2)*depth**2  # Cross sectional area of flow in the culvert
-        perimeter = 2.0*width+(z1+z2)*depth + (depth**2+(z1*depth)**2)**0.5 + (depth**2+(z2*depth)**2)**0.5
-        case = 'Inlet CTRL Outlet unsubmerged PIPE PART FULL'
-    else:
-        flow_area = width*outlet_culvert_depth+0.5*(z1+z2)*outlet_culvert_depth**2
-        perimeter = 2.0*width+(z1+z2)*outlet_culvert_depth + (outlet_culvert_depth**2+(z1*outlet_culvert_depth)**2)**0.5 + (outlet_culvert_depth**2+(z2*outlet_culvert_depth)**2)**0.5
-        case = 'INLET CTRL Culvert is open channel flow we will for now assume critical depth'
-    # Initial Estimate of Flow for Outlet Control using energy slope
-    #( may need to include Culvert Bed Slope Comparison)
-    hyd_rad = flow_area/perimeter
-    culvert_velocity = math.sqrt(delta_total_energy/((sum_loss/2/anuga.g) \
-                                                          +(manning**2*length)/hyd_rad**1.33333))
-    Q_outlet_tailwater = flow_area * culvert_velocity
-
-
-    if delta_total_energy < driving_energy:
-        # Calculate flows for outlet control
-
-        # Determine the depth at the outlet relative to the depth of flow in the Culvert
-        if outlet_enquiry_depth > depth:        # The Outlet is Submerged
-            outlet_culvert_depth=depth
-            flow_area=width*depth+0.5*(z1+z2)*depth**2  # Cross sectional area of flow in the culvert
-            perimeter=width + (depth**2+(z1*depth)**2)**0.5 + (depth**2+(z2*depth)**2)**0.5
-            case = 'Outlet submerged'
-        else:   # Here we use the Culvert Slope to calculate Actual Culvert Depth & Velocity             
-            # Iterate to determine Normal Depth...........
+        if outlet_culvert_depth > depth:
+            outlet_culvert_depth = depth  # Once again the pipe is flowing full not partfull
+            flow_area = b_width*b_depth+0.5*(z1+z2)*(b_depth)**2  # Cross sectional area of flow in the culvert
+            perimeter = 2.0*b_width+(z1+z2)*b_depth + ((b_depth)**2+(z1*(b_depth))**2)**0.5 + ((b_depth)**2+(z2*(b_depth))**2)**0.5
+            case = 'Inlet CTRL Outlet unsubmerged PIPE PART FULL'
+        else:
+            flow_area = b_width*outlet_culvert_depth+0.5*(z1+z2)*outlet_culvert_depth**2
+            perimeter = 2.0*b_width+(z1+z2)*outlet_culvert_depth + (outlet_culvert_depth**2+(z1*outlet_culvert_depth)**2)**0.5 + (outlet_culvert_depth**2+(z2*outlet_culvert_depth)**2)**0.5
+            case = 'INLET CTRL Culvert is open channel flow we will for now assume critical depth'
+        # Initial Estimate of Flow for Outlet Control using energy slope
+        #( may need to include Culvert Bed Slope Comparison)
+        hyd_rad = flow_area/perimeter
+        culvert_velocity = math.sqrt(delta_total_energy/((sum_loss/2/anuga.g) \
+                                                              +(manning**2*length)/hyd_rad**1.33333))
+        Q_outlet_tailwater = flow_area * culvert_velocity
+    
+    
+        if delta_total_energy < driving_energy:
+            # Calculate flows for outlet control
+    
+            # Determine the depth at the outlet relative to the depth of flow in the Culvert
+            if outlet_enquiry_depth > depth:        # The Outlet is Submerged
+                outlet_culvert_depth=depth
+                flow_area=b_width*b_depth+0.5*(z1+z2)*(b_depth)**2  # Cross sectional area of flow in the culvert
+                perimeter=b_width + ((b_depth)**2+(z1*(b_depth))**2)**0.5 + ((b_depth)**2+(z2*(b_depth))**2)**0.5
+                case = 'Outlet submerged'
+            else:   # Here we use the Culvert Slope to calculate Actual Culvert Depth & Velocity             
+                # Iterate to determine Normal Depth...........
+                
+                Q = min(Q, Q_outlet_tailwater)
+                dnorm=0.00001
+                idn=1
+                dyn=0.001
+                
+                while abs(dyn)>0.0001:
+                    Tn=b_width+(z1+z2)*dnorm
+                    An=0.5*dnorm*(b_width+Tn)
+                    Pn=b_width+2*dnorm*((z1**2+z2**2)**0.5)
+                    Rn=An/Pn
+                    fn=((culvert_slope**0.5*An*Rn**0.67)/manning) - Q
+                    ffn=((culvert_slope**0.5)/manning)*(((Rn**0.67)*Tn)+(Tn/Pn)-(2*dnorm*Rn/Pn))
+                    dnorm=dnorm-fn/ffn        
+                    dyn=-fn/ffn
+                    idn=idn+1
+                outlet_culvert_depth = dnorm
+                if outlet_culvert_depth > depth:
+                    outlet_culvert_depth=depth
+                    flow_area=b_width*b_depth+0.5*(z1+z2)*(b_depth)**2
+                    perimeter=b_width + ((b_depth)**2+(z1*(b_depth))**2)**0.5 + ((b_depth)**2+(z2*(b_depth))**2)**0.5
+                    case = 'Outlet is Flowing Full'
+                else:
+                    flow_area=b_width*outlet_culvert_depth+0.5*(z1+z2)*outlet_culvert_depth**2
+                    perimeter=b_width + (outlet_culvert_depth**2+(z1*outlet_culvert_depth)**2)**0.5 + (outlet_culvert_depth**2+(z2*outlet_culvert_depth)**2)**0.5
+                    case = 'Outlet is open channel flow'
+    
+            hyd_rad = flow_area/perimeter
+    
+            # Final Outlet control velocity using tail water
+            culvert_velocity = math.sqrt(delta_total_energy/((sum_loss/2/anuga.g)\
+                                                              +(manning**2*length)/hyd_rad**1.33333))
+            Q_outlet_tailwater = flow_area * culvert_velocity
+    
+            Q = min(Q, Q_outlet_tailwater)
+        else:
+    
+            pass
+            #FIXME(Ole): What about inlet control?
+    
+        if  flow_area <= 0.0 :
+            culv_froude = 0.0
+        else:
+            culv_froude=math.sqrt(Q**2*flow_width/(anuga.g*flow_area**3))
             
-            Q = min(Q, Q_outlet_tailwater)
-            dnorm=0.00001
-            idn=1
-            dyn=0.001
-            
-            while abs(dyn)>0.0001:
-                Tn=width+(z1+z2)*dnorm
-                An=0.5*dnorm*(width+Tn)
-                Pn=width+2*dnorm*((z1**2+z2**2)**0.5)
-                Rn=An/Pn
-                fn=((culvert_slope**0.5*An*Rn**0.67)/manning) - Q
-                ffn=((culvert_slope**0.5)/manning)*(((Rn**0.67)*Tn)+(Tn/Pn)-(2*dnorm*Rn/Pn))
-                dnorm=dnorm-fn/ffn        
-                dyn=-fn/ffn
-                idn=idn+1
-            outlet_culvert_depth = dnorm
-            if outlet_culvert_depth > depth:
-                outlet_culvert_depth=depth
-                flow_area=width*depth+0.5*(z1+z2)*depth**2
-                perimeter=width + (depth**2+(z1*depth)**2)**0.5 + (depth**2+(z2*depth)**2)**0.5
-                case = 'Outlet is Flowing Full'
-            else:
-                flow_area=width*outlet_culvert_depth+0.5*(z1+z2)*outlet_culvert_depth**2
-                perimeter=width + (outlet_culvert_depth**2+(z1*outlet_culvert_depth)**2)**0.5 + (outlet_culvert_depth**2+(z2*outlet_culvert_depth)**2)**0.5
-                case = 'Outlet is open channel flow'
-
-        hyd_rad = flow_area/perimeter
-
-        # Final Outlet control velocity using tail water
-        culvert_velocity = math.sqrt(delta_total_energy/((sum_loss/2/anuga.g)\
-                                                          +(manning**2*length)/hyd_rad**1.33333))
-        Q_outlet_tailwater = flow_area * culvert_velocity
-
-        Q = min(Q, Q_outlet_tailwater)
-    else:
-
-        pass
-        #FIXME(Ole): What about inlet control?
-
-    if  flow_area <= 0.0 :
-        culv_froude = 0.0
-    else:
-        culv_froude=math.sqrt(Q**2*flow_width/(anuga.g*flow_area**3))
-        
-    if local_debug:
-        anuga.log.critical('FLOW AREA = %s' % str(flow_area))
-        anuga.log.critical('PERIMETER = %s' % str(perimeter))
-        anuga.log.critical('Q final = %s' % str(Q))
-        anuga.log.critical('FROUDE = %s' % str(culv_froude))
-
-    # Determine momentum at the outlet
-    barrel_velocity = Q/(flow_area + anuga.velocity_protection/flow_area)
-
-    # END CODE BLOCK for DEPTH  > Required depth for CULVERT Flow
-
-    return Q, barrel_velocity, outlet_culvert_depth, flow_area, case
+        if local_debug:
+            anuga.log.critical('FLOW AREA = %s' % str(flow_area))
+            anuga.log.critical('PERIMETER = %s' % str(perimeter))
+            anuga.log.critical('Q final = %s' % str(Q))
+            anuga.log.critical('FROUDE = %s' % str(culv_froude))
+    
+        # Determine momentum at the outlet
+        barrel_velocity = Q/(flow_area + anuga.velocity_protection/flow_area)
+    
+        # END CODE BLOCK for DEPTH  > Required depth for CULVERT Flow
+    
+        return Q, barrel_velocity, outlet_culvert_depth, flow_area, case