Changeset 6079

Timestamp:

Dec 15, 2008, 7:15:59 PM (16 years ago)

Author:

ole

Message:

Fixed dry bed culvert behaviour

Location:

anuga_core/source/anuga/culvert_flows

Files:

• culvert_class.py (modified) (2 diffs)
• test_culvert_class.py (modified) (3 diffs)

anuga_core/source/anuga/culvert_flows/culvert_class.py

@@ -308 +308 @@
                     
                 # Store current average stage and depth with each opening object
+                opening.depth = stage - opening.elevation
                 opening.stage = stage
 
@@ -327 +328 @@
 
 
-            # Calculate discharge for one barrel and set inlet.rate and outlet.rate
-            try:
-                Q = interpolate_linearly(delta_w, self.rating_curve[:,0], self.rating_curve[:,1]) 
-            except Below_interval, e:
-                Q = self.rating_curve[0,1]             
-                msg = '%.2fs: Delta head smaller than rating curve minimum: ' %time
-                msg += str(e)
-                msg += '\n        I will use minimum discharge %.2f m^3/s for culvert "%s"'\
-                    %(Q, self.label)
-                if hasattr(self, 'log_filename'):                    
-                    log_to_file(self.log_filename, msg)
-            except Above_interval, e:
-                Q = self.rating_curve[-1,1]             
-                msg = '%.2fs: Delta head greater than rating curve maximum: ' %time
-                msg += str(e)
-                msg += '\n        I will use maximum discharge %.2f m^3/s for culvert "%s"'\
-                    %(Q, self.label)
-                if hasattr(self, 'log_filename'):                    
-                    log_to_file(self.log_filename, msg)                    
+            if self.inlet.depth <= 0.01:
+                Q = 0.0
+            else:
+                # Calculate discharge for one barrel and set inlet.rate and outlet.rate
+                try:
+                    Q = interpolate_linearly(delta_w, self.rating_curve[:,0], self.rating_curve[:,1]) 
+                except Below_interval, e:
+                    Q = self.rating_curve[0,1]             
+                    msg = '%.2fs: Delta head smaller than rating curve minimum: ' %time
+                    msg += str(e)
+                    msg += '\n        I will use minimum discharge %.2f m^3/s for culvert "%s"'\
+                        %(Q, self.label)
+                    if hasattr(self, 'log_filename'):                    
+                        log_to_file(self.log_filename, msg)
+                except Above_interval, e:
+                    Q = self.rating_curve[-1,1]             
+                    msg = '%.2fs: Delta head greater than rating curve maximum: ' %time
+                    msg += str(e)
+                    msg += '\n        I will use maximum discharge %.2f m^3/s for culvert "%s"'\
+                        %(Q, self.label)
+                    if hasattr(self, 'log_filename'):                    
+                        log_to_file(self.log_filename, msg)                    
 
                 

anuga_core/source/anuga/culvert_flows/test_culvert_class.py

@@ -12 +12 @@
     Transmissive_boundary, Time_boundary
 
-from anuga.culvert_flows.culvert_class import Culvert_flow_rating
+from anuga.culvert_flows.culvert_class import Culvert_flow_rating, Culvert_flow_energy
 from anuga.culvert_flows.culvert_routines import boyd_generalised_culvert_model
      
@@ -124 +124 @@
 
 
-    def test_that_culvert_dry_bed(self):
+    def test_that_culvert_dry_bed_rating(self):
         """test_that_culvert_in_dry_bed_does_not_produce_flow(self):
         
         Test that culvert on a sloping dry bed doesn't produce flows
         although there will be a 'pressure' head due to delta_w > 0
+
+        This one is using the rating curve variant
         """
 
@@ -214 +216 @@
         ref_volume = domain.get_quantity('stage').get_integral()
         for t in domain.evolve(yieldstep = 1, finaltime = 25):
-            print domain.timestepping_statistics()
+            #print domain.timestepping_statistics()
+            new_volume = domain.get_quantity('stage').get_integral()
+            
+            msg = 'Total volume has changed'
+            assert allclose(new_volume, ref_volume), msg
+            pass
+    
+
+
+    def test_that_culvert_dry_bed_boyd(self):
+        """test_that_culvert_in_dry_bed_does_not_produce_flow(self):
+        
+        Test that culvert on a sloping dry bed doesn't produce flows
+        although there will be a 'pressure' head due to delta_w > 0
+
+        This one is using the 'Boyd' variant        
+        """
+
+        path = get_pathname_from_package('anuga.culvert_flows')    
+        
+        length = 40.
+        width = 5.
+
+        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('Test_culvert_dry') # Output name
+        domain.set_default_order(2)
+
+
+        #----------------------------------------------------------------------
+        # Setup initial conditions
+        #----------------------------------------------------------------------
+
+        def topography(x, y):
+            """Set up a weir
+            
+            A culvert will connect either side
+            """
+            # General Slope of Topography
+            z=-x/1000
+            
+            N = len(x)
+            for i in range(N):
+
+               # Sloping Embankment Across Channel
+                if 5.0 < x[i] < 10.1:
+                    # Cut Out Segment for Culvert FACE                
+                    if  1.0+(x[i]-5.0)/5.0 <  y[i]  < 4.0 - (x[i]-5.0)/5.0: 
+                       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:
+                    # Cut Out Segment for Culvert FACE                
+                    if  2.0-(x[i]-12.0)/2.5 <  y[i]  < 3.0 + (x[i]-12.0)/2.5:
+                       z[i]=z[i]
+                    else:
+                       z[i] +=  2.5-1.0*(x[i] -12.0) # Sloping D/S Face
+                                   
+                        
+            return z
+
+
+        domain.set_quantity('elevation', topography) 
+        domain.set_quantity('friction', 0.01)         # Constant friction 
+        domain.set_quantity('stage',
+                            expression='elevation')   # Dry initial condition
+
+
+        filename = os.path.join(path, 'example_rating_curve.csv')
+
+
+        culvert = Culvert_flow_energy(domain,
+                                      label='Culvert No. 1',
+                                      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],
+                                      width=1.20,height=0.75,
+                                      culvert_routine=boyd_generalised_culvert_model,        
+                                      number_of_barrels=1,
+                                      update_interval=2,
+                                      verbose=True)
+        
+        domain.forcing_terms.append(culvert)
+        
+
+        #-----------------------------------------------------------------------
+        # Setup boundary conditions
+        #-----------------------------------------------------------------------
+
+        # Inflow based on Flow Depth and Approaching Momentum
+
+        Br = Reflective_boundary(domain)              # Solid reflective wall
+        domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
+
+
+        #-----------------------------------------------------------------------
+        # Evolve system through time
+        #-----------------------------------------------------------------------
+
+        ref_volume = domain.get_quantity('stage').get_integral()
+        for t in domain.evolve(yieldstep = 1, finaltime = 25):
+            #print domain.timestepping_statistics()
             new_volume = domain.get_quantity('stage').get_integral()