Changeset 6143

Timestamp:

Jan 13, 2009, 11:12:52 AM (16 years ago)

Author:

ole

Message:

Simplified flow protection at the inlet and wrote better diagnostics.

Location:

anuga_core/source/anuga/culvert_flows

Files:

• Test_Culvert_Flat_Water_Lev.py (modified) (3 diffs)
• culvert_class.py (modified) (4 diffs)

anuga_core/source/anuga/culvert_flows/Test_Culvert_Flat_Water_Lev.py

@@ -21 +21 @@
      Transmissive_boundary, Time_boundary
 
-from anuga.culvert_flows.culvert_class import Culvert_flow_rating, Culvert_flow_energy
+from anuga.culvert_flows.culvert_class import Culvert_flow
 from anuga.culvert_flows.culvert_routines import boyd_generalised_culvert_model
      
@@ -104 +104 @@
 
 
-culvert_rating = Culvert_flow_rating(domain,
+culvert_rating = Culvert_flow(domain,
                        culvert_description_filename='example_rating_curve.csv',
                        end_point0=[9.0, 2.5], 
@@ -111 +111 @@
 
 
-culvert_energy = Culvert_flow_energy(domain,
+culvert_energy = Culvert_flow(domain,
                        label='Culvert No. 1',
                        description='This culvert is a test unit 1.2m Wide by 0.75m High',   

anuga_core/source/anuga/culvert_flows/culvert_class.py

@@ -26 +26 @@
 class Above_interval(Exception): pass
 
+# FIXME(Ole): Take a good hard look at logging here
 
 
@@ -404 +405 @@
         """Adjust Q downwards depending on available water at inlet
         """
+    
+        if delta_t < epsilon:
+            # No need to adjust if time step is very small or zero
+            # In this case the possible flow will be very large
+            # anyway.
+            return Q
         
         # Short hands
@@ -410 +417 @@
         time = domain.get_time()
         I = self.inlet
-
-        
-        # Find triangle with the smallest available flow
-        max_Q = sys.maxint
-        if delta_t > 0:
-            for i in I.exchange_indices:
-                stage = dq['stage'].get_values(location='centroids', 
-                                               indices=[i])[0]
-                elevation = dq['elevation'].get_values(location='centroids', 
-                                                       indices=[i])[0]        
-                depth = stage-elevation
-                area = domain.areas[i]
-
-                # Possible rate based on this triangle
-                Q_possible = depth*I.exchange_area/delta_t
-                                
-                # Use triangle with least depth to determine total flow
-                max_Q = min(max_Q, Q_possible)
-
-            
+        idx = I.exchange_indices    
+
+        # Find triangle with the smallest depth
+        stage = dq['stage'].get_values(location='centroids', 
+                                               indices=[idx])
+        elevation = dq['elevation'].get_values(location='centroids', 
+                                               indices=[idx])        
+        depth = stage-elevation
+        min_depth = min(depth.flat)
+
+        # Compute possible flow for exchange region based on
+        # triangle with smallest depth
+        max_Q = min_depth*I.exchange_area/delta_t        
+        
+        # Calculate the minimum in absolute terms of
+        # the requsted flow and the possible flow
         Q_reduced = sign(Q)*min(abs(Q), abs(max_Q))
         
         if abs(Q_reduced) < abs(Q): 
-            msg = '%.2fs: Computed extraction for this time interval (Q*dt) is ' % time
-            msg += 'greater than current volume (V) at inlet.\n'
-            msg += ' Q will be reduced from %.2f m^3/s to %.2f m^3/s.' % (Q, Q_reduced)
+            msg = '%.2fs: Requested flow is ' % time
+            msg += 'greater than what is supported by the smallest '
+            msg += 'depth at inlet exchange area:\n        '
+            msg += 'h_min*inlet_area/delta_t = %.2f*%.2f/%.2f '\
+                % (min_depth, I.exchange_area, delta_t)
+            msg += ' = %.2f m^3/s\n        ' % Q_reduced
+            msg += 'Q will be reduced from %.2f m^3/s to %.2f m^3/s.' % (Q, Q_reduced)
             if self.verbose is True:
                 print msg
@@ -521 +528 @@
             # Flow will be purely controlled by uphill outlet face
             if self.verbose is True:
-                print '%.2f - WARNING: Flow is running uphill.' % time
+                print '%.2fs - WARNING: Flow is running uphill.' % time
             
         if hasattr(self, 'log_filename'):