Changeset 4697 for anuga_work/development/Rudy_vandrie_2007/culvert_flow.py

Timestamp:

Sep 3, 2007, 6:29:38 PM (17 years ago)

Author:

ole

Message:

Culvert flow prototype for Rudy

File:

• anuga_work/development/Rudy_vandrie_2007/culvert_flow.py (modified) (5 diffs)

anuga_work/development/Rudy_vandrie_2007/culvert_flow.py

@@ -12 +12 @@
 # Import necessary modules
 #------------------------------------------------------------------------------
-from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular
+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
 
+
 #------------------------------------------------------------------------------
 # Setup computational domain
 #------------------------------------------------------------------------------
-
-#N = 120
-N = 40
-points, vertices, boundary = rectangular(N, N)
+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 = .1           # Resolution: Length of subdivisions on both axes
+
+points, vertices, boundary = rectangular_cross(int(length/dx), int(width/dy),
+                                               len1=length, len2=width)
 domain = Domain(points, vertices, boundary)   
-domain.set_name('kitchen_sink')                 # Output name
+domain.set_name('culvert_example')                 # Output name
 print 'Size', len(domain)
 
-
 #------------------------------------------------------------------------------
 # Setup initial conditions
 #------------------------------------------------------------------------------
 
-domain.set_quantity('elevation', 0.0)  # Zero bed elevation
-domain.set_quantity('stage', 0.015)    # Constant stage
-domain.set_quantity('friction', 0.005) # Constant friction 
+def topography(x, y):
+    """Set up a weir
+    
+    A culvert will connect either side
+    """
+
+    z = -x/10                                
+
+    N = len(x)
+    for i in range(N):
+
+        # Weir from 10 to 12 m
+        if 10 < x[i] < 12:
+            z[i] += 0.4 - 0.05*y[i]        
+        
+        # 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
+
+    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
+
+
 
 
@@ -67 +101 @@
                                                      # specified area 
     """
-    
+
     # FIXME (OLE): Add a polygon as an alternative.
     # FIXME (OLE): Generalise to all quantities
-
-    def __init__(self, 
+    # FIXME (OLE): 
+
+    def __init__(self,
+                 domain,
                  center=None, radius=None,
                  flow=0.0,
                  quantity_name = 'stage'):
-                 
+
+        from anuga.utilities.numerical_tools import ensure_numeric
+    
         if center is not None and radius is not None:
             assert len(center) == 2
@@ -81 +119 @@
             msg = 'Both center and radius must be specified'
             raise Exception, msg
-    
-        self.center = center
+
+        self.domain = domain
+        self.center = ensure_numeric(center)
         self.radius = radius
         self.flow = flow
+        
         self.quantity = domain.quantities[quantity_name].explicit_update
         
+        # Determine indices in flow area
+        N = len(domain)    
+        self.indices = []
+        coordinates = domain.get_centroid_coordinates() 
+        for k in range(N):
+            x, y = coordinates[k,:] # Centroid
+            if ((x-center[0])**2+(y-center[1])**2) < radius**2:
+                self.indices.append(k)    
+
     
     def __call__(self, domain):
     
-        # Determine indices in flow area
-        if not hasattr(self, 'indices'):
-            center = self.center
-            radius = self.radius
-            
-            N = len(domain)    
-            self.indices = []
-            coordinates = domain.get_centroid_coordinates()     
-            for k in range(N):
-                x, y = coordinates[k,:] # Centroid
-                if ((x-center[0])**2+(y-center[1])**2) < radius**2:
-                    self.indices.append(k)    
-
         # Update inflow
         if callable(self.flow):
@@ -117 +153 @@
     """ 
 
-    def __init__(self, 
-                 inflow_center=None, inflow_radius=None,                    
-                 outflow_center=None, outflow_radius=None,
-                 transfer_flow=None)
-                 
-        self.inflow=Inflow(center=inflow_center,
-                           radius=inflow_radius,
-                           flow=transfer_flow)
+    def __init__(self,
+                 domain,
+                 center0=None, radius0=None,                
+                 center1=None, radius1=None,
+                 transfer_flow=None):
+        
+        from Numeric import sqrt, sum
+
+        self.openings = []
+        self.openings.append(Inflow(domain,
+                                    center=center0,
+                                    radius=radius0,
+                                    flow=None))
                            
