Changeset 5288

Timestamp:

May 7, 2008, 4:01:52 PM (17 years ago)

Author:

ole

Message:

Implemeted and tested ticket:175 - flow_through_cross_section.
Also updated manual

Location:

anuga_core

Files:

• documentation/release_information/release_notes.txt (modified) (1 diff)
• documentation/user_manual/anuga_user_manual.tex (modified) (1 diff)
• source/anuga/abstract_2d_finite_volumes/neighbour_mesh.py (modified) (3 diffs)
• source/anuga/abstract_2d_finite_volumes/test_neighbour_mesh.py (modified) (3 diffs)
• source/anuga/abstract_2d_finite_volumes/util.py (modified) (1 diff)
• source/anuga/coordinate_transforms/geo_reference.py (modified) (1 diff)
• source/anuga/fit_interpolate/interpolate.py (modified) (4 diffs)
• source/anuga/shallow_water/data_manager.py (modified) (1 diff)
• source/anuga/shallow_water/test_data_manager.py (modified) (6 diffs)
• source/anuga/utilities/test_polygon.py (modified) (1 diff)

anuga_core/documentation/release_information/release_notes.txt

@@ -9 +9 @@
 More validation examples. 
 (up to jan 2008)
+
+
+time, Q = get_flow_through_cross_section(swwfile,
+                                         cross_section,
+                                         verbose=False)
+
+      
   
 

anuga_core/documentation/user_manual/anuga_user_manual.tex

@@ -2580 +2580 @@
   See doc string for general function sww2timeseries for details.
 
+\end{funcdesc}
+
+
+\begin{funcdesc}{get\_flow\_through\_cross\_section}{filename, cross\_section, verbose=False}
+  Module: \module{shallow\_water.data\_manager}
+
+  Obtain flow $[m^2]$ perpendicular to specified cross section.
+
+  Inputs:
+  \begin{itemize} 
+        \item filename: Name of sww file containing ANUGA model output.
+        \item polyline: Representation of desired cross section - it may contain multiple
+          sections allowing for complex shapes. Assume absolute UTM coordinates.
+  \end{itemize} 
+
+  Output:
+  \begin{itemize}
+    \item time: All stored times in sww file
+    \item Q: Hydrograph of total flow across given segments for all stored times.
+  \end{itemize} 
+  
+  The normal flow is computed for each triangle intersected by the polyline and
+  added up.  Multiple segments at different angles are specified the normal flows
+  may partially cancel each other.
+  
+  Example to find flow through cross section:
+  
+  \begin{verbatim}  
+  cross_section = [[x, 0], [x, width]]
+  time, Q = get_flow_through_cross_section(filename,
+                                           cross_section,
+                                           verbose=False)
+  \end{verbatim} 
+
+
+  See doc string for general function get_maximum_inundation_data for details.
+  
 \end{funcdesc}
 

anuga_core/source/anuga/abstract_2d_finite_volumes/neighbour_mesh.py

@@ -672 +672 @@
                 # Normalise
                 l_u = sqrt(u[0]*u[0] + u[1]*u[1])
-                l_v = sqrt(v[0]*v[0] + v[1]*v[1])                
-
-                x = (u[0]*v[0] + u[1]*v[1])/l_u/l_v # Inner product
+                l_v = sqrt(v[0]*v[0] + v[1]*v[1])
+
+                msg = 'Normal vector in triangle %d does not have unit length' %i
+                assert allclose(l_v, 1), msg
+
+                x = (u[0]*v[0] + u[1]*v[1])/l_u # Inner product
                 
                 msg = 'Normal vector (%f,%f) is not perpendicular to' %tuple(v)
@@ -948 +951 @@
 
               status, value = intersection(line, edge)
+              #if value is not None: print 'Triangle %d, Got' %i, status, value
+                  
               if status == 1:
                   # Normal intersection of one edge or vertex
