Changeset 5726

Timestamp:

Sep 3, 2008, 4:44:23 PM (16 years ago)

Author:

ole

Message:

Energy through cross section and simple tests done for sww file input.
This is work towards ticket:295

Location:

anuga_core/source/anuga/shallow_water

Files:

• data_manager.py (modified) (1 diff)
• test_data_manager.py (modified) (5 diffs)

anuga_core/source/anuga/shallow_water/data_manager.py

@@ -6466 +6466 @@
     
 
+def get_energy_through_cross_section(filename,
+                                     polyline,
+                                     kind = 'total',
+                                     verbose=False):
+    """Obtain flow (m^3/s) perpendicular to specified cross section.
+
+    Inputs:
+        filename: Name of sww file
+        polyline: Representation of desired cross section - it may contain 
+                  multiple sections allowing for complex shapes. Assume 
+                  absolute UTM coordinates.
+                  Format [[x0, y0], [x1, y1], ...]
+        kind:     Select which energy to compute. 
+                  Options are 'specific' and 'total' (default)
+
+    Output:
+        time: All stored times in sww file
+        Q: Average energy [m] across given segments for all stored times.
+
+    The average velocity is computed for each triangle intersected by the polyline 
+    and averaged weigted by segment lengths. 
+
+    The typical usage of this function would be to get average energy of flow in a channel,
+    and the polyline would then be a cross section perpendicular to the flow.
+
+    
+    #FIXME (Ole) - need name for this energy reflecting that its dimension is [m].
+    """
+
+    from anuga.config import g, epsilon, velocity_protection as h0        
+            
+    quantity_names =['elevation',
+                     'stage',
+                     'xmomentum',
+                     'ymomentum']
+
+
+    # Get values for quantities at each midpoint of poly line from sww file
+    X = get_interpolated_quantities_at_polyline_midpoints(filename,
+                                                          quantity_names=quantity_names,
+                                                          polyline=polyline,
+                                                          verbose=verbose)    
+    segments, interpolation_function = X
+
+    
+    # Get vectors for time and interpolation_points
+    time = interpolation_function.time
+    interpolation_points = interpolation_function.interpolation_points    
+
+    if verbose: print 'Computing %s energy' %kind
+    
+    # Compute total length of polyline for use with weighted averages
+    total_line_length = 0.0
+    for segment in segments:
+        total_line_length += segment.length
+        
+    # Compute energy
+    E = []
+    for t in time:
+        average_energy=0.0
+        for i, p in enumerate(interpolation_points):
+            elevation, stage, uh, vh = interpolation_function(t, point_id=i)
+            
+            # Depth
+            h = depth = stage-elevation
+            
+            # Average velocity across this segment
+            if h > epsilon:
+                # Use protection against degenerate velocities
+                u = uh/(h + h0/h)
+                v = vh/(h + h0/h)
+            else:
+                u = v = 0.0
+                
+            speed_squared = u*u + v*v    
+            kinetic_energy = 0.5*speed_squared/g
+            
+            if kind == 'specific':
+                segment_energy = depth + kinetic_energy
+            elif kind == 'total':
+                segment_energy = stage + kinetic_energy                
+            else:
+                msg = 'Energy kind must be either "specific" or "total".'
+                msg += ' I got %s' %kind
+                
+
+            # Add to weighted average
+            weigth = segments[i].length/total_line_length
+            average_energy += segment_energy*weigth
+             
+
+        # Store energy at this timestep    
+        E.append(average_energy)
+
+
+    return time, E
+    
+
 
 

anuga_core/source/anuga/shallow_water/test_data_manager.py

@@ -17 +17 @@
 from anuga.shallow_water import *
 from anuga.shallow_water.data_manager import *
-from anuga.config import epsilon
+from anuga.config import epsilon, g
 from anuga.utilities.anuga_exceptions import ANUGAError
 from anuga.utilities.numerical_tools import ensure_numeric
@@ -10198 +10198 @@
         The specifics are
         u = 2 m/s
-        h = 1 m
+        h = 2 m
         w = 5 m (width of channel)
 
