Changeset 5436

Timestamp:

Jun 25, 2008, 5:57:01 PM (16 years ago)

Author:

ole

Message:

Moved culvert routines into own area

Files:

• anuga_core/source/anuga/config.py (modified) (1 diff)
• anuga_core/source/anuga/culvert_flows/culvert_flow.py (deleted)
• anuga_core/source/anuga/culvert_flows/culvert_routines.py (modified) (1 diff)
• anuga_core/source/anuga/shallow_water/shallow_water_domain.py (modified) (8 diffs)
• anuga_core/source/anuga/shallow_water/test_shallow_water_domain.py (modified) (1 diff)
• anuga_core/source/anuga/utilities/system_tools.py (modified) (1 diff)
• anuga_work/development/Rudy_vandrie_2007/momentum_jet.py (added)
• anuga_work/development/culvert_flow/culvert_boyd_channel.py (modified) (7 diffs)

anuga_core/source/anuga/config.py

@@ -9 +9 @@
 default_smoothing_parameter = 0.001 # Default alpha for penalised
                                     # least squares fitting
+
+velocity_protection = 1.0e-6                                     
 
 #-------------------------------------------

anuga_core/source/anuga/culvert_flows/culvert_routines.py

@@ +1 @@
+"""Collection of culvert routines for use with Culvert_flow in culvert_class
+
+Usage:
+
+   
+
+"""
+
+#   NEW DEFINED CULVERT FLOW---- Flow from INLET 1 ------> INLET 2 (Outlet)
+# 
+# The First Attempt has a Simple Horizontal Circle as a Hole on the Bed
+# Flow Is Removed at a Rate of INFLOW 
+# Downstream there is a similar Circular Hole on the Bed where INFLOW effectively Surcharges
+#
+# This SHould be changed to a Vertical Opening Both BOX and Circular
+# There will be several Culvert Routines such as:
+# CULVERT_Boyd_Channel
+# CULVERT_Orifice_and_Weir
+# CULVERT_Simple_FLOOR
+# CULVERT_Simple_WALL
+# CULVERT_Eqn_Floor
+# CULVERT_Eqn_Wall
+# CULVERT_Tab_Floor
+# CULVERT_Tab_Wall
+# BRIDGES.....
+# NOTE NEED TO DEVELOP CONCEPT 1D Model for Linked Pipe System !!!!
+
+#  COULD USE EPA SWMM Model !!!!
+
+
+from math import pi, sqrt, sin, cos
+
+
+def boyd_generalised_culvert_model(culvert, inlet, outlet, delta_Et, g):
+
+    """Boyd's generalisation of the US department of transportation culvert model
+        # == The quantity of flow passing through a culvert is controlled by many factors
+        # == It could be that the culvert is controled by the inlet only, with it being Un submerged this is effectively equivalent to the WEIR Equation
+        # == Else the culvert could be controlled by the inlet, with it being Submerged, this is effectively the Orifice Equation
+        # == Else it may be controlled by Down stream conditions where depending on the down stream depth, the momentum in the culvert etc. flow is controlled 
+    """
+
+    from anuga.utilities.system_tools import log_to_file
+    from anuga.config import velocity_protection
+    from anuga.utilities.numerical_tools import safe_acos as acos
+
+        
+    Q_outlet_tailwater = 0.0
+    inlet.rate = 0.0
+    outlet.rate = 0.0
+    Q_inlet_unsubmerged = 0.0
+    Q_inlet_submerged = 0.0
+    Q_outlet_critical_depth = 0.0
+
+    log_filename = culvert.log_filename
+
+    manning = culvert.manning
+    sum_loss = culvert.sum_loss
+    length = culvert.length
+    
+    if inlet.depth >= 0.01:
+        # Calculate driving energy
+        E = inlet.specific_energy
+
+        s = 'Driving energy  = %f m' %E
+        log_to_file(log_filename, s)
+        
+        msg = 'Driving energy is negative'
+        assert E >= 0.0, msg
+                      
+        
+        # Water has risen above inlet
+        if culvert.culvert_type == 'circle':
+            # Round culvert
+
+            # Calculate flows for inlet control            
+            diameter = culvert.diameter
+
+            Q_inlet_unsubmerged = 0.421*g**0.5*diameter**0.87*E**1.63    # Inlet Ctrl Inlet Unsubmerged 
+            Q_inlet_submerged = 0.530*g**0.5*diameter**1.87*E**0.63      # Inlet Ctrl Inlet Submerged
+
+            s = 'Q_inlet_unsubmerged = %.6f, Q_inlet_submerged = %.6f' %(Q_inlet_unsubmerged, Q_inlet_submerged)
+            log_to_file(log_filename, s)
+
+            case = ''
+            if Q_inlet_unsubmerged < Q_inlet_submerged:
+                Q = Q_inlet_unsubmerged
+
+                alpha = acos(1 - inlet.depth/diameter)
+                flow_area = diameter**2 * (alpha - sin(alpha)*cos(alpha))
+                outlet_culvert_depth = inlet.depth
+                width = diameter*sin(alpha)
+                perimeter = alpha*diameter                
+                case = 'Inlet unsubmerged'
+            else:    
+                Q = Q_inlet_submerged
+                flow_area = (diameter/2)**2 * pi
+                outlet_culvert_depth = diameter
+                width = diameter
+                perimeter = diameter
+                case = 'Inlet submerged'                    
+
+
+
+            if delta_Et < E:
+                # Calculate flows for outlet control
+                # Determine the depth at the outlet relative to the depth of flow in the Culvert
+
+                if outlet.depth > diameter:        # The Outlet is Submerged
+                    outlet_culvert_depth=diameter
+                    flow_area = (diameter/2)**2 * pi  # Cross sectional area of flow in the culvert
+                    perimeter = diameter * pi
+                    width = diameter
+                    case = 'Outlet submerged'
+                elif outlet.depth==0.0: 
+                    outlet_culvert_depth=inlet.depth   # For purpose of calculation assume the outlet depth = the inlet depth
+                    alpha = acos(1 - inlet.depth/diameter)
+                    flow_area = diameter**2 * (alpha - sin(alpha)*cos(alpha))
+                    perimeter = alpha*diameter
+                    width = diameter*sin(alpha)                    
+
+                    case = 'Outlet depth is zero'                        
+                else:   # Here really should use the Culvert Slope to calculate Actual Culvert Depth & Velocity
+                    outlet_culvert_depth=outlet.depth   # For purpose of calculation assume the outlet depth = the inlet depth
+                    alpha = acos(1 - outlet.depth/diameter)
+                    flow_area = diameter**2 * (alpha - sin(alpha)*cos(alpha))
+                    perimeter = alpha*diameter
+                    width = diameter*sin(alpha)                    
+                    case = 'Outlet is open channel flow'
+
+
+        else:
+            # Box culvert (rectangle or square)
+
+            # Calculate flows for inlet control
+            height = culvert.height
+            width = culvert.width            
+            
+            Q_inlet_unsubmerged = 0.540*g**0.5*width*E**1.50 # Flow based on Inlet Ctrl Inlet Unsubmerged
+            Q_inlet_submerged = 0.702*g**0.5*width*height**0.89*E**0.61  # Flow based on Inlet Ctrl Inlet Submerged
+
+            s = 'Q_inlet_unsubmerged = %.6f, Q_inlet_submerged = %.6f' %(Q_inlet_unsubmerged, Q_inlet_submerged)
+            log_to_file(log_filename, s)
+
+            case = ''
+            if Q_inlet_unsubmerged < Q_inlet_submerged:
+                Q = Q_inlet_unsubmerged
+                flow_area = width*inlet.depth
+                outlet_culvert_depth = inlet.depth
+                perimeter=(width+2.0*inlet.depth)                
+                case = 'Inlet unsubmerged'
+            else:    
+                Q = Q_inlet_submerged
+                flow_area = width*height
+                outlet_culvert_depth = height
+                perimeter=2.0*(width+height)                
+                case = 'Inlet submerged'                    
+
+            if delta_Et < E:
+                # Calculate flows for outlet control
+                # Determine the depth at the outlet relative to the depth of flow in the Culvert
+
+                if outlet.depth > height:        # The Outlet is Submerged
+                    outlet_culvert_depth=height
+                    flow_area=width*height       # Cross sectional area of flow in the culvert
+                    perimeter=2.0*(width+height)
+                    case = 'Outlet submerged'
+                elif outlet.depth==0.0: 
+                    outlet_culvert_depth=inlet.depth   # For purpose of calculation assume the outlet depth = the inlet depth
+                    flow_area=width*inlet.depth
+                    perimeter=(width+2.0*inlet.depth)
+                    case = 'Outlet depth is zero'                        
+                else:   # Here really should use the Culvert Slope to calculate Actual Culvert Depth & Velocity
+                    outlet_culvert_depth=outlet.depth
+                    flow_area=width*outlet.depth
+                    perimeter=(width+2.0*outlet.depth)
+                    case = 'Outlet is open channel flow'
+                    
+
+
+
+        # Common code for rectangular and circular culvert types    
+        hyd_rad = flow_area/perimeter
+        s = 'hydraulic radius at outlet = %f' %hyd_rad
+        log_to_file(log_filename, s)
+
+        # Outlet control velocity using tail water
+        culvert_velocity = sqrt(delta_Et/((sum_loss/2*g)+(manning**2*length)/hyd_rad**1.33333)) 
+        Q_outlet_tailwater = flow_area * culvert_velocity
+
+        s = 'Q_outlet_tailwater = %.6f' %Q_outlet_tailwater
+        log_to_file(log_filename, s)
+        Q = min(Q, Q_outlet_tailwater)
+
+        log_to_file(log_filename, 'Case: "%s"' %case)
+    
+        flow_rate_control=Q
+
+        s = 'Flow Rate Control = %f' %flow_rate_control
+        log_to_file(log_filename, s)
+
+        inlet.rate = -flow_rate_control
+        outlet.rate = flow_rate_control                
+            
+        culv_froude=sqrt(flow_rate_control**2*width/(g*flow_area**3))
+        s = 'Froude in Culvert = %f' %culv_froude
+        log_to_file(log_filename, s)
+
+        # Determine momentum at the outlet 
+        barrel_velocity = Q/(flow_area + velocity_protection/flow_area)
+    else: #inlet.depth < 0.01:
+        Q = barrel_velocity = outlet_culvert_depth = 0.0
+        
+    return Q, barrel_velocity, outlet_culvert_depth 
+
+

