Changeset 5736

Timestamp:

Sep 5, 2008, 11:30:19 AM (16 years ago)

Author:

ole

Message:

Implemented run-time version of get_energy_through_cross_section and tested it.

Location:

anuga_core/source/anuga

Files:

• abstract_2d_finite_volumes/quantity.py (modified) (2 diffs)
• coordinate_transforms/geo_reference.py (modified) (1 diff)
• shallow_water/data_manager.py (modified) (3 diffs)
• shallow_water/shallow_water_domain.py (modified) (15 diffs)
• shallow_water/test_shallow_water_domain.py (modified) (1 diff)

anuga_core/source/anuga/abstract_2d_finite_volumes/quantity.py

@@ -1012 +1012 @@
         associated matrices will be reused to save time.
         
-        The argument interpolation points must be given as either a list of absolute UTM coordinates or
-        a geospatial data object.
-        """
-
-        # FIXME (Ole): Could do with an input check (should be generalised and used widely)
-        # That will check that interpolation points is either a list of points, Nx2 array, or geospatial
+        The argument interpolation points must be given as either a 
+        list of absolute UTM coordinates or a geospatial data object.
+        """
+
+        # FIXME (Ole): Could do with an input check (should be generalised 
+        # and used widely)
+        # That will check that interpolation points is either a list of 
+        # points, Nx2 array, or geospatial
         
         # Ensure points are converted to coordinates relative to mesh origin
-        # FIXME (Ole): This could all be refactored using the 'change_points_geo_ref' method 
-        # of Class geo_reference. The purpose is to make interpolation points relative
+        # FIXME (Ole): This could all be refactored using the 
+        # 'change_points_geo_ref' method of Class geo_reference. 
+        # The purpose is to make interpolation points relative
         # to the mesh origin.
         #
@@ -1034 +1037 @@
         interpolation_points = ensure_numeric(interpolation_points)
         
-        # Get internal (discontinuous) triangles for use with the interpolation object.
+        
+        
+        # Get internal (discontinuous) triangles for use with the 
+        # interpolation object.
         x, y, vertex_values, triangles = self.get_vertex_values(xy=True,
                                                                 smooth=False)

anuga_core/source/anuga/coordinate_transforms/geo_reference.py

@@ -212 +212 @@
         assert points.shape[1] == 2, msg                
 
-            
+        
+        # FIXME (Ole): Could also check if zone, xllcorner, yllcorner 
+        # are identical in the two geo refs.    
         if points_geo_ref is not self:
-            #add point geo ref to points
+            # If georeferences are different
+        
             if not points_geo_ref is None:
+                # Convert points to absolute coordinates
                 points[:,0] += points_geo_ref.xllcorner 
                 points[:,1] += points_geo_ref.yllcorner 
         
-            #subtract primary geo ref from points
+            # Make points relative to primary geo reference
             points[:,0] -= self.xllcorner 
             points[:,1] -= self.yllcorner

anuga_core/source/anuga/shallow_water/data_manager.py

@@ -6466 +6466 @@
                                      kind = 'total',
                                      verbose=False):
-    """Obtain flow (m^3/s) perpendicular to specified cross section.
+    """Obtain average energy head [m] across 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 
@@ -6478 +6477 @@
 
     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].
+        E: Average energy [m] across given segments for all stored times.
+
+    The average velocity is computed for each triangle intersected by 
+    the polyline and averaged weighted 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].
     """
 
@@ -6501 +6500 @@
     # 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,
+                                                          quantity_names=\
+                                                          quantity_names,
                                                           polyline=polyline,
                                                           verbose=verbose)    

anuga_core/source/anuga/shallow_water/shallow_water_domain.py

@@ -384 +384 @@
         """Get the total flow through an arbitrary poly line.        
         
-        This is a run-time equivalent of the function with same name in data_manager.py
+        This is a run-time equivalent of the function with same name 
+        in data_manager.py
         
         Input:
@@ -410 +411 @@
         midpoints = segment_midpoints(segments)       
         
-        # FIXME (Ole): HACK - need to make midpoints Geospatial instances
-        #midpoints = self.geo_reference.get_absolute(midpoints)        
-        
         # Make midpoints Geospatial instances
         midpoints = ensure_geospatial(midpoints, self.geo_reference)        