-        q = u*h*w = 10 m^3/s
-
-
-        This run tries it with georeferencing
+        q = u*h*w = 20 m^3/s
+
+
+        This run tries it with georeferencing and with elevation = -1
         
         """
@@ -10236 +10236 @@
         domain.smooth = True
 
-        h = 1.0
+        e = -1.0
+        w = 1.0
+        h = w-e
         u = 2.0
         uh = u*h
 
         Br = Reflective_boundary(domain)     # Side walls
-        Bd = Dirichlet_boundary([h, uh, 0])  # 2 m/s across the 5 m inlet: 
-
-
-
-        
-        domain.set_quantity('elevation', 0.0)
-        domain.set_quantity('stage', h)
+        Bd = Dirichlet_boundary([w, uh, 0])  # 2 m/s across the 5 m inlet: 
+
+
+
+        
+        domain.set_quantity('elevation', e)
+        domain.set_quantity('stage', w)
         domain.set_quantity('xmomentum', uh)
         domain.set_boundary( {'left': Bd, 'right': Bd, 'top': Br, 'bottom': Br})
@@ -10268 +10270 @@
         for t in range(t_end+1):
             for i in range(3):
-                assert allclose(f(t, i), [1, 2, 0])
+                #print i, t, f(t, i)            
+                assert allclose(f(t, i), [w, uh, 0])
             
 
@@ -10284 +10287 @@
 
 
+            
+    def test_get_energy_through_cross_section(self):
+        """test_get_energy_through_cross_section(self):
+
+        Test that the specific and total energy through a cross section can be
+        correctly obtained from an sww file.
+        
+        This test creates a flat bed with a known flow through it and tests
+        that the function correctly returns the expected energies.
+
+        The specifics are
+        u = 2 m/s
+        h = 1 m
+        w = 5 m (width of channel)
+
+        q = u*h*w = 10 m^3/s
+        Es = h + 0.5*v*v/g  # Specific energy head [m]
+        Et = w + 0.5*v*v/g  # Total energy head [m]        
+
+
+        This test uses georeferencing
+        
+        """
+
+        import time, os
+        from Numeric import array, zeros, allclose, Float, concatenate
+        from Scientific.IO.NetCDF import NetCDFFile
+
+        # Setup
+        from mesh_factory import rectangular
+
+        # Create basic mesh (20m x 3m)
+        width = 3
+        length = 20
+        t_end = 1
+        points, vertices, boundary = rectangular(length, width,
+                                                 length, width)
+
+        # Create shallow water domain
+        domain = Domain(points, vertices, boundary,
+                        geo_reference = Geo_reference(56,308500,6189000))
+
+        domain.default_order = 2
+        domain.set_minimum_storable_height(0.01)
+
+        domain.set_name('flowtest')
+        swwfile = domain.get_name() + '.sww'
+
+        domain.set_datadir('.')
+        domain.format = 'sww'
+        domain.smooth = True
+
+        e = -1.0
+        w = 1.0
+        h = w-e
+        u = 2.0
+        uh = u*h
+
+        Br = Reflective_boundary(domain)     # Side walls
+        Bd = Dirichlet_boundary([w, uh, 0])  # 2 m/s across the 5 m inlet: 
+
+        
+        domain.set_quantity('elevation', e)
+        domain.set_quantity('stage', w)
+        domain.set_quantity('xmomentum', uh)
+        domain.set_boundary( {'left': Bd, 'right': Bd, 'top': Br, 'bottom': Br})
+
+        for t in domain.evolve(yieldstep=1, finaltime = t_end):
+            pass
+
+        # Check that momentum is as it should be in the interior
+
+        I = [[0, width/2.],
+             [length/2., width/2.],
+             [length, width/2.]]
+        
+        I = domain.geo_reference.get_absolute(I)
+        f = file_function(swwfile,
+                          quantities=['stage', 'xmomentum', 'ymomentum'],
+                          interpolation_points=I,
+                          verbose=False)
+
+        for t in range(t_end+1):
+            for i in range(3):
+                #print i, t, f(t, i)
+                assert allclose(f(t, i), [w, uh, 0])
+            
+
+        # Check energies through the middle
+        for i in range(5):
+            x = length/2. + i*0.23674563 # Arbitrary
+            cross_section = [[x, 0], [x, width]]
+
+            cross_section = domain.geo_reference.get_absolute(cross_section)            
+            
+            time, Es = get_energy_through_cross_section(swwfile,
+                                                       cross_section,
+                                                       kind='specific',
+                                                       verbose=False)
+            assert allclose(Es, h + 0.5*u*u/g)
+            
+            time, Et = get_energy_through_cross_section(swwfile,
+                                                        cross_section,
+                                                        kind='total',
+                                                        verbose=False)
+            assert allclose(Et, w + 0.5*u*u/g)            
+
+            
         
     def test_get_all_swwfiles(self):