anuga_core/source/anuga/shallow_water/shallow_water_domain.py

@@ -1527 +1527 @@
                 negative values indicate decreases.
                 Rate can be None at initialisation but must be specified
-                before forcting term is applied (i.e. simulation has started).
+                before forcing term is applied (i.e. simulation has started).
 
     center [m]: Coordinates at center of flow point
@@ -1565 +1565 @@
 
         # Update area if applicable
-        self.area = None        
+        self.exchange_area = None        
         if center is not None and radius is not None:
             assert len(center) == 2
@@ -1571 +1571 @@
             assert polygon is None, msg
 
-            self.area = radius**2*pi
+            self.exchange_area = radius**2*pi
         
         if polygon is not None:
-            self.area = polygon_area(self.polygon)
+            self.exchange_area = polygon_area(self.polygon)
 
 
@@ -1584 +1584 @@
         points = domain.get_centroid_coordinates(absolute=True)
 
-        self.indices = None
+        # Calculate indices in exchange area for this forcing term
+        self.exchange_indices = None
         if self.center is not None and self.radius is not None:
             # Inlet is circular
             
-            self.indices = []
+            self.exchange_indices = []
             for k in range(N):
                 x, y = points[k,:] # Centroid
                 if ((x-self.center[0])**2+(y-self.center[1])**2) < self.radius**2:
