Changeset 7035 for branches/numpy/anuga/culvert_flows/culvert_class.py

Timestamp:

May 14, 2009, 2:24:54 PM (15 years ago)

Author:

rwilson

Message:

Back-merge from Numeric trunk.

File:

• branches/numpy/anuga/culvert_flows/culvert_class.py (modified) (6 diffs)

branches/numpy/anuga/culvert_flows/culvert_class.py

@@ -15 +15 @@
 
 import numpy as num
-
+from math import sqrt
 
 class Below_interval(Exception): pass 
@@ -137 +137 @@
                  manning=None,          # Mannings Roughness for Culvert 
                  sum_loss=None,
-                 use_velocity_head=False, 
+                 use_velocity_head=False, # FIXME(Ole): Get rid of - always True
                  use_momentum_jet=False, # FIXME(Ole): Not yet implemented
                  label=None,
@@ -464 +464 @@
             log_filename = self.log_filename
             
-        # Compute stage and energy at the 
+        # Compute stage, energy and velocity at the 
         # enquiry points at each end of the culvert
         openings = self.openings
@@ -474 +474 @@
             depth = h = stage-opening.elevation
                                                            
-                                               
+            
+            # Get velocity                                 
+            xmomentum = dq['xmomentum'].get_values(location='centroids',
+                                                   indices=[idx])[0]
+            ymomentum = dq['xmomentum'].get_values(location='centroids',
+                                                   indices=[idx])[0]
+            
+            if h > minimum_allowed_height:
+                u = xmomentum/(h + velocity_protection/h)
+                v = ymomentum/(h + velocity_protection/h)
+            else:
+                u = v = 0.0
+                
+            v_squared = u*u + v*v
+            
             if self.use_velocity_head is True:
-                xmomentum = dq['xmomentum'].get_values(location='centroids',
-                                                       indices=[idx])[0]
-                ymomentum = dq['xmomentum'].get_values(location='centroids',
-                                                       indices=[idx])[0]
-            
-                if h > minimum_allowed_height:
-                    u = xmomentum/(h + velocity_protection/h)
-                    v = ymomentum/(h + velocity_protection/h)
-                else:
-                    u = v = 0.0
-                
-                velocity_head = 0.5*(u*u + v*v)/g    
+                velocity_head = 0.5*v_squared/g    
             else:
                 velocity_head = 0.0
@@ -495 +498 @@
             opening.stage = stage
             opening.depth = depth
+            opening.velocity = sqrt(v_squared)
             
 
@@ -583 +587 @@
                     self.culvert_routine(inlet.depth,
                                          outlet.depth,
+                                         inlet.velocity,
+                                         outlet.velocity,
                                          inlet.specific_energy, 
                                          delta_total_energy,