@@ -960 +965 @@
                   # Edge is sharing a segment with line
 
-                  # This is currently covered by the two
+                  # This is usually covered by the two
                   # vertices that would have been picked up
-                  # under status == 1
-                  pass
+                  # under status == 1.
+                  # However, if coinciding line stops partway
+                  # along this edge, it will be recorded here.
+                  intersections[tuple(value[0,:])] = i
+                  intersections[tuple(value[1,:])] = i                                    
+
                   
-
-
           if len(intersections) == 1:
               # Check if either line end point lies fully within this triangle

anuga_core/source/anuga/abstract_2d_finite_volumes/test_neighbour_mesh.py

@@ -1338 +1338 @@
             
 
-        # Check all normals point straight down
+        # Check all 
         total_length = 0
-        for i, x in enumerate(L): 
+        for x in L:
             assert allclose(x.length, 1.0)
             assert allclose(x.normal, [0,-1])
@@ -1355 +1355 @@
 
         msg = 'Segments do not add up'
-        assert allclose(total_length, 2), msg            
+        assert allclose(total_length, 2), msg
+        
+
+    def test_get_intersecting_segments_partially_coinciding(self):
+        """test_get_intersecting_segments_partially_coinciding(self):
+
+        Test that line coinciding with triangle edges work.
+        But this ones only coincide with parts of the edge. 
+        """
+
+        # Build test mesh
+        
+        # Create basic mesh
+        # 9 points at (0,0), (0, 1), ..., (2,2)
+        # 8 triangles enumerated from left bottom to right top.
+        points, vertices, boundary = rectangular(2, 2, 2, 2)
+        mesh = Mesh(points, vertices, boundary)
+        intersected_triangles = [1,5]
+
+        # Horizontal line intersecting along center but stopping
+        # parway through second triangle's edge
+        #
+
+        y_line = 1.0
+        
+        #line = [[0, y_line], [2, y_line]]
+        line = [[0, y_line], [1.5, y_line]]
+
+        L = mesh.get_intersecting_segments(line)
+        #for x in L:
+        #    print x
+
+        msg = 'Two triangles should be returned'    
+        assert len(L) == 2, msg    
+            
+
+        # Check all 
+        total_length = 0
+        for x in L:
+            if x.triangle_id == 1:
+                assert allclose(x.length, 1)        
+                assert allclose(x.normal, [0, -1])
+                
+            if x.triangle_id == 5:
+                assert allclose(x.length, 0.5)
+                assert allclose(x.normal, [0, -1])
+
+
+            assert x.triangle_id in intersected_triangles
+
+            total_length += x.length
+
+        msg = 'Segments do not add up'
+        assert allclose(total_length, 1.5), msg            
 
 
@@ -1747 +1800 @@
 
 
+
+
+    def test_get_intersecting_segments7(self):
+        """test_get_intersecting_segments(self):
+
+        Check that line can stop inside a triangle - this is from
+        flow throug a cross sections example in test_datamanager.
+        
+        """
+
+        # Build test mesh
+        width = 5
+        length = 100
+        t_end = 1
+        points, vertices, boundary = rectangular(length, width,
+                                                 length, width)
+
+        mesh = Mesh(points, vertices, boundary)
+        
+
+        # A range of partial lines
+        x = length/2.
+        for i in range(10):
+            start_point = [length/2., i*width/10.]
+            #print 
+            #print start_point
+                            
+            line = [start_point, [length/2., width]]
+ 
+            L = mesh.get_intersecting_segments(line)
+
+            if start_point[1] < 1:
+                assert len(L) == 5
+                
+            
+            total_length = 0    
+            for x in L:
+                total_length += x.length
+                
+
+            ref_length = line[1][1] - line[0][1]
+            #print ref_length, total_length
+            assert allclose(total_length, ref_length)
+
+
+
         
 #-------------------------------------------------------------
 if __name__ == "__main__":
     #suite = unittest.makeSuite(Test_Mesh,'test_two_triangles')