-                    self.indices.append(k)
+                    self.exchange_indices.append(k)
                     
         if self.polygon is not None:                    
             # Inlet is polygon
-            self.indices = inside_polygon(points, self.polygon)
+            self.exchange_indices = inside_polygon(points, self.polygon)
             
 
@@ -1621 +1622 @@
 
 
-        if self.indices is None:
+        if self.exchange_indices is None:
             self.update[:] += rate
         else:
             # Brute force assignment of restricted rate
-            for k in self.indices:
+            for k in self.exchange_indices:
                 self.update[k] += rate
 
@@ -1645 +1646 @@
         """Return values for specified quantity restricted to opening 
         """
-        return self.domain.quantities[self.quantity_name].get_values(indices=self.indices)
+        return self.domain.quantities[self.quantity_name].get_values(indices=self.exchange_indices)
     
 
@@ -1651 +1652 @@
         """Set values for specified quantity restricted to opening 
         """
-        self.domain.quantities[self.quantity_name].set_values(val, indices=self.indices)    
+        self.domain.quantities[self.quantity_name].set_values(val, indices=self.exchange_indices)    
 
 
@@ -1811 +1812 @@
         
         if callable(self.rate):
-            _rate = self.rate(t)/self.area
+            _rate = self.rate(t)/self.exchange_area
         else:
-            _rate = self.rate/self.area
+            _rate = self.rate/self.exchange_area
 
         return _rate

anuga_core/source/anuga/shallow_water/test_shallow_water_domain.py

@@ -2231 +2231 @@
         # Setup only one forcing term, constant rainfall restricted to a polygon enclosing triangle #1 (bce)
         domain.forcing_terms = []
