Changeset 5340

Timestamp:

May 16, 2008, 5:12:37 PM (17 years ago)

Author:

ole

Message:

Work on momentum jet

File:

• anuga_work/development/Rudy_vandrie_2007/culvert_development.py (modified) (6 diffs)

anuga_work/development/Rudy_vandrie_2007/culvert_development.py

@@ -286 +286 @@
         openings = self.openings
         for i, opening in enumerate(openings):
-            stage = domain.quantities['stage'].get_values(location='centroids',
+            # Quantities corresponding to fluid exchange field for this opening            
+            stage_o = domain.quantities['stage'].get_values(location='centroids',
                                                           indices=opening.indices)
-            elevation = domain.quantities['elevation'].get_values(location='centroids',
+            elevation_o = domain.quantities['elevation'].get_values(location='centroids',
                                                                   indices=opening.indices)
+            xmomentum_o = domain.quantities['xmomentum'].get_values(location='centroids',
+                                                                  indices=opening.indices)
+            ymomentum_o = domain.quantities['xmomentum'].get_values(location='centroids',
+                                                                  indices=opening.indices)
+
+            # Compute mean velocity in the exchange area in front of the culvert (taking zero depths into account)
+            depth_o = stage_o - elevation_o
+            
+            ux_o = xmomentum_o/(depth_o+epsilon)
+            uy_o = ymomentum_o/(depth_o+epsilon)
+            v_o = mean(sqrt(ux_o**2+uy_o**2))
+            d_o = mean(depth_o)
+            w_o = mean(stage_o)
+            z_o = mean(elevation_o)
+            self.mean_xmomentum_o = mean(xmomentum_o)
+            self.mean_ymomentum_o = mean(ymomentum_o)             
 
             # Indices corresponding to energy enquiry field for this opening
@@ -317 +334 @@
             w = mean(stage)
             z = mean(elevation)
-            
+            self.mean_xmomentum = mean(xmomentum)
+            self.mean_ymomentum = mean(ymomentum)             
+
             # Ratio of depth to culvert height
             ratio = d/(2*self.height)
@@ -436 +455 @@
                 delta_t = time - self.last_time
                 if delta_t > 0.0:
-                    xmomentum_rate = outlet_mom_x - self.xmom_forcing1.value
-                    xmomentum_rate /= delta_t
+                    xmomentum_rate = outlet_mom_x - self.mean_xmomentum
+                    if xmomentum_rate > 0:
+                        xmomentum_rate /= delta_t
+                    else:
+                        xmomentum_rate = 0.0
                     
-                    ymomentum_rate = outlet_mom_y - self.ymom_forcing1.value
-                    ymomentum_rate /= delta_t                        
+                    ymomentum_rate = outlet_mom_y - self.mean_ymomentum
+                    if ymomentum_rate > 0:
+                        ymomentum_rate /= delta_t
+                    else:
+                        ymomentum_rate = 0.0                    
                 else:
                     xmomentum_rate = ymomentum_rate = 0.0
@@ -453 +478 @@
 
 
-                if int(domain.time*100) % 100 == 0:
-                    s = 'T=%.2f, Culvert Discharge = %.3f  Culv. Vel. %0.3f'\
+                if int(time*100) % 100 == 0:
+                    s = 'T=%.3f, Culvert Discharge = %.3f  Culv. Vel. %0.3f'\
                         %(time, flow_rate_control, outlet_vel)
-                    s += ' Depth= %0.3f  Momentum = (%0.3f, %0.3f)'\
-                         %(outlet_dep, outlet_mom_x,outlet_mom_y)
-                    s += ' Momentum rate: (%.4f, %.4f)'\
+                    s += ' Depth= %0.3f\n    Outlet Momentum = (%0.3f, %0.3f)\n'\
+                         %(outlet_dep, outlet_mom_x, outlet_mom_y)
+                    s += '    Avg Momentum at opening = (%0.3f, %0.3f)\n'\
+                         %(self.mean_xmomentum_o, self.mean_ymomentum_o)
+                    s += '    Avg Momentum in enquiry = (%0.3f, %0.3f)\n'\
+                         %(self.mean_xmomentum, self.mean_ymomentum)                    
+                    s += '    Momentum rate: (%.4f, %.4f)'\
                          %(xmomentum_rate, ymomentum_rate)                    
                     fid.write(s)
@@ -488 +517 @@
         # Print out flows every 1 seconds
         if int(time*100) % 100 == 0:
-            s = 'Time %.2f\n' %time
+            s = 'Time %.3f\n' %time
             s += '    Opening[0]: d=%.3f, A=%.3f, E=%.3f, r=%.3f\n'\
                  %(openings[0].depth,
@@ -506 +535 @@
             
             print s
+
             fid.write(s)