-    #suite = unittest.makeSuite(Test_Mesh,'test_get_intersecting_segments_coinciding')
-    suite = unittest.makeSuite(Test_Mesh,'test') #_get_intersecting_segments6')
+    #suite = unittest.makeSuite(Test_Mesh,'test_get_intersecting_segments_partially_coinciding')
+    #suite = unittest.makeSuite(Test_Mesh,'test_get_intersecting_segments7')
+    suite = unittest.makeSuite(Test_Mesh,'test')
     runner = unittest.TextTestRunner()
     runner.run(suite)

anuga_core/source/anuga/abstract_2d_finite_volumes/util.py

@@ -305 +305 @@
         if interpolation_points is not None:
             # Adjust for georef
+
+            # FIXME (Ole): Use geo_reference.get_relative
             interpolation_points[:,0] -= xllcorner
             interpolation_points[:,1] -= yllcorner        

anuga_core/source/anuga/coordinate_transforms/geo_reference.py

@@ -274 +274 @@
 
 
+    def get_relative(self, points):
+        """
+        Given a set of points in absolute UTM coordinates,
+        make them relative to this geo_reference instance,
+        return the points as relative values.
+
+        This is the inverse of get_absolute.
+        """
+
+        is_list = False
+        if type(points) == types.ListType:
+            is_list = True
+
+        points = ensure_numeric(points, Float)
+        if len(points.shape) == 1:
+            #One point has been passed
+            msg = 'Single point must have two elements'
+            if not len(points) == 2:
+                raise ShapeError, msg    
+                #points = reshape(points, (1,2))
+
+
+        msg = 'Input must be an N x 2 array or list of (x,y) values. '
+        msg += 'I got an %d x %d array' %points.shape    
+        if not points.shape[1] == 2:
+            raise ShapeError, msg    
+            
+        
+        # Subtract geo ref from points
+        if not self.is_absolute():
+            points[:,0] -= self.xllcorner 
+            points[:,1] -= self.yllcorner
+
+        
+        if is_list:
+            points = points.tolist()
+             
+        return points
+
+
+
     def reconcile_zones(self, other):
         if other is None:

anuga_core/source/anuga/fit_interpolate/interpolate.py

@@ -439 +439 @@
 
     Let m be the number of vertices, n the number of triangles
-    and p the number of timesteps. 
+    and p the number of timesteps.
+    Also, let N be the number of interpolation points.
 
     Mandatory input
@@ -446 +447 @@
                             The arrays must either have dimensions pxm or mx1.
                             The resulting function will be time dependent in
-                            the former case while it will be constan with
+                            the former case while it will be constant with
                             respect to time in the latter case.
         
     Optional input:
-        quantity_names:     List of keys into the quantities dictionary 
+        quantity_names:     List of keys into the quantities dictionary for
+                            imposing a particular order on the output vector.
         vertex_coordinates: mx2 array of coordinates (Float)
         triangles:          nx3 array of indices into vertex_coordinates (Int)
@@ -470 +472 @@
     geospatial data objects
 
