Index: /anuga_work/development/culvert_flow/culvert_boyd_channel.py
===================================================================
--- /anuga_work/development/culvert_flow/culvert_boyd_channel.py	(revision 5428)
+++ /anuga_work/development/culvert_flow/culvert_boyd_channel.py	(revision 5429)
@@ -53,4 +53,7 @@
 from anuga.culvert_flows.culvert_polygons import create_culvert_polygons
 from anuga.utilities.polygon import plot_polygons
+from anuga.utilities.system_tools import log_to_file
+
+
 
 from math import pi,sqrt
@@ -60,9 +63,4 @@
 #------------------------------------------------------------------------------
 
-
-def log(fid, s):
-    print s
-    fid.write(s + '\n')
-    
 
 # Open file for storing some specific results...
@@ -89,7 +87,126 @@
 
 s='Size'+str(len(domain))
-log(fid, s)
-
-velocity_protection = 1.0e-4
+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 
+
 
 
@@ -137,21 +254,4 @@
                z[i] += embank_hgt-embank_dnslope*(x[i] -12.0)       # Sloping D/S Face
       
-        
-        # Constriction
-        #if 27 < x[i] < 29 and y[i] > 3:
-        #    z[i] += 2        
-        
-        # Pole
-        #if (x[i] - 34)**2 + (y[i] - 2)**2 < 0.4**2:
-        #    z[i] += 2
-
-        # HOLE For Pit at Opening[0]
-        #if (x[i]-4)**2 + (y[i]-1.5)**2<0.75**2:
-        #  z[i]-=1
-
-        # HOLE For Pit at Opening[1]
-        #if (x[i]-20)**2 + (y[i]-3.5)**2<0.5**2:
-        #  z[i]-=1
-
     return z
 
@@ -237,4 +337,5 @@
                  blockage_topdwn=None,
                  blockage_bottup=None,
+                 culvert_routine=None,
                  verbose=False):
         
@@ -252,4 +353,5 @@
         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        
         
 
@@ -303,4 +405,7 @@
         self.blockage_bottup = blockage_bottup
 
+        # Store culvert routine
+        self.culvert_routine = culvert_routine
+
         
         # Create objects to update momentum (a bit crude at this stage)
@@ -324,8 +429,8 @@
         # Print something
         s = 'Culvert Effective Length = %.2f m' %(self.distance)
-        log(fid, s)
+        log_to_file(fid, s)
    
         s = 'Culvert Direction is %s\n' %str(self.vector)
-        log(fid, s)        
+        log_to_file(fid, s)        
         
     def __call__(self, domain):
@@ -340,5 +445,5 @@
         delta_t = time-self.last_time
         s = '\nTime = %.2f, delta_t = %f' %(time, delta_t)
-        log(fid, s)
+        log_to_file(fid, s)
         
         msg = 'Time did not advance'
@@ -358,9 +463,4 @@
             enquiry_indices = inside_polygon(coordinates, self.enquiry_polygons[i]) 
             
-            if self.verbose:
-                pass
-                #print 'Opening %d: Got %d points in enquiry polygon:\n%s'\
-                #      %(i, len(idx), self.enquiry_polygons[i])
-                
             
             # Get model values for points in enquiry polygon for this opening
@@ -434,5 +534,6 @@
 
             inlet.momentum = self.opening_momentum[0]
-            outlet.momentum = self.opening_momentum[1]            
+            outlet.momentum = self.opening_momentum[1]
+
         else:
             #print 'Flow D/S ---> U/S'
@@ -459,164 +560,62 @@
 
         s = 'Delta total energy = %.3f' %(delta_Et)
-        log(fid, s)
-        
-
-        # ====================== NOW ENTER INTO THE CULVERT EQUATIONS AS DEFINED BY BOYD GENERALISED METHOD
-        # == 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 
-        
-        
-        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 self.width==self.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
+        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
                 
