Changeset 7962

Timestamp:

Aug 20, 2010, 5:27:03 PM (15 years ago)

Author:

steve

Message:

Got a flow from one region to another.

Location:

trunk/anuga_core/source/anuga/structures

Files:

• culvert_operator.py (modified) (4 diffs)
• testing_culvert.py (modified) (7 diffs)

trunk/anuga_core/source/anuga/structures/culvert_operator.py

@@ -46 +46 @@
         self.domain = domain
 
+        self.domain.set_fractional_step_operator(self)
+
         self.end_points= [end_point0, end_point1]
         self.enquiry_gap_factor = enquiry_gap_factor
@@ -62 +64 @@
         self.compute_enquiry_indices()
         self.check_culvert_inside_domain()
-
-        # Establish initial values at each enquiry point
-        self.enquiry_quantity_values = []
-        dq = domain.quantities 
-        for i in [0, 1]:
-            idx = self.enquiry_indices[i]
-            elevation = dq['elevation'].get_values(location='centroids',
-                                                   indices=[idx])[0]
-            stage = dq['stage'].get_values(location='centroids',
-                                           indices=[idx])[0]
-            depth = stage - elevation
-            
-            quantity_values = {'stage' : stage, 'elevation' : elevation, 'depth' : depth }
-            self.enquiry_quantity_values.append(quantity_values)
+        self.compute_exchange_triangle_indices()
+
+
+        #self.print_stats()
+
+
+
+    def __call__(self):
+
+
+        # Time stuff
+        time     = self.domain.get_time()
+        timestep = self.domain.get_timestep()
+
+
+        inlet_indices  = self.exchange_triangle_indices[0]
+        outlet_indices = self.exchange_triangle_indices[1]
+
+        areas      = self.domain.areas
+        stage     = self.domain.quantities['stage'].centroid_values
+        elevation = self.domain.quantities['elevation'].centroid_values
+        
+        xmom = self.domain.quantities['xmomentum'].centroid_values
+        ymom = self.domain.quantities['ymomentum'].centroid_values
+
+        # Inlet averages
+        inlet_heights  = stage[inlet_indices]-elevation[inlet_indices]
+        inlet_areas = areas[inlet_indices]
+
+        inlet_water = num.sum(inlet_heights*inlet_areas)
+
+        average_inlet_water = inlet_water/self.exchange_areas[0]
+
+        # Outlet averages
+        outlet_heights  = stage[outlet_indices]-elevation[outlet_indices]
+        outlet_areas = areas[outlet_indices]
+
+        outlet_water = num.sum(outlet_heights*outlet_areas)
+
+        average_outlet_water = outlet_water/self.exchange_areas[1]
+
+
+        # Transfer
+        transfer_water = timestep*inlet_water
+
+        stage[inlet_indices] = elevation[inlet_indices] + average_inlet_water - transfer_water
+        xmom[inlet_indices] = 0.0
+        ymom[inlet_indices] = 0.0
+
+
+        stage[outlet_indices] = elevation[outlet_indices] + average_outlet_water + transfer_water
+        xmom[outlet_indices] = 0.0
+        ymom[outlet_indices] = 0.0
+
+
+    def print_stats(self):
+
+        print '====================================='
+        print 'Generic Culvert Operator'
+        print '====================================='
+        print "enquiry_gap_factor"
+        print self.enquiry_gap_factor
+        
+        for i in [0,1]:
+            print '-------------------------------------'
+            print 'exchange_region %i' % i
+            print '-------------------------------------'
+
+            print 'exchange triangle indices and centres'
+            print self.exchange_triangle_indices[i]
+            print self.domain.get_centroid_coordinates()[self.exchange_triangle_indices[i]]
+        
+            print 'end_point'
+            print self.end_points[i]
+
+
+            print 'exchange_polygon'
+            print self.exchange_polygons[i]
+
+            print 'enquiry_point'
+            print self.enquiry_points[i]
+
+        print '====================================='
+
+
+ 
+
+
+    def compute_exchange_triangle_indices(self):
+
+        # Get boundary (in absolute coordinates)
+        domain = self.domain
+        bounding_polygon = domain.get_boundary_polygon()
+        centroids = domain.get_centroid_coordinates(absolute=True)
+        self.exchange_triangle_indices = []
+        self.exchange_areas = []
+
+        for i in [0,1]:
+            exchange_polygon = self.exchange_polygons[i]
+
+            # Check that polygon lies within the mesh.
+            for point in exchange_polygon:
+                msg = 'Point %s in polygon for forcing term' % str(point)
+                msg += ' did not fall within the domain boundary.'
+                assert is_inside_polygon(point, bounding_polygon), msg
+
+            
+            exchange_triangle_indices = inside_polygon(centroids, exchange_polygon)
+
+            if len(exchange_triangle_indices) == 0:
+                region = 'polygon=%s' % (exchange_polygon)
+                msg = 'No triangles have been identified in '
+                msg += 'specified region: %s' % region
+                raise Exception, msg
+
+            # Compute exchange area as the sum of areas of triangles identified
+            # by polygon
+            exchange_area = 0.0
+            for j in exchange_triangle_indices:
+                exchange_area += domain.areas[j]
+
+
+            msg = 'Exchange area %f in culvert' % i
+            msg += ' has area = %f' % exchange_area
+            assert exchange_area > 0.0
+
+            self.exchange_triangle_indices.append(exchange_triangle_indices)
+            self.exchange_areas.append(exchange_area)
+
+
 
     def set_store_hydrograph_discharge(self, filename=None):