@@ -439 +437 @@
             # Accumulate
             total_flow += segment_flow
-
             
         return total_flow
-        
+
+        
+               
+    def get_energy_through_cross_section(self, polyline,
+                                         kind = 'total',
+                                         verbose=False):               
+        """Obtain average energy head [m] across specified cross section.
+
+        Inputs:
+            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:
+            E: Average energy [m] across given segments for all stored times.
+
+        The average velocity is computed for each triangle intersected by 
+        the polyline and averaged weighted 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        
+                                         
+        
+        # Adjust polyline to mesh spatial origin
+        polyline = self.geo_reference.get_relative(polyline)
+        
+        # Find all intersections and associated triangles.
+        segments = self.get_intersecting_segments(polyline, verbose=verbose)
+
+        msg = 'No segments found'
+        assert len(segments) > 0, msg
+        
+        # Get midpoints
+        midpoints = segment_midpoints(segments)       
+        
+        # Make midpoints Geospatial instances
+        midpoints = ensure_geospatial(midpoints, self.geo_reference)        
+        
+        # Compute energy        
+        if verbose: print 'Computing %s energy' %kind        
+        
+        # Get interpolated values
+        stage = self.get_quantity('stage')        
+        elevation = self.get_quantity('elevation')                
+        xmomentum = self.get_quantity('xmomentum')
+        ymomentum = self.get_quantity('ymomentum')        
+
+        w = stage.get_values(interpolation_points=midpoints)
+        z = elevation.get_values(interpolation_points=midpoints)        
+        uh = xmomentum.get_values(interpolation_points=midpoints)
+        vh = ymomentum.get_values(interpolation_points=midpoints)        
+        h = w-z # Depth
+        
+        # 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 and sum flows across each segment
+        average_energy=0.0
+        for i in range(len(w)):
+            
+            # Average velocity across this segment
+            if h[i] > epsilon:
+                # Use protection against degenerate velocities
+                u = uh[i]/(h[i] + h0/h[i])
+                v = vh[i]/(h[i] + h0/h[i])
+            else:
+                u = v = 0.0
+                
+            speed_squared = u*u + v*v    
+            kinetic_energy = 0.5*speed_squared/g
+            
+            if kind == 'specific':
+                segment_energy = h[i] + kinetic_energy
+            elif kind == 'total':
+                segment_energy = w[i] + 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
+             
+            
+        return average_energy
+        
+
+                        
 
     def check_integrity(self):
@@ -1533 +1630 @@
 
 class General_forcing:
-    """Class General_forcing - general explicit forcing term for update of quantity
+    """General explicit forcing term for update of quantity
     
     This is used by Inflow and Rainfall for instance
@@ -1541 +1638 @@
 
     domain:     ANUGA computational domain
-    quantity_name: Name of quantity to update. It must be a known conserved quantity.
+    quantity_name: Name of quantity to update. 
+                   It must be a known conserved quantity.
+                   
     rate [?/s]: Total rate of change over the specified area.  
                 This parameter can be either a constant or a
@@ -1625 +1724 @@
 
             for point in periphery_points:
-                msg = 'Point %s on periphery for forcing term did not' %(str(point))
-                msg += ' fall within the domain boundary.'
+                msg = 'Point %s on periphery for forcing term' %(str(point))
+                msg += ' did not fall within the domain boundary.'
                 assert is_inside_polygon(point, bounding_polygon), msg
 
@@ -1635 +1734 @@
             # Check that polygon lies within the mesh.
             for point in self.polygon:
-                msg = 'Point %s in polygon for forcing term did not' %(str(point))
-                msg += 'fall within the domain boundary.'
+                msg = 'Point %s in polygon for forcing term' %(str(point))
+                msg += ' did not fall within the domain boundary.'
                 assert is_inside_polygon(point, bounding_polygon), msg
-        
-
-
-            
 
 
         # Pointer to update vector
-        self.update = domain.quantities[self.quantity_name].explicit_update            
+        self.update = domain.quantities[self.quantity_name].explicit_update
 
         # Determine indices in flow area
@@ -1661 +1756 @@
             for k in range(N):
                 x, y = points[k,:] # Centroid
-                if ((x-self.center[0])**2+(y-self.center[1])**2) < self.radius**2:
+                
+                c = self.center
+                if ((x-c[0])**2+(y-c[1])**2) < self.radius**2:
                     self.exchange_indices.append(k)
                     
@@ -1676 +1773 @@
         
             if len(self.exchange_indices) == 0:
-                msg = 'No triangles have been identified in specified region: %s' %inlet_region
+                msg = 'No triangles have been identified in '
+                msg += 'specified region: %s' %inlet_region
                 raise Exception, msg
 
@@ -1723 +1821 @@
         """Return values for specified quantity restricted to opening 
         """
-        return self.domain.quantities[self.quantity_name].get_values(indices=self.exchange_indices)
+        
+        q = self.domain.quantities[self.quantity_name]
+        return q.get_values(indices=self.exchange_indices)
     
 
@@ -1729 +1829 @@
         """Set values for specified quantity restricted to opening 
         """
-        self.domain.quantities[self.quantity_name].set_values(val, indices=self.exchange_indices)    
+
+        q = self.domain.quantities[self.quantity_name]                
+        q.set_values(val, indices=self.exchange_indices)    
 
 
@@ -1766 +1868 @@
     How to put them in a run File...
         
-    #--------------------------------------------------------------------------
+    #------------------------------------------------------------------------
     # Setup specialised forcing terms
-    #--------------------------------------------------------------------------
+    #------------------------------------------------------------------------
     # This is the new element implemented by Ole and Rudy to allow direct
     # input of Inflow in mm/s
@@ -1842 +1944 @@
 
 
-    #--------------------------------------------------------------------------
+    #------------------------------------------------------------------------
     # Setup specialised forcing terms
-    #--------------------------------------------------------------------------
+    #------------------------------------------------------------------------
     # This is the new element implemented by Ole to allow direct input
     # of Inflow in m^3/s
@@ -1862 +1964 @@
                  polygon=None,
                  verbose=False):                 
-
-
-        #msg = 'Class Inflow must have either center & radius or a polygon specified.'
-        #assert center is not None and radius is not None or\
-        #       polygon is not None, msg
 
 
@@ -1882 +1979 @@
         overriding version in descendant
 
-        This one converts m^3/s to m/s which can be added directly to 'stage' in ANUGA
+        This one converts m^3/s to m/s which can be added directly 
+        to 'stage' in ANUGA
         """
 

