Index: anuga_work/development/Rudy_vandrie_2007/culvert_development.py
===================================================================
--- anuga_work/development/Rudy_vandrie_2007/culvert_development.py	(revision 5305)
+++ anuga_work/development/Rudy_vandrie_2007/culvert_development.py	(revision 5305)
@@ -0,0 +1,591 @@
+""" 
+
+Ole Check Culvert Routine from Line 258
+
+Although it is Setup as a Culvert with the Opening presented vertically, 
+for now the calculation of flow rate is assuming a horizontal hole in the
+ground (Fix this Later)
+
+MOST importantly 2 things...
+1. How to use the Create Polygon Routine to enquire Depth ( or later energy)
+infront of the Culvert
+
+Done (Ole)
+
+2. How to apply the Culvert velocity and thereby Momentum to the Outlet
+Ject presented at the Outlet
+
+
+
+Testing CULVERT (Changing from Horizontal Abstraction to Vertical Abstraction
+
+This Version CALCULATES the Culvert Velocity and Uses it to establish 
+The Culvert Outlet Momentum
+
+The Aim is to define a Flow Transfer function that Simulates a Culvert
+by using the Total Available Energy to Drive the Culvert
+as Derived by determining the Difference in Total Energy
+between 2 Points, Just Up stream and Just Down Stream of the Culvert
+away from the influence of the Flow Abstraction etc..
+
+This example includes a Model Topography that shows a
+TYPICAL Headwall Configuration
+
+The aim is to change the Culvert Routine to Model more precisely the
+abstraction
+from a vertical face.
+
+The inflow must include the impact of Approach velocity.
+Similarly the Outflow has MOMENTUM Not just Up welling as in the
+Horizontal Style
+abstraction
+
+"""
+
+#------------------------------------------------------------------------------
+# Import necessary modules
+#------------------------------------------------------------------------------
+from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular_cross
+from anuga.shallow_water import Domain
+from anuga.shallow_water.shallow_water_domain import Reflective_boundary
+from anuga.shallow_water.shallow_water_domain import Dirichlet_boundary
+from anuga.shallow_water.shallow_water_domain import Inflow, General_forcing
+from anuga.culvert_flows.culvert_polygons import create_culvert_polygons
+from anuga.utilities.polygon import plot_polygons
+
+from math import pi,sqrt
+
+#------------------------------------------------------------------------------
+# Setup computational domain
+#------------------------------------------------------------------------------
+
+
+# Open file for storing some specific results...
+fid = open('Culvert_Headwall_VarM.txt', 'w')
+
+length = 40.
+width = 5.
+
+#dx = dy = 1           # Resolution: Length of subdivisions on both axes
+#dx = dy = .5           # Resolution: Length of subdivisions on both axes
+dx = dy = .25           # Resolution: Length of subdivisions on both axes
+#dx = dy = .1           # Resolution: Length of subdivisions on both axes
+
+# OLE.... How do I refine the resolution... in the area where I have the Culvert Opening ???......
+#  Can I refine in a X & Y Range ???
+points, vertices, boundary = rectangular_cross(int(length/dx), int(width/dy),
+                                               len1=length, len2=width)
+domain = Domain(points, vertices, boundary)   
+domain.set_name('culvert_HW_Var_Mom')                 # Output name
+domain.set_default_order(2)
+domain.H0 = 0.01
+domain.tight_slope_limiters = True
+domain.set_minimum_storable_height(0.02)
+print 'Size', len(domain)
+
+
+#------------------------------------------------------------------------------
+# Setup initial conditions
+#------------------------------------------------------------------------------
+
+# Define the topography (land scape)
+def topography(x, y):
+    """Set up a weir
+    
+    A culvert will connect either side of an Embankment with a Headwall type structure
+    The aim is for the Model to REALISTICALY model flow through the Culvert
+    """
+    # General Slope of Topography
+    z=-x/100
+    
+    # Add bits and Pieces to topography
+    N = len(x)
+    for i in range(N):
+
+        # Sloping Embankment Across Channel
+        if 5.0 < x[i] < 10.1:
+            if  1.0+(x[i]-5.0)/5.0 <  y[i]  < 4.0 - (x[i]-5.0)/5.0: # Cut Out Segment for Culvert FACE
+               z[i]=z[i]
+            else:
+               z[i] +=  0.5*(x[i] -5.0)    # Sloping Segment  U/S Face
+        if 10.0 < x[i] < 12.1:
+           z[i] +=  2.5              # Flat Crest of Embankment
+        if 12.0 < x[i] < 14.5:
+            if  2.0-(x[i]-12.0)/2.5 <  y[i]  < 3.0 + (x[i]-12.0)/2.5: # Cut Out Segment for Culvert FACE
+               z[i]=z[i]
+            else:
+               z[i] += 2.5-1.0*(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
+
+
+domain.set_quantity('elevation', topography)  # Use function for elevation
+domain.set_quantity('friction', 0.01)         # Constant friction 
+domain.set_quantity('stage',
+                    expression='elevation')   # Dry initial condition
+
+
+
+
+#------------------------------------------------------------------------------
+# 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 Verical Opening Both BOX and Circular
+ # There will be several Culvert Routines such as:
+ # 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,
+                 end_point0=None, 
+                 end_point1=None,
+                 width=None,
+                 height=None,
+                 verbose=False):
+        
+        from Numeric import sqrt, sum
+
+        if height is None: height = width
+
+        # Create the fundamental culvert polygons
+        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.enquiry_polygons = [P['enquiry_polygon0'], P['enquiry_polygon1']]
+        self.vector = P['vector']
+        self.distance = P['length']
+        self.verbose = verbose
+        self.width = width
+        self.height = height
+        self.last_time = 0.0        
+
+        # Create objects to update momentum (a bit crude at this stage)
+        self.xmom_forcing0 = General_forcing(domain,
+                                             'xmomentum',
+                                             polygon=P['exchange_polygon0'])
+
+        self.xmom_forcing1 = General_forcing(domain,
+                                             'xmomentum',
+                                             polygon=P['exchange_polygon1'])
+
+        self.ymom_forcing0 = General_forcing(domain,
+                                             'ymomentum',
+                                             polygon=P['exchange_polygon0'])
+
+        self.ymom_forcing1 = General_forcing(domain,
+                                             'ymomentum',
+                                             polygon=P['exchange_polygon1'])                
+
+        # Print something
+        
+        print 'Culvert Effective Length =', self.distance
+        print 'Culvert Slope is Delta Z / Dist'
+        print 'Culvert Direction is ', self.vector
+        print 'Point 1m Up Stream is X,Y ='
+        print 'Point 1m Down Stream is X,Y ='
+       
+
+        
+    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
+        
+    
+        # Get average water depths at each opening
+
+        
+        time = domain.get_time()
+        openings = self.openings
+        for i, opening in enumerate(openings):
+            stage = domain.quantities['stage'].get_values(location='centroids',
+                                                          indices=opening.indices)
+            elevation = domain.quantities['elevation'].get_values(location='centroids',
+                                                                  indices=opening.indices)
+
+            # Indices corresponding to energy enquiry field for this opening
+            coordinates = domain.get_centroid_coordinates() # Get all centroid points (x,y)
+            idx = 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 average model values for points in
+            # enquiry polygon for this opening
+            dq = domain.quantities
+            stage = dq['stage'].get_values(location='centroids', indices=idx)
+            xmomentum = dq['xmomentum'].get_values(location='centroids', indices=idx)
+            ymomentum = dq['ymomentum'].get_values(location='centroids', indices=idx)                        
+            elevation = dq['elevation'].get_values(location='centroids', indices=idx)
+            depth = stage - elevation            
+
+            # Compute mean velocity in the area in front of the culvert (taking zero depths into account)
+            ux = xmomentum/(depth+epsilon)
+            uy = ymomentum/(depth+epsilon)
+            v = mean(sqrt(ux**2+uy**2))
+            d = mean(depth)
+            w = mean(stage)
+            z = mean(elevation)
+            
+            # Ratio of depth to culvert height
+            ratio = d/(2*self.height)
+            if ratio > 1.0:    # Assume culvert is running full & 
+                ratio = 1.0    # under pressure. Note this is usually ~ 1.35
+            opening.ratio = ratio
+                
+
+            # Average measure of total energy (D + V^2/2g) in enquiry field in front of this opening
+            E = d + 0.5*v**2/g  # RUDY - please check this
+            opening.energy = E
+            
+            
+
+            # Store current average stage and depth at enquiry field with each opening object 
+            opening.depth = d
+            opening.stage = w
+            opening.elevation = z
+            
+
+        # Now we are done calculating energy etc for each opening.
+        # At this point we can work on the transfer functions etc.
+        
+        # Handy values (all calculated at enquiry polygon -
+        # if you need them at exchange polygons we can easily do that.)
+        delta_h = openings[1].stage - openings[0].stage
+        #delta_h = openings[1].depth - openings[0].depth
+        delta_z = openings[1].elevation - openings[0].elevation
+
+        # Ideas.....
+        if openings[0].depth > 0 and openings[1].depth > 0:
+            flow_rate = delta_h
+        else:
+            flow_rate = 0.0
+
+        if openings[0].depth > self.height:
+            # This could be usefull mayby
+            #print 'Inlet has been overflowed'
+            pass
+
+        if openings[1].depth > self.height:
+            # This could be usefull mayby
+            #print 'Outlet has been overflowed'            
+            pass
+
+		
+        if delta_h > 0:
+            # Water Level U/S is higher than DS
+            flow_rate_orifice = 0.6*openings[0].area*(2*g*delta_h)**0.5      # Orifice Eqn Q= cA(2gh)^0.5
+            flow_rate_weir = 1.69*(openings[0].area)*openings[0].depth**1.5  # WEIR Eqn   Q= CLH^1.5
+
+            # Does Weir or Orifice Control Flow Rate ?
+            if flow_rate_weir > flow_rate_orifice:   
+                flow_rate_control = flow_rate_orifice
+            else:
+                flow_rate_control = flow_rate_weir
+            
+        elif delta_h < 0:
+            # Water Level D/S is higher than US
+            # That is reverse flow in culvert
+            flow_rate_orifice = 0.6*openings[0].area*(2*g*-delta_h)**0.5     # Orifice Eqn Q= cA(2gh)^0.5
+            flow_rate_weir = 1.69*openings[0].area*openings[0].depth**1.5    # WEIR Eqn   Q= CLH^1.5
+
+            # Does Weir or Orifice Control Flow Rate ?
+            if flow_rate_weir > flow_rate_orifice:  
+                flow_rate_control = flow_rate_orifice
+            else:
+                flow_rate_control = flow_rate_weir
+
+        if openings[0].depth < epsilon:
+            # Do nothing because No water over Opening	That is Set Flow to Zero!!
+            self.openings[0].flow = 0
+            self.openings[1].flow = 0
+  
+        elif openings[0].depth > 0:
+            # Only calculate fllow if there is some depth over the inlet
+            #if delta_h > 0:
+            if 1:
+                # FIXME (OLE): We never get in here. Why?
+                
+                #print 'We got water at inlet and dh > 0'
+            
+                # Flow will go from opening 0 to opening 1
+                #- that is abstract from [0] and add to [1]
+                self.openings[0].flow = -flow_rate_control
+                self.openings[1].flow = flow_rate_control
+
+                # Add jet in the form of absolute momentum to opening 1
+                #speed = 20.0 # m/s
+                # This Should be Based on the VELOCITY in the Culvert 
+                # Based on the Flow Depth in the Culvert for part full
+                # or Flow Jet of the Pressurised Culvert if FUll
+                # If Part Full Flow Calculate Part Full Depth
+                # Based on Depth Calculate Area... the Vel = flow_rate_control / Area
+                # Need to Break Velocity into X & Y Components
+                #self.openings[1].set_quantity_values(delta_h*speed, 'xmomentum')
+                # Previous Calculated Depth/ Culvert Heigh Ratio Use it to Determine Velocity !!!! FIX LAter
+                
+                #if (ratio*self.area)== 0: # Don't Really need this already established water depth here
+                #   outlet_vel=0.0
+                #else:
+
+
+                # FIXME (Ole): Shouldn't this be openings[1]??
+                outlet_vel =(flow_rate_control/(openings[0].ratio*openings[0].area))
+                outlet_dep = 2.0*openings[0].depth*openings[0].ratio #?????
+                outlet_mom = outlet_vel*outlet_dep
+
+                
+                # Eventually will Need Momentum in X & Y Components based on the orientation of
+                # the culvert from X0,Y0, X1, Y1 from Create Polygon Routine
+                # YES - use self.vector which is a unit vector in the direction of the culvert.
+                    
+                outlet_mom_x, outlet_mom_y = self.vector * outlet_mom
+                #print 'Directional momentum', outlet_mom_x, outlet_mom_y
+                    
+                # Update the momentum forcing terms with directional momentum at the outlet
+                delta_t = time - self.last_time
+                if delta_t > 0.0:
+                    xmomentum_rate = outlet_mom_x - self.xmom_forcing1.value
+                    xmomentum_rate /= delta_t
+                    
+                    ymomentum_rate = outlet_mom_y - self.ymom_forcing1.value
+                    ymomentum_rate /= delta_t                        
+                else:
+                    xmomentum_rate = ymomentum_rate = 0.0
+
+                # Set momentum rates for outlet jet
+                self.xmom_forcing1.rate = xmomentum_rate
+                self.ymom_forcing1.rate = ymomentum_rate
+
+                # Remember this value for next step (IMPORTANT)    
+                self.xmom_forcing1.value = outlet_mom_x
+                self.ymom_forcing1.value = outlet_mom_y                    
+
+
+                if int(domain.time*100) % 100 == 0:
+                    s = 'T=%.2f, Culvert Discharge = %.3f  Culv. Vel. %0.3f'\
+                        %(time, flow_rate_control, outlet_vel)
+                    s += ' Depth= %0.3f  Momentum = %0.3f' %(outlet_dep, outlet_mom)
+                    fid.write(s)
+                    print s
+                    
+
+            else:
+                # Opening 1 has the greatest depth. Therefore Reverse the Flow !!!
+                # Flow will go from opening 1 to opening 0, That is Abstract from [1] and add to [0]
+                self.openings[0].flow = flow_rate_control        # Else it will be Orifice Flow (Going US)
+                self.openings[1].flow = -flow_rate_control
+
+        # Second Else.... if water at outlet before at inlet !!!
+
+
+       
+        # 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
+        self.xmom_forcing0(domain)
+        self.ymom_forcing0(domain)                
+        self.xmom_forcing1(domain)
+        self.ymom_forcing1(domain)        
+
+            
+
+        # Print out flows every 1 seconds
+        if int(time*100) % 100 == 0:
+            s = 'Time %.2f\n' %time
+            s += '    Opening[0]: d=%.3f, A=%.3f, E=%.3f, r=%.3f\n'\
+                 %(openings[0].depth,
+                   openings[0].area,
+                   openings[0].energy,
+                   openings[0].ratio)                   
+            s += '    Opening[1]: d=%.3f, A=%.3f, E=%.3f, r=%.3f\n'\
+                 %(openings[1].depth,
+                   openings[1].area,
+                   openings[1].energy,
+                   openings[1].ratio)                  
+            s += '    Distance=%.2f, W=%.3f, O=%.3f, C=%.3f\n'\
+                 %(self.distance,
+                   flow_rate_weir,
+                   flow_rate_orifice,
+                   flow_rate_control)
+            
+            print s
+            fid.write(s)
+
+        # Store value of time    
+        self.last_time = time
+
+#------------------------------------------------------------------------------
+# Setup culvert inlets and outlets in current topography
+#------------------------------------------------------------------------------
+
+# Define culvert inlet and outlets
+culvert = Culvert_flow(domain,
+                       end_point0=[9.0, 2.5], 
+                       end_point1=[13.0, 2.5],
+                       width=1.00,                       
+                       verbose=True)
+                       
+domain.forcing_terms.append(culvert)
+
+
+#------------------------------------------------------------------------------
+# Setup boundary conditions
+#------------------------------------------------------------------------------
+#Bi = Dirichlet_boundary([0.5, 0.0, 0.0])          # Inflow based on Flow Depth (0.5m) and Approaching Momentum !!!
+Bi = Dirichlet_boundary([0.0, 0.0, 0.0])          # Inflow based on Flow Depth and Approaching Momentum !!!
+Br = Reflective_boundary(domain)              # Solid reflective wall
+Bo = Dirichlet_boundary([-5, 0, 0])           # Outflow
+
+domain.set_boundary({'left': Br, 'right': Bo, 'top': Br, 'bottom': Br})
+
+#------------------------------------------------------------------------------
+# Setup Application of  specialised forcing terms
+#------------------------------------------------------------------------------
+
+# This is the new element implemented by Ole to allow direct input of Inflow in m^3/s
+fixed_flow = Inflow(domain,
+                    rate=6.00,
+                    center=(2.1, 2.1),
+                    radius=1.261566)   #   Fixed Flow Value Over Area of 5m2 at 1m/s = 5m^3/s
+
+#	     flow=file_function('Q/QPMF_Rot_Sub13.tms'))        # Read Time Series in  from File
+#             flow=lambda t: min(0.01*t, 0.01942))                             # Time Varying Function   Tap turning up  	
+
+domain.forcing_terms.append(fixed_flow)
+
+
+#------------------------------------------------------------------------------
+# Evolve system through time
+#------------------------------------------------------------------------------
+
+
+
+for t in domain.evolve(yieldstep = 0.1, finaltime = 10):
+    pass
+    #if int(domain.time*100) % 100 == 0:
+    #    domain.write_time()
+    
+    #if domain.get_time() >= 4 and tap.flow != 0.0:
+    #    print 'Turning tap off'
+    #    tap.flow = 0.0
+	
+    #if domain.get_time() >= 3 and sink.flow < 0.0:
+    #    print 'Turning drain on'
+    #    sink.flow = -0.8	
+    # Close
+
+fid.close()
+
+
+#------------------------------------------------------------------------------
+# Query output
+#------------------------------------------------------------------------------
+
+from anuga.shallow_water.data_manager import get_flow_through_cross_section
+
+swwfilename = domain.get_name()+'.sww'  # Output name from script
+print swwfilename
+
+polyline = [[17., 0.], [17., 5.]]
+
+time, Q = get_flow_through_cross_section(swwfilename, polyline, verbose=True)
+
+from pylab import ion, plot
+ion()
+plot(time, Q)
+raw_input('done')