@@ -123 +241 @@
         # Build exchange polygon and enquiry points 0 and 1
         for i in [0, 1]:
+            i0 = (2*i-1)
             p0 = self.end_points[i] + w
             p1 = self.end_points[i] - w
-            p2 = p1 - h
-            p3 = p0 - h
+            p2 = p1 + i0*h
+            p3 = p0 + i0*h
             self.exchange_polygons.append(num.array([p0, p1, p2, p3]))
-            self.enquiry_points.append(self.end_points[i] + (2*i-1)*gap)
-
-        self.polygon_areas = []
+            self.enquiry_points.append(self.end_points[i] + i0*gap)
+
+
+
 
         # Check that enquiry points are outside exchange polygons
         for i in [0,1]:
             polygon = self.exchange_polygons[i]
-            # FIXME(SR) should use area of triangles associated with polygon
+            # FIXME (SR) Probably should calculate the area of all the triangles
+            # associated with this polygon, as there is likely to be some
+            # inconsistency between triangles and ploygon
             area = polygon_area(polygon)
-            self.polygon_areas.append(area)
+            
 
             msg = 'Polygon %s ' %(polygon)
@@ -150 +272 @@
 
         
-    def __call__(self, domain):
-
-        # Time stuff
-        time = domain.get_time()
-        
-        
-        update = False
-        if self.update_interval is None:
-            # Use next timestep as has been computed in domain.py        
-            delta_t = domain.timestep            
-            update = True
-        else:    
-            # Use update interval 
-            delta_t = self.update_interval            
-            if time - self.last_update > self.update_interval or time == 0.0:
-                update = True
-            
-        if self.log_filename is not None:        
-            s = '\nTime = %.2f, delta_t = %f' %(time, delta_t)
-            log_to_file(self.log_filename, s)
-                
-                                
-        if update is True:
-            self.compute_rates(delta_t)
-            
-    
-        # Execute flow term for each opening
-        # This is where Inflow objects are evaluated using the last rate 
-        # that has been calculated
-        # 
-        # This will take place at every internal timestep and update the domain
-        for opening in self.openings:
-            opening(domain)
-            
 
 

trunk/anuga_core/source/anuga/structures/testing_culvert.py

@@ -33 +33 @@
 width = 5.
 
-dx = dy = 1           # Resolution: Length of subdivisions on both axes
+dx = dy = 0.5          # Resolution: Length of subdivisions on both axes
 
 points, vertices, boundary = rectangular_cross(int(length/dx),
@@ -85 +85 @@
 
 filename=os.path.join(path, 'example_rating_curve.csv')
-culvert = Generic_box_culvert(domain,
+culvert1 = Generic_box_culvert(domain,
                               end_point0=[9.0, 2.5], 
                               end_point1=[13.0, 2.5],
@@ -91 +91 @@
                               verbose=False)
 
-#domain.forcing_terms.append(culvert)
 
+#culvert2 = Generic_box_culvert(domain,
+#                              end_point0=[19.0, 2.5],
+#                              end_point1=[25.0, 2.5],
+#                              width=1.00,
+#                              verbose=False)
+
+
+
+
+#print domain.fractional_step_operators
+
+#domain.apply_fractional_steps()
 
 ##-----------------------------------------------------------------------
@@ -99 +110 @@
 
 ## Inflow based on Flow Depth and Approaching Momentum
-#Bi = anuga.Dirichlet_boundary([0.0, 0.0, 0.0])
-#Br = anuga.Reflective_boundary(domain)              # Solid reflective wall
+Bi = anuga.Dirichlet_boundary([1.0, 0.0, 0.0])
+Br = anuga.Reflective_boundary(domain)              # Solid reflective wall
 #Bo = anuga.Dirichlet_boundary([-5, 0, 0])           # Outflow
 
@@ -111 +122 @@
             #lambda t: [-5*(num.cos(2*pi*(t-4)/20)), 0.0, 0.0])
 #domain.set_boundary({'left': Btus, 'right': Btds, 'top': Br, 'bottom': Br})
+domain.set_boundary({'left': Bi, 'right': Br, 'top': Br, 'bottom': Br})
 
 
@@ -119 +131 @@
 #min_delta_w = sys.maxint 
 #max_delta_w = -min_delta_w
-#for t in domain.evolve(yieldstep = 1, finaltime = 25):
+for t in domain.evolve(yieldstep = 1.0, finaltime = 300):
+    domain.write_time()
     #delta_w = culvert.inlet.stage - culvert.outlet.stage
     
@@ -125 +138 @@
     #if delta_w < min_delta_w: min_delta_w = delta_w            
     
-    #pass
+    pass
 
 ## Check that extreme values in rating curve have been exceeded