-        self.outflow=Inflow(center=outflow_center,
-                           radius=outflow_radius,
-                           flow=-transfer_flow) #!!!                       
+        self.openings.append(Inflow(domain,
+                                    center=center1,
+                                    radius=radius1,
+                                    flow=None))                    
                            
-                                
+        # Compute distance between opening centers. I assume there are only two.
+        self.distance = sqrt(sum( (self.openings[0].center - self.openings[1].center)**2 ))
+        print 'Distance between openings is', self.distance
+        
                  
     def __call__(self, domain):
-        self.inflow(domain)              
-        self.outflow(domain)                    
+        from anuga.utilities.numerical_tools import mean
+
+    
+        # Get average water depths at each opening
+        for opening in self.openings:
+            stage = domain.quantities['stage'].get_values(location='centroids',
+                                                          indices=opening.indices)
+            elevation = domain.quantities['elevation'].get_values(location='centroids',
+                                                              indices=opening.indices)
+
+            # Store current avg depth with opening object
+            opening.depth = mean(stage-elevation)
+            # print 'Depth', opening.depth
+
+
+        # Here's a simple transfer function hardwired into this Culvert class:
+        # Let flow rate depend on difference in depth
+        # Flow is instantaneous for now but one might use
+        # the precomputed distance between the two openings
+        # to somehow introduce a delay.
+
+        delta_h = self.openings[0].depth - self.openings[1].depth
+        flow_rate = 0.5*9.81*delta_h**2                          # Just an example of flow rate
+        flow_rate /= (self.openings[0].radius + self.openings[1].radius)/2 #  
+
+        if delta_h > 0:
+            # Opening 0 has the greatest depth.
+            # Flow will go from opening 0 to opening 1
+
+            self.openings[0].flow = -flow_rate
+            self.openings[1].flow = flow_rate
+        else:
+            # Opening 1 has the greatest depth.
+            # Flow will go from opening 1 to opening 0
+
+            self.openings[0].flow = flow_rate
+            self.openings[1].flow = -flow_rate
+            
+
+        # Execute flow term for each opening 
+        for opening in self.openings:
+            opening(domain)
+
                                                          
                     
-sink = Inflow(center=[0.7, 0.4], radius=0.0707, flow = 0.0)             
-tap = Inflow(center=(0.5, 0.5), radius=0.0316,
-             flow=lambda t: min(4*t, 5)) # Tap turning up       
-        
-
-domain.forcing_terms.append(tap)
-domain.forcing_terms.append(sink)
+#sink = Inflow(domain, center=[0.7, 0.4], radius=0.0707, flow = 0.0)            
+#tap = Inflow(domain, center=(0.5, 0.5), radius=0.0316,
+#            flow=lambda t: min(4*t, 5)) # Tap turning up       
+#domain.forcing_terms.append(tap)
+#domain.forcing_terms.append(sink)
+
+culvert = Culvert_flow(domain,
+                       center0=[5.0, 2.5], radius0=0.3,
+                       center1=[20., 2.5], radius1=0.3,
+                       transfer_flow=None) # Hardwired for now
+                       
+domain.forcing_terms.append(culvert)
 
 #------------------------------------------------------------------------------
 # Setup boundary conditions
 #------------------------------------------------------------------------------
-
+Bi = Dirichlet_boundary([0.4, 0, 0])          # Inflow
 Br = Reflective_boundary(domain)              # Solid reflective wall
-domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
+Bo = Dirichlet_boundary([-5, 0, 0])           # Outflow
+
+domain.set_boundary({'left': Bi, 'right': Bo, 'top': Br, 'bottom': Br})
+
 
 #------------------------------------------------------------------------------
@@ -158 +252 @@
     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() >= 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        
-    
+    #if domain.get_time() >= 3 and sink.flow < 0.0:
+    #    print 'Turning drain on'
+    #    sink.flow = -0.8       
+