-                sum_loss=self.loss_exit+self.loss_entry+self.loss_bend+self.loss_special 
-                manning=self.manning
-
-                # Calculate flows for inlet control
-                E = inlet.specific_energy
-                #E = min(inlet.specific_energy, delta_Et)
-
-                s = 'Driving energy  = %f m' %E
-                log(fid, s)
-
-                msg = 'Driving energy is negative'
-                assert E >= 0.0, msg
-                          
-                Q_inlet_unsubmerged = 0.540*g**0.5*self.width*E**1.50 # Flow based on Inlet Ctrl Inlet Unsubmerged
-                Q_inlet_submerged = 0.702*g**0.5*self.width*self.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(fid, s)
-
-                case = ''
-                if Q_inlet_unsubmerged < Q_inlet_submerged:
-                    Q = Q_inlet_unsubmerged
-                    flow_area = self.width*inlet.depth
-                    outlet_culv_depth = inlet.depth
-                    case = 'Inlet unsubmerged'
-                else:    
-                    Q = Q_inlet_submerged
-                    flow_area = self.width*self.height
-                    outlet_culv_depth = self.height
-                    case = 'Inlet submerged'                    
-
-                if delta_Et < Es:
-                    # Calculate flows for outlet control
-                    # Determine the depth at the outlet relative to the depth of flow in the Culvert
-
-                    sum_loss = self.sum_loss
+        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
                 
-                    if outlet.depth > self.height:  # The Outlet is Submerged
-                        outlet_culv_depth=self.height
-                        flow_area=self.width*self.height # Cross sectional area of flow in the culvert
-                        perimeter=2.0*(self.width+self.height)
-                        case = 'Outlet submerged'
-                    elif outlet.depth==0.0: 
-                        outlet_culv_depth=inlet.depth   # For purpose of calculation assume the outlet depth = the inlet depth
-                        flow_area=self.width*inlet.depth
-                        perimeter=(self.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_culv_depth=outlet.depth
-                        flow_area=self.width*outlet.depth
-                        perimeter=(self.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(fid, s)
+            ymomentum_rate = outlet_mom_y - outlet.momentum[1].value
+            ymomentum_rate /= delta_t
                     
-
-                    # Outlet control velocity using tail water
-                    culvert_velocity = sqrt(delta_Et/((sum_loss/2*g)+(manning**2*self.distance)/hyd_rad**1.33333)) 
-                    Q_outlet_tailwater = flow_area * culvert_velocity
-
-                    s = 'Q_outlet_tailwater = %.6f' %Q_outlet_tailwater
-                    log(fid, s)
-
-                    Q = min(Q, Q_outlet_tailwater)
-
-
-                    
-                log(fid, 'Case: "%s"' %case)
+            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)
+            
                 
-                flow_rate_control=Q
-
-                s = 'Flow Rate Control = %f' %flow_rate_control
-                log(fid, s)
-
-                inlet.rate = -flow_rate_control
-                outlet.rate = flow_rate_control                
-                
-                culv_froude=sqrt(flow_rate_control**2*self.width/(g*flow_area**3))
-                s = 'Froude in Culvert = %f' %culv_froude
-                log(fid, s)
-
-                
-
-                # Determine momentum at the outlet 
-                barrel_velocity = Q/(flow_area + velocity_protection/flow_area)
-                barrel_momentum = barrel_velocity*outlet_culv_depth
-                
-                outlet_E_Loss= 0.5*0.5*barrel_velocity**2/g   # Ke v^2/2g
-                s = 'Barrel velocity = %f' %barrel_velocity
-                log(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(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(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  Culv. Vel. %0.3f'\
-                        %(time, flow_rate_control, barrel_velocity)
-                    s += ' Depth= %0.3f  Momentum = (%0.3f, %0.3f)'\
-                         %(outlet_culv_depth, outlet_mom_x,outlet_mom_y)
-                    s += ' Momentum rate: (%.4f, %.4f)'\
-                         %(xmomentum_rate, ymomentum_rate)                    
-                    s+='Outlet Vel= %.3f'\
-                         %(barrel_velocity)
-                    log(fid, s)
-         
+
+
        
         # Execute flow term for each opening
@@ -650,5 +649,6 @@
                        end_point0=[9.0, 2.5], 
                        end_point1=[13.0, 2.5],
-                       width=1.20,height=0.75,                       
+                       width=1.20,height=0.75,
+                       culvert_routine=boyd_generalised_culvert_model,        
                        verbose=True)
                        