-    Time assumed to be relative to starttime    
+    Time assumed to be relative to starttime (FIXME (Ole): This comment should be removed)
+    All coordinates assume origin of (0,0) - e.g. georeferencing must be taken care of
+    outside this function
     """
   
@@ -522 +526 @@
         if vertex_coordinates is None:
             self.spatial = False
-        else:    
+        else:
+            # FIXME (Ole): Try ensure_numeric here -
+            #this function knows nothing about georefering.
             vertex_coordinates = ensure_absolute(vertex_coordinates)
 

anuga_core/source/anuga/shallow_water/data_manager.py

@@ -5621 +5621 @@
         filename: Name of sww file
         polyline: Representation of desired cross section - it may contain multiple
-                  sections allowing for complex shapes.
+                  sections allowing for complex shapes. Assume absolute UTM coordinates.
+                  Format [[x0, y0], [x1, y1], ...]
 
     Output:
-        time: All stored times
+        time: All stored times in sww file
         Q: Hydrograph of total flow across given segments for all stored times.
 
     The normal flow is computed for each triangle intersected by the polyline and
-    added up.  multiple segments at different angles are specified the normal flows
+    added up.  Multiple segments at different angles are specified the normal flows
     may partially cancel each other.
 
-    """
+    The typical usage of this function would be to get flow through a channel,
+    and the polyline would then be a cross section perpendicular to the flow.
+
+    """
+
+    from anuga.fit_interpolate.interpolate import Interpolation_function
+
+    quantity_names =['elevation',
+                     'stage',
+                     'xmomentum',
+                     'ymomentum']
 
     # Get mesh and quantities from sww file
     X = get_mesh_and_quantities_from_file(filename,
-                                          quantities=['elevation',
-                                                      'stage',
-                                                      'xmomentum',
-                                                      'ymomentum'], 
+                                          quantities=quantity_names,
                                           verbose=verbose)
     mesh, quantities, time = X
 
 
+    # Adjust polyline to mesh spatial origin
+    polyline = mesh.geo_reference.get_relative(polyline)
     
     # Find all intersections and associated triangles.
-    
-    mesh.get_intersecting_segments(polyline)
+    segments = mesh.get_intersecting_segments(polyline)
+    #print
+    #for x in segments:
+    #    print x
+
     
     # Then store for each triangle the length of the intersecting segment(s),
-    # right hand normal(s) and midpoints.
+    # right hand normal(s) and midpoints as interpolation_points.
+    interpolation_points = []
+    for segment in segments:
+        midpoint = sum(array(segment.segment))/2
+        interpolation_points.append(midpoint)
+
+
+    I = Interpolation_function(time,
+                               quantities,
+                               quantity_names=quantity_names,
+                               vertex_coordinates=mesh.nodes,
+                               triangles=mesh.triangles,
+                               interpolation_points=interpolation_points,
+                               verbose=verbose)
+
+    # Compute hydrograph
+    Q = []
+    for t in time:
+        total_flow=0
+        for i, p in enumerate(interpolation_points):
+            elevation, stage, uh, vh = I(t, point_id=i)
+            normal = segments[i].normal
+
+            # Inner product of momentum vector with segment normal [m^2/s] 
+            normal_momentum = uh*normal[0] + vh*normal[1] 
+
+            # Flow across this segment [m^3/s]
+            segment_flow = normal_momentum*segments[i].length
+
+            # Accumulate
+            total_flow += segment_flow
+             
+
+        # Store flow at this timestep    
+        Q.append(total_flow)
+    
+
 
     return time, Q

anuga_core/source/anuga/shallow_water/test_data_manager.py

@@ -7493 +7493 @@
         
 
-    def NOtest_get_flow_through_cross_section(self):
+    def test_get_flow_through_cross_section(self):
         """test_get_flow_through_cross_section(self):
 
@@ -7507 +7507 @@
         w = 5 m (width of channel)
 