-        domain.forcing_terms.append( Rainfall(domain, rate=2.0, polygon = [[1,1], [2,1], [2,2], [1,2]]))
+        R = Rainfall(domain, rate=2.0, polygon = [[1,1], [2,1], [2,2], [1,2]])
+
+        assert allclose(R.exchange_area, 1)
+        
+        domain.forcing_terms.append(R)
 
         domain.compute_forcing_terms()

anuga_core/source/anuga/utilities/system_tools.py

@@ -6 +6 @@
 import sys
 import os
+
+def log_to_file(filename, s, verbose=True):
+    """Log string to open file descriptor
+    """
+
+    fid = open(filename, 'a')
+    if verbose: print s
+    fid.write(s + '\n')
+    fid.close()
+
 
 def get_user_name():

anuga_work/development/culvert_flow/culvert_boyd_channel.py

@@ -51 +51 @@
 from anuga.shallow_water.shallow_water_domain import Dirichlet_boundary
 from anuga.shallow_water.shallow_water_domain import Inflow, General_forcing
-from anuga.culvert_flows.culvert_polygons import create_culvert_polygons
+
+from anuga.culvert_flows.culvert_class import Culvert_flow
+from anuga.culvert_flows.culvert_routines import boyd_generalised_culvert_model
+
+
 from anuga.utilities.polygon import plot_polygons
-from anuga.utilities.system_tools import log_to_file
-
 
 
@@ -62 +64 @@
 # Setup computational domain
 #------------------------------------------------------------------------------
-
-
-# Open file for storing some specific results...
-fid = open('Culvert_Headwall_VarM.txt', 'w')
 
 length = 40.
@@ -86 +84 @@
 domain.set_minimum_storable_height(0.001)
 