anuga_core/source/anuga/shallow_water/test_shallow_water_domain.py

@@ -1502 +1502 @@
                 assert allclose(Q, uh*width)
 
+
+        
+    def test_get_energy_through_cross_section_with_geo(self):
+        """test_get_energy_through_cross_section(self):
+
+        Test that the total and specific energy through a cross section can be
+        correctly obtained at run-time from the ANUGA domain.
+        
+        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
+        e = -1 m
+        u = 2 m/s
+        h = 2 m
+        w = 3 m (width of channel)
+
+        q = u*h*w = 12 m^3/s
+
+
+        This run tries it with georeferencing and with elevation = -1
+        
+        """
+
+        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_quantities_to_be_stored(None)
+
+
+        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 3 m inlet: 
+
+
+        # Initial conditions
+        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})
+        
+        
+        # Interpolation points down the middle
+        I = [[0, width/2.],
+             [length/2., width/2.],
+             [length, width/2.]]
+        interpolation_points = domain.geo_reference.get_absolute(I)        
+        
+        # Shortcuts to quantites
+        stage = domain.get_quantity('stage')        
+        xmomentum = domain.get_quantity('xmomentum')       
+        ymomentum = domain.get_quantity('ymomentum')               
+
+        for t in domain.evolve(yieldstep=0.1, finaltime = t_end):
+            # Check that quantities are they should be in the interior
+
+            w_t = stage.get_values(interpolation_points)            
+            uh_t = xmomentum.get_values(interpolation_points)
+            vh_t = ymomentum.get_values(interpolation_points)            
+            
+            assert allclose(w_t, w)
+            assert allclose(uh_t, uh)            
+            assert allclose(vh_t, 0.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)    
+                Es = domain.get_energy_through_cross_section(cross_section,
+                                                             kind='specific',
+                                                             verbose=False)
+                                                      
+                assert allclose(Es, h + 0.5*u*u/g)
+            
+                Et = domain.get_energy_through_cross_section(cross_section,
+                                                             kind='total',
+                                                             verbose=False)
+                assert allclose(Et, w + 0.5*u*u/g)