-        q = u*h*w = 10 m^3/s 
+        q = u*h*w = 10 m^3/s
+
+        #---------- First run without geo referencing        
         
         """
@@ -7521 +7523 @@
         width = 5
         length = 100
-        t_end = 10
+        t_end = 3
         points, vertices, boundary = rectangular(length, width,
                                                  length, width)
@@ -7545 +7547 @@
 
 
-        #---------- First run without geo referencing
+
         
         domain.set_quantity('elevation', 0.0)
@@ -7571 +7573 @@
             
 
-
-        # Check flow through the middle
-        cross_section = [[length/2., 0], [length/2., width]]
-        Q = get_flow_through_cross_section(swwfile,
-                                           cross_section,
-                                           verbose=True)
-
-        assert allclose(Q, uh*width)  
+        # Check flows through the middle
+        for i in range(10):
+            x = length/2. + i*0.23674563 # Arbitrary
+            cross_section = [[x, 0], [x, width]]
+            time, Q = get_flow_through_cross_section(swwfile,
+                                                     cross_section,
+                                                     verbose=False)
+
+            assert allclose(Q, uh*width)
+
+
+       
+        # Try the same with partial lines
+        x = length/2.
+        for i in range(10):
+            start_point = [length/2., i*width/10.]
+            #print start_point
+                            
+            cross_section = [start_point, [length/2., width]]
+            time, Q = get_flow_through_cross_section(swwfile,
+                                                     cross_section,
+                                                     verbose=False)
+
+            #print i, Q, (width-start_point[1])
+            assert allclose(Q, uh*(width-start_point[1]))
+
+
+        # Verify no flow when line is parallel to flow
+        cross_section = [[length/2.-10, width/2.], [length/2.+10, width/2.]]
+        time, Q = get_flow_through_cross_section(swwfile,
+                                                 cross_section,
+                                                 verbose=False)
+
+        #print i, Q
+        assert allclose(Q, 0)        
+
+
+        # Try with lines on an angle (all flow still runs through here)
+        cross_section = [[length/2., 0], [length/2.+width, width]]
+        time, Q = get_flow_through_cross_section(swwfile,
+                                                 cross_section,
+                                                 verbose=False)
+
+        assert allclose(Q, uh*width)        
+        
+
+
                                       
-
-
-
-
-        
+    def test_get_flow_through_cross_section_with_geo(self):
+        """test_get_flow_through_cross_section(self):
+
+        Test that the total flow 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 flow.
+
+        The specifics are
+        u = 2 m/s
+        h = 1 m
+        w = 5 m (width of channel)
+
+        q = u*h*w = 10 m^3/s
+
+
+        This run tries it with 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 (100m x 5m)
+        width = 5
+        length = 100
+        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
+
+        h = 1.0
+        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)
+        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):
+                assert allclose(f(t, i), [1, 2, 0])
+            
+
+        # Check flows through the middle
+        for i in range(10):
+            x = length/2. + i*0.23674563 # Arbitrary
+            cross_section = [[x, 0], [x, width]]
+
+            cross_section = domain.geo_reference.get_absolute(cross_section)            
+            time, Q = get_flow_through_cross_section(swwfile,
+                                                     cross_section,
+                                                     verbose=False)
+
+            assert allclose(Q, uh*width)
+
+
         
     def test_get_all_swwfiles(self):
@@ -7694 +7827 @@
 #    suite = unittest.makeSuite(Test_Data_Manager,'test_screen_catcher')
     suite = unittest.makeSuite(Test_Data_Manager,'test')
-    #suite = unittest.makeSuite(Test_Data_Manager,'test_get_flow_through_cross_section')
+    #suite = unittest.makeSuite(Test_Data_Manager,'test_get_flow_through_cross_section_with_geo')
     #suite = unittest.makeSuite(Test_Data_Manager,'covered_')
 

anuga_core/source/anuga/utilities/test_polygon.py

@@ -627 +627 @@
         assert allclose(value, [14.068965517, 7.0344827586])        
 
+
+    def test_intersection_endpoints(self):
+        """test_intersection_endpoints(self):
+
+        Test that coinciding endpoints are picked up
+        """
+        line0 = [[0,0], [1,1]]
+        line1 = [[1,1], [2,1]]
+
+        status, value = intersection(line0, line1)
+        assert status == 1
+        assert allclose(value, [1.0, 1.0])
+
+
+        line0 = [[1,1], [2,0]]
+        line1 = [[1,1], [2,1]]
+
+        status, value = intersection(line0, line1)
+        assert status == 1
+        assert allclose(value, [1.0, 1.0])        
+        
 
     def test_intersection_direction_invariance(self):