-s='Size'+str(len(domain))
-log_to_file(fid, s)
-
-
-
-velocity_protection = 1.0e-5
-
-
-
-
-
-
-def boyd_generalised_culvert_model(inlet, outlet, delta_Et, width, height, length, sum_loss, manning, g, fid):
-    """Boyd's generalisation of the US department of transportation culvert model
-        # == The quantity of flow passing through a culvert is controlled by many factors
-        # == It could be that the culvert is controled by the inlet only, with it being Un submerged this is effectively equivalent to the WEIR Equation
-        # == Else the culvert could be controlled by the inlet, with it being Submerged, this is effectively the Orifice Equation
-        # == Else it may be controlled by Down stream conditions where depending on the down stream depth, the momentum in the culvert etc. flow is controlled 
-    """
-
-    from anuga.utilities.system_tools import log_to_file
-
-        
-    Q_outlet_tailwater = 0.0
-    inlet.rate = 0.0
-    outlet.rate = 0.0
-    Q_inlet_unsubmerged = 0.0
-    Q_inlet_submerged = 0.0
-    Q_outlet_critical_depth = 0.0
-        
-    if inlet.depth >= 0.01: 
-        # Water has risen above inlet
-        if width==height:  # Need something better to Flag Diamater !!!!!!! Can't we just have DIAMETER as well as Width & Height ???
-            pass
-            #Q1[t]= 0.421*g**0.5*Diam**0.87*HW**1.63    # Inlet Ctrl Inlet Unsubmerged
-            #Q2[t]= 0.530*g**0.5*Diam**1.87*HW**0.63    # Inlet Ctrl Inlet Submerged
-            #PipeDcrit=
-            #Delta_E=HW-TW
-        else:
-            # Box culvert
-            
-            # Calculate flows for inlet control
-            E = inlet.specific_energy
-
-            s = 'Driving energy  = %f m' %E
-            log_to_file(fid, s)
-
-            msg = 'Driving energy is negative'
-            assert E >= 0.0, msg
-                      
-            Q_inlet_unsubmerged = 0.540*g**0.5*width*E**1.50 # Flow based on Inlet Ctrl Inlet Unsubmerged
-            Q_inlet_submerged = 0.702*g**0.5*width*height**0.89*E**0.61  # Flow based on Inlet Ctrl Inlet Submerged
-
-            s = 'Q_inlet_unsubmerged = %.6f, Q_inlet_submerged = %.6f' %(Q_inlet_unsubmerged, Q_inlet_submerged)
-            log_to_file(fid, s)
-
-            case = ''
-            if Q_inlet_unsubmerged < Q_inlet_submerged:
-                Q = Q_inlet_unsubmerged
-                flow_area = width*inlet.depth
-                outlet_culvert_depth = inlet.depth
-                case = 'Inlet unsubmerged'
-            else:    
-                Q = Q_inlet_submerged
-                flow_area = width*height
-                outlet_culvert_depth = height
-                case = 'Inlet submerged'                    
-
-            if delta_Et < E:
-                # Calculate flows for outlet control
-                # Determine the depth at the outlet relative to the depth of flow in the Culvert
-
-                if outlet.depth > height:        # The Outlet is Submerged
-                    outlet_culvert_depth=height
-                    flow_area=width*height       # Cross sectional area of flow in the culvert
-                    perimeter=2.0*(width+height)
-                    case = 'Outlet submerged'
-                elif outlet.depth==0.0: 
-                    outlet_culvert_depth=inlet.depth   # For purpose of calculation assume the outlet depth = the inlet depth
-                    flow_area=width*inlet.depth
-                    perimeter=(width+2.0*inlet.depth)
-                    case = 'Outlet depth is zero'                        
-                else:   # Here really should use the Culvert Slope to calculate Actual Culvert Depth & Velocity
-                    outlet_culvert_depth=outlet.depth
-                    flow_area=width*outlet.depth
-                    perimeter=(width+2.0*outlet.depth)
-                    case = 'Outlet is open channel flow'
-                    
-                hyd_rad = flow_area/perimeter
-                s = 'hydraulic radius at outlet = %f' %hyd_rad
-                log_to_file(fid, s)
-                
-
-                # Outlet control velocity using tail water
-                culvert_velocity = sqrt(delta_Et/((sum_loss/2*g)+(manning**2*length)/hyd_rad**1.33333)) 
-                Q_outlet_tailwater = flow_area * culvert_velocity
-
-                s = 'Q_outlet_tailwater = %.6f' %Q_outlet_tailwater
-                log_to_file(fid, s)
-
-                Q = min(Q, Q_outlet_tailwater)
-
-
-                
-            log_to_file(fid, 'Case: "%s"' %case)
-            
-            flow_rate_control=Q
-
-            s = 'Flow Rate Control = %f' %flow_rate_control
-            log_to_file(fid, s)
-
-            inlet.rate = -flow_rate_control
-            outlet.rate = flow_rate_control                
-            
-            culv_froude=sqrt(flow_rate_control**2*width/(g*flow_area**3))
-            s = 'Froude in Culvert = %f' %culv_froude
-            log_to_file(fid, s)
-
-            # Determine momentum at the outlet 
-            barrel_velocity = Q/(flow_area + velocity_protection/flow_area)
-
-            return Q, barrel_velocity, outlet_culvert_depth 
 
 
@@ -264 +140 @@
 
 
-
-#------------------------------------------------------------------------------
-# Setup specialised forcing terms
-#------------------------------------------------------------------------------
-
- #   NEW DEFINED CULVERT FLOW---- Flow from INLET 1 ------> INLET 2 (Outlet)
- # 
- # The First Attempt has a Simple Horizontal Circle as a Hole on the Bed
- # Flow Is Removed at a Rate of INFLOW 
- # Downstream there is a similar Circular Hole on the Bed where INFLOW effectively Surcharges
- #
- # This SHould be changed to a Vertical Opening Both BOX and Circular
- # There will be several Culvert Routines such as:
- # CULVERT_Boyd_Channel
- # CULVERT_Orifice_and_Weir
- # CULVERT_Simple_FLOOR
- # CULVERT_Simple_WALL
- # CULVERT_Eqn_Floor
- # CULVERT_Eqn_Wall
- # CULVERT_Tab_Floor
- # CULVERT_Tab_Wall
- # BRIDGES.....
- # NOTE NEED TO DEVELOP CONCEPT 1D Model for Linked Pipe System !!!!
- 
- #  COULD USE EPA SWMM Model !!!!
- 
- 
- 
-class Culvert_flow:
-    """Culvert flow - transfer water from one hole to another
-    
-    Using Momentum as Calculated by Culvert Flow !!
-    Could be Several Methods Investigated to do This !!!
-
-    2008_May_08
-    To Ole:
-    OK so here we need to get the Polygon Creating code to create polygons for the culvert Based on
-    the two input Points (X0,Y0) and (X1,Y1) - need to be passed to create polygon
-
-    The two centers are now passed on to create_polygon.
-    
-
-    Input: Two points, pipe_size (either diameter or width, height), mannings_rougness,
-    inlet/outlet energy_loss_coefficients, internal_bend_coefficent,
-    top-down_blockage_factor and bottom_up_blockage_factor
-    
-    
-    And the Delta H enquiry should be change from Openings in line 412 to the enquiry Polygons infront
-    of the culvert
-    At the moment this script uses only Depth, later we can change it to Energy...
-
-    Once we have Delta H can calculate a Flow Rate and from Flow Rate an Outlet Velocity
-    The Outlet Velocity x Outlet Depth = Momentum to be applied at the Outlet...
-        
-    """ 
-
-    def __init__(self,
-                 domain,
-                 label=None, 
-                 description=None,
-                 end_point0=None, 
-                 end_point1=None,
-                 width=None,
-                 height=None,
-                 manning=None,          # Mannings Roughness for Culvert
-                 invert_level0=None,    # Invert level if not the same as the Elevation on the Domain
-                 invert_level1=None,    # Invert level if not the same as the Elevation on the Domain
-                 loss_exit=None,
-                 loss_entry=None,
-                 loss_bend=None,
-                 loss_special=None,
-                 blockage_topdwn=None,
-                 blockage_bottup=None,
-                 culvert_routine=None,
-                 verbose=False):
-        
-        from Numeric import sqrt, sum
-
-        # Input check
-        if height is None: height = width
-
-        # Set defaults
-        if manning is None: manning = 0.012   # Set a Default Mannings Roughness for Pipe
-        if loss_exit is None: loss_exit = 1.00
-        if loss_entry is None: loss_entry = 0.50
-        if loss_bend is None: loss_bend=0.00
-        if loss_special is None: loss_special=0.00
-        if blockage_topdwn is None: blockage_topdwn=0.00
-        if blockage_bottup is None: blockage_bottup=0.00
-        if culvert_routine is None: culvert_routine=boyd_generalised_culvert_model        
-        
-
-        # Create the fundamental culvert polygons from POLYGON 
-        P = create_culvert_polygons(end_point0,
-                                    end_point1,
-                                    width=width,   
-                                    height=height)
-        
-        if verbose is True:
-            pass
-            #plot_polygons([[end_point0, end_point1],
-            #               P['exchange_polygon0'],
-            #               P['exchange_polygon1'],
-            #               P['enquiry_polygon0'],
-            #               P['enquiry_polygon1']],
-            #              figname='culvert_polygon_output')
-        
-        self.openings = []
-        self.openings.append(Inflow(domain,
-                                    polygon=P['exchange_polygon0']))
-
-        self.openings.append(Inflow(domain,
-                                    polygon=P['exchange_polygon1']))                                    
-
-
-        # Assume two openings for now: Referred to as 0 and 1
-        assert len(self.openings) == 2
-        
-        # Store basic geometry 
-        self.end_points = [end_point0, end_point1]
-        self.invert_levels = [invert_level0, invert_level1]                
-        self.enquiry_polygons = [P['enquiry_polygon0'], P['enquiry_polygon1']]
-        self.vector = P['vector']
-        self.distance = P['length']; assert self.distance > 0.0
-        self.verbose = verbose
-        self.width = width
-        self.height = height
-        self.last_time = 0.0        
-        self.temp_keep_delta_t = 0.0
-        
-
-        # Store hydraulic parameters 
-        self.manning = manning
-        self.loss_exit = loss_exit
-        self.loss_entry = loss_entry
-        self.loss_bend = loss_bend
-        self.loss_special = loss_special
-        self.sum_loss = loss_exit + loss_entry + loss_bend + loss_special
-        self.blockage_topdwn = blockage_topdwn
-        self.blockage_bottup = blockage_bottup
-
-        # Store culvert routine
-        self.culvert_routine = culvert_routine
-
-        
-        # Create objects to update momentum (a bit crude at this stage)
-
-        
-        xmom0 = General_forcing(domain, 'xmomentum',
-                                polygon=P['exchange_polygon0'])
-
-        xmom1 = General_forcing(domain, 'xmomentum',
-                                polygon=P['exchange_polygon1'])
-
-        ymom0 = General_forcing(domain, 'ymomentum',
-                                polygon=P['exchange_polygon0'])
-
-        ymom1 = General_forcing(domain, 'ymomentum',
-                                polygon=P['exchange_polygon1'])
-
-        self.opening_momentum = [ [xmom0, ymom0], [xmom1, ymom1] ]
-        
-
-        # Print something
-        s = 'Culvert Effective Length = %.2f m' %(self.distance)
-        log_to_file(fid, s)
-   
-        s = 'Culvert Direction is %s\n' %str(self.vector)
-        log_to_file(fid, s)        
-        
-    def __call__(self, domain):
-        from anuga.utilities.numerical_tools import mean
-        from anuga.utilities.polygon import inside_polygon
-        from anuga.config import g, epsilon
-        from Numeric import take, sqrt
-
-         
-        # Time stuff
-        time = domain.get_time()
-        delta_t = time-self.last_time
-        s = '\nTime = %.2f, delta_t = %f' %(time, delta_t)
-        log_to_file(fid, s)
-        
-        msg = 'Time did not advance'
-        if time > 0.0: assert delta_t > 0.0, msg
-
-
-        # Get average water depths at each opening        
-        openings = self.openings   # There are two Opening [0] and [1]
-        for i, opening in enumerate(openings):
-            stage = domain.quantities['stage'].get_values(location='centroids',
-                                                          indices=opening.exchange_indices)
-            elevation = domain.quantities['elevation'].get_values(location='centroids',
-                                                                  indices=opening.exchange_indices)
-
-            # Indices corresponding to energy enquiry field for this opening
-            coordinates = domain.get_centroid_coordinates() # Get all centroid points (x,y)
-            enquiry_indices = inside_polygon(coordinates, self.enquiry_polygons[i]) 
-            
-            
-            # Get model values for points in enquiry polygon for this opening
-            dq = domain.quantities
-            stage = dq['stage'].get_values(location='centroids', indices=enquiry_indices)
-            xmomentum = dq['xmomentum'].get_values(location='centroids', indices=enquiry_indices)
-            ymomentum = dq['ymomentum'].get_values(location='centroids', indices=enquiry_indices)                        
-            elevation = dq['elevation'].get_values(location='centroids', indices=enquiry_indices)
-            depth = stage - elevation
-
-            # Compute mean values of selected quantitites in the enquiry area in front of the culvert
-            # Epsilon handles a dry cell case
-            ux = xmomentum/(depth+velocity_protection/depth)   # Velocity (x-direction)
-            uy = ymomentum/(depth+velocity_protection/depth)   # Velocity (y-direction)
-            v = mean(sqrt(ux**2+uy**2))      # Average velocity
-            w = mean(stage)                  # Average stage
-
-            # Store values at enquiry field
-            opening.velocity = v
-
-
-            # Compute mean values of selected quantitites in the exchange area in front of the culvert
-            # Stage and velocity comes from enquiry area and elevation from exchange area
-            
-            # Use invert level instead of elevation if specified
-            invert_level = self.invert_levels[i]
-            if invert_level is not None:
-                z = invert_level
-            else:
-                elevation = dq['elevation'].get_values(location='centroids', indices=opening.exchange_indices)
-                z = mean(elevation)                   # Average elevation
-                
-            # Estimated depth above the culvert inlet
-            d = w - z
-
-            if d < 0.0:
-                # This is possible since w and z are taken at different locations
-                #msg = 'D < 0.0: %f' %d
-                #raise msg
-                d = 0.0
-            
-            # Ratio of depth to culvert height.
-            # If ratio > 1 then culvert is running full
-            ratio = d/self.height  
-            opening.ratio = ratio
-                
-            # Average measures of energy in front of this opening
-            Es = d + 0.5*v**2/g  #  Specific energy in exchange area
-            Et = w + 0.5*v**2/g  #  Total energy in the enquiry area
-            opening.total_energy = Et
-            opening.specific_energy = Es            
-            
-            # Store current average stage and depth with each opening object 
-            opening.depth = d
-            opening.stage = w
-            opening.elevation = z
-            
-
-        #################  End of the FOR loop ################################################
-
-            
-        # We now need to deal with each opening individually
-
-        # Determine flow direction based on total energy difference
-        delta_Et = openings[0].total_energy - openings[1].total_energy
-
-        if delta_Et > 0:
-            #print 'Flow U/S ---> D/S'
-            inlet=openings[0]
-            outlet=openings[1]
-
-            inlet.momentum = self.opening_momentum[0]
-            outlet.momentum = self.opening_momentum[1]
-
-        else:
-            #print 'Flow D/S ---> U/S'
-            inlet=openings[1]
-            outlet=openings[0]
-
-            inlet.momentum = self.opening_momentum[1]
-            outlet.momentum = self.opening_momentum[0]
-            
-            delta_Et = -delta_Et
-
-        msg = 'Total energy difference is negative'
-        assert delta_Et > 0.0, msg
-
-        delta_h = inlet.stage - outlet.stage
-        delta_z = inlet.elevation - outlet.elevation
-        culvert_slope = (delta_z/self.distance)
-
-        if culvert_slope < 0.0:
-            # Adverse gradient - flow is running uphill
-            # Flow will be purely controlled by uphill outlet face
-            print 'WARNING: Flow is running uphill. Watch Out!', inlet.elevation, outlet.elevation
-
-
-        s = 'Delta total energy = %.3f' %(delta_Et)
-        log_to_file(fid, s)
-        
-        sum_loss = self.sum_loss
-        manning=self.manning
-        length = self.distance
-        height = self.height
-        width = self.width
-                    
-        Q, barrel_velocity, culvert_outlet_depth = self.culvert_routine(inlet, outlet, delta_Et, width, height, length, sum_loss, manning, g, fid)
-
-        #####################################################
-        barrel_momentum = barrel_velocity*culvert_outlet_depth
-                
-        s = 'Barrel velocity = %f' %barrel_velocity
-        log_to_file(fid, s)
-
-        # Compute momentum vector at outlet
-        outlet_mom_x, outlet_mom_y = self.vector * barrel_momentum
-            
-        s = 'Directional momentum = (%f, %f)' %(outlet_mom_x, outlet_mom_y)
-        log_to_file(fid, s)
-
-        delta_t = time - self.last_time
-        if delta_t > 0.0:
-            xmomentum_rate = outlet_mom_x - outlet.momentum[0].value
-            xmomentum_rate /= delta_t
-                
-            ymomentum_rate = outlet_mom_y - outlet.momentum[1].value
-            ymomentum_rate /= delta_t
-                    
-            s = 'X Y MOM_RATE = (%f, %f) ' %(xmomentum_rate, ymomentum_rate)
-            log_to_file(fid, s)                    
-        else:
-            xmomentum_rate = ymomentum_rate = 0.0
-
-
-        # Set momentum rates for outlet jet
-        outlet.momentum[0].rate = xmomentum_rate
-        outlet.momentum[1].rate = ymomentum_rate
-
-        # Remember this value for next step (IMPORTANT)
-        outlet.momentum[0].value = outlet_mom_x
-        outlet.momentum[1].value = outlet_mom_y                    
-
-        if int(domain.time*100) % 100 == 0:
-            s = 'T=%.5f, Culvert Discharge = %.3f f'\
-                %(time, Q)
-            s += ' Depth= %0.3f  Momentum = (%0.3f, %0.3f)'\
-                 %(culvert_outlet_depth, outlet_mom_x,outlet_mom_y)
-            s += ' Momentum rate: (%.4f, %.4f)'\
-                 %(xmomentum_rate, ymomentum_rate)                    
-            s+='Outlet Vel= %.3f'\
-                %(barrel_velocity)
-            log_to_file(fid, s)
-            
-                
-
-
-       
-        # Execute flow term for each opening
-        # This is where Inflow objects are evaluated and update the domain
-        for opening in self.openings:
-            opening(domain)
-
-        # Execute momentum terms
-        # This is where Inflow objects are evaluated and update the domain
-        outlet.momentum[0](domain)
-        outlet.momentum[1](domain)        
-            
-        # Store value of time    
-        self.last_time = time
-
-
-
-
 #------------------------------------------------------------------------------
 # Setup culvert inlets and outlets in current topography
@@ -644 +149 @@
 # Also Set the Shape & Gap Factors fo rthe Enquiry PolyGons
 # ALSO Allow the Invert Level to be provided by the USER
-culvert = Culvert_flow(domain,
+
+
+culvert1 = Culvert_flow(domain,
                        label='Culvert No. 1',
-                       description=' This culvert is a test unit 1.2m Wide by 0.75m High',   
+                       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],
@@ -652 +159 @@
                        culvert_routine=boyd_generalised_culvert_model,        
                        verbose=True)
+
+culvert2 = Culvert_flow(domain,
+                       label='Culvert No. 2',
+                       description='This culvert is a circular test with d=1.2m',   
+                       end_point0=[9.0, 1.5], 
+                       end_point1=[30.0, 3.5],
+                       diameter=1.20,
+                       invert_level0=7,
+                       culvert_routine=boyd_generalised_culvert_model,        
+                       verbose=True)
                        
-domain.forcing_terms.append(culvert)
+domain.forcing_terms.append(culvert1)
+domain.forcing_terms.append(culvert2)
 
 
@@ -686 +204 @@
 #------------------------------------------------------------------------------
 
-temp_keep_delta_t=0.0
 
 for t in domain.evolve(yieldstep = 0.1, finaltime = 20):