Changeset 5288


Ignore:
Timestamp:
May 7, 2008, 4:01:52 PM (16 years ago)
Author:
ole
Message:

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

Location:
anuga_core
Files:
10 edited

Legend:

Unmodified
Added
Removed
  • anuga_core/documentation/release_information/release_notes.txt

    r5237 r5288  
    99More validation examples.
    1010(up to jan 2008)
     11
     12
     13time, Q = get_flow_through_cross_section(swwfile,
     14                                         cross_section,
     15                                         verbose=False)
     16
     17     
    1118 
    1219
  • anuga_core/documentation/user_manual/anuga_user_manual.tex

    r5246 r5288  
    25802580  See doc string for general function sww2timeseries for details.
    25812581
     2582\end{funcdesc}
     2583
     2584
     2585\begin{funcdesc}{get\_flow\_through\_cross\_section}{filename, cross\_section, verbose=False}
     2586  Module: \module{shallow\_water.data\_manager}
     2587
     2588  Obtain flow $[m^2]$ perpendicular to specified cross section.
     2589
     2590  Inputs:
     2591  \begin{itemize}
     2592        \item filename: Name of sww file containing ANUGA model output.
     2593        \item polyline: Representation of desired cross section - it may contain multiple
     2594          sections allowing for complex shapes. Assume absolute UTM coordinates.
     2595  \end{itemize}
     2596
     2597  Output:
     2598  \begin{itemize}
     2599    \item time: All stored times in sww file
     2600    \item Q: Hydrograph of total flow across given segments for all stored times.
     2601  \end{itemize}
     2602 
     2603  The normal flow is computed for each triangle intersected by the polyline and
     2604  added up.  Multiple segments at different angles are specified the normal flows
     2605  may partially cancel each other.
     2606 
     2607  Example to find flow through cross section:
     2608 
     2609  \begin{verbatim} 
     2610  cross_section = [[x, 0], [x, width]]
     2611  time, Q = get_flow_through_cross_section(filename,
     2612                                           cross_section,
     2613                                           verbose=False)
     2614  \end{verbatim}
     2615
     2616
     2617  See doc string for general function get_maximum_inundation_data for details.
     2618 
    25822619\end{funcdesc}
    25832620
  • anuga_core/source/anuga/abstract_2d_finite_volumes/neighbour_mesh.py

    r5287 r5288  
    672672                # Normalise
    673673                l_u = sqrt(u[0]*u[0] + u[1]*u[1])
    674                 l_v = sqrt(v[0]*v[0] + v[1]*v[1])               
    675 
    676                 x = (u[0]*v[0] + u[1]*v[1])/l_u/l_v # Inner product
     674                l_v = sqrt(v[0]*v[0] + v[1]*v[1])
     675
     676                msg = 'Normal vector in triangle %d does not have unit length' %i
     677                assert allclose(l_v, 1), msg
     678
     679                x = (u[0]*v[0] + u[1]*v[1])/l_u # Inner product
    677680               
    678681                msg = 'Normal vector (%f,%f) is not perpendicular to' %tuple(v)
     
    948951
    949952              status, value = intersection(line, edge)
     953              #if value is not None: print 'Triangle %d, Got' %i, status, value
     954                 
    950955              if status == 1:
    951956                  # Normal intersection of one edge or vertex
     
    960965                  # Edge is sharing a segment with line
    961966
    962                   # This is currently covered by the two
     967                  # This is usually covered by the two
    963968                  # vertices that would have been picked up
    964                   # under status == 1
    965                   pass
     969                  # under status == 1.
     970                  # However, if coinciding line stops partway
     971                  # along this edge, it will be recorded here.
     972                  intersections[tuple(value[0,:])] = i
     973                  intersections[tuple(value[1,:])] = i                                   
     974
    966975                 
    967 
    968 
    969976          if len(intersections) == 1:
    970977              # Check if either line end point lies fully within this triangle
  • anuga_core/source/anuga/abstract_2d_finite_volumes/test_neighbour_mesh.py

    r5287 r5288  
    13381338           
    13391339
    1340         # Check all normals point straight down
     1340        # Check all
    13411341        total_length = 0
    1342         for i, x in enumerate(L):
     1342        for x in L:
    13431343            assert allclose(x.length, 1.0)
    13441344            assert allclose(x.normal, [0,-1])
     
    13551355
    13561356        msg = 'Segments do not add up'
    1357         assert allclose(total_length, 2), msg           
     1357        assert allclose(total_length, 2), msg
     1358       
     1359
     1360    def test_get_intersecting_segments_partially_coinciding(self):
     1361        """test_get_intersecting_segments_partially_coinciding(self):
     1362
     1363        Test that line coinciding with triangle edges work.
     1364        But this ones only coincide with parts of the edge.
     1365        """
     1366
     1367        # Build test mesh
     1368       
     1369        # Create basic mesh
     1370        # 9 points at (0,0), (0, 1), ..., (2,2)
     1371        # 8 triangles enumerated from left bottom to right top.
     1372        points, vertices, boundary = rectangular(2, 2, 2, 2)
     1373        mesh = Mesh(points, vertices, boundary)
     1374        intersected_triangles = [1,5]
     1375
     1376        # Horizontal line intersecting along center but stopping
     1377        # parway through second triangle's edge
     1378        #
     1379
     1380        y_line = 1.0
     1381       
     1382        #line = [[0, y_line], [2, y_line]]
     1383        line = [[0, y_line], [1.5, y_line]]
     1384
     1385        L = mesh.get_intersecting_segments(line)
     1386        #for x in L:
     1387        #    print x
     1388
     1389        msg = 'Two triangles should be returned'   
     1390        assert len(L) == 2, msg   
     1391           
     1392
     1393        # Check all
     1394        total_length = 0
     1395        for x in L:
     1396            if x.triangle_id == 1:
     1397                assert allclose(x.length, 1)       
     1398                assert allclose(x.normal, [0, -1])
     1399               
     1400            if x.triangle_id == 5:
     1401                assert allclose(x.length, 0.5)
     1402                assert allclose(x.normal, [0, -1])
     1403
     1404
     1405            assert x.triangle_id in intersected_triangles
     1406
     1407            total_length += x.length
     1408
     1409        msg = 'Segments do not add up'
     1410        assert allclose(total_length, 1.5), msg           
    13581411
    13591412
     
    17471800
    17481801
     1802
     1803
     1804    def test_get_intersecting_segments7(self):
     1805        """test_get_intersecting_segments(self):
     1806
     1807        Check that line can stop inside a triangle - this is from
     1808        flow throug a cross sections example in test_datamanager.
     1809       
     1810        """
     1811
     1812        # Build test mesh
     1813        width = 5
     1814        length = 100
     1815        t_end = 1
     1816        points, vertices, boundary = rectangular(length, width,
     1817                                                 length, width)
     1818
     1819        mesh = Mesh(points, vertices, boundary)
     1820       
     1821
     1822        # A range of partial lines
     1823        x = length/2.
     1824        for i in range(10):
     1825            start_point = [length/2., i*width/10.]
     1826            #print
     1827            #print start_point
     1828                           
     1829            line = [start_point, [length/2., width]]
     1830 
     1831            L = mesh.get_intersecting_segments(line)
     1832
     1833            if start_point[1] < 1:
     1834                assert len(L) == 5
     1835               
     1836           
     1837            total_length = 0   
     1838            for x in L:
     1839                total_length += x.length
     1840               
     1841
     1842            ref_length = line[1][1] - line[0][1]
     1843            #print ref_length, total_length
     1844            assert allclose(total_length, ref_length)
     1845
     1846
     1847
    17491848       
    17501849#-------------------------------------------------------------
    17511850if __name__ == "__main__":
    17521851    #suite = unittest.makeSuite(Test_Mesh,'test_two_triangles')
    1753     #suite = unittest.makeSuite(Test_Mesh,'test_get_intersecting_segments_coinciding')
    1754     suite = unittest.makeSuite(Test_Mesh,'test') #_get_intersecting_segments6')
     1852    #suite = unittest.makeSuite(Test_Mesh,'test_get_intersecting_segments_partially_coinciding')
     1853    #suite = unittest.makeSuite(Test_Mesh,'test_get_intersecting_segments7')
     1854    suite = unittest.makeSuite(Test_Mesh,'test')
    17551855    runner = unittest.TextTestRunner()
    17561856    runner.run(suite)
  • anuga_core/source/anuga/abstract_2d_finite_volumes/util.py

    r5221 r5288  
    305305        if interpolation_points is not None:
    306306            # Adjust for georef
     307
     308            # FIXME (Ole): Use geo_reference.get_relative
    307309            interpolation_points[:,0] -= xllcorner
    308310            interpolation_points[:,1] -= yllcorner       
  • anuga_core/source/anuga/coordinate_transforms/geo_reference.py

    r4839 r5288  
    274274
    275275
     276    def get_relative(self, points):
     277        """
     278        Given a set of points in absolute UTM coordinates,
     279        make them relative to this geo_reference instance,
     280        return the points as relative values.
     281
     282        This is the inverse of get_absolute.
     283        """
     284
     285        is_list = False
     286        if type(points) == types.ListType:
     287            is_list = True
     288
     289        points = ensure_numeric(points, Float)
     290        if len(points.shape) == 1:
     291            #One point has been passed
     292            msg = 'Single point must have two elements'
     293            if not len(points) == 2:
     294                raise ShapeError, msg   
     295                #points = reshape(points, (1,2))
     296
     297
     298        msg = 'Input must be an N x 2 array or list of (x,y) values. '
     299        msg += 'I got an %d x %d array' %points.shape   
     300        if not points.shape[1] == 2:
     301            raise ShapeError, msg   
     302           
     303       
     304        # Subtract geo ref from points
     305        if not self.is_absolute():
     306            points[:,0] -= self.xllcorner
     307            points[:,1] -= self.yllcorner
     308
     309       
     310        if is_list:
     311            points = points.tolist()
     312             
     313        return points
     314
     315
     316
    276317    def reconcile_zones(self, other):
    277318        if other is None:
  • anuga_core/source/anuga/fit_interpolate/interpolate.py

    r5222 r5288  
    439439
    440440    Let m be the number of vertices, n the number of triangles
    441     and p the number of timesteps.
     441    and p the number of timesteps.
     442    Also, let N be the number of interpolation points.
    442443
    443444    Mandatory input
     
    446447                            The arrays must either have dimensions pxm or mx1.
    447448                            The resulting function will be time dependent in
    448                             the former case while it will be constan with
     449                            the former case while it will be constant with
    449450                            respect to time in the latter case.
    450451       
    451452    Optional input:
    452         quantity_names:     List of keys into the quantities dictionary
     453        quantity_names:     List of keys into the quantities dictionary for
     454                            imposing a particular order on the output vector.
    453455        vertex_coordinates: mx2 array of coordinates (Float)
    454456        triangles:          nx3 array of indices into vertex_coordinates (Int)
     
    470472    geospatial data objects
    471473
    472     Time assumed to be relative to starttime   
     474    Time assumed to be relative to starttime (FIXME (Ole): This comment should be removed)
     475    All coordinates assume origin of (0,0) - e.g. georeferencing must be taken care of
     476    outside this function
    473477    """
    474478 
     
    522526        if vertex_coordinates is None:
    523527            self.spatial = False
    524         else:   
     528        else:
     529            # FIXME (Ole): Try ensure_numeric here -
     530            #this function knows nothing about georefering.
    525531            vertex_coordinates = ensure_absolute(vertex_coordinates)
    526532
  • anuga_core/source/anuga/shallow_water/data_manager.py

    r5276 r5288  
    56215621        filename: Name of sww file
    56225622        polyline: Representation of desired cross section - it may contain multiple
    5623                   sections allowing for complex shapes.
     5623                  sections allowing for complex shapes. Assume absolute UTM coordinates.
     5624                  Format [[x0, y0], [x1, y1], ...]
    56245625
    56255626    Output:
    5626         time: All stored times
     5627        time: All stored times in sww file
    56275628        Q: Hydrograph of total flow across given segments for all stored times.
    56285629
    56295630    The normal flow is computed for each triangle intersected by the polyline and
    5630     added up.  multiple segments at different angles are specified the normal flows
     5631    added up.  Multiple segments at different angles are specified the normal flows
    56315632    may partially cancel each other.
    56325633
    5633     """
     5634    The typical usage of this function would be to get flow through a channel,
     5635    and the polyline would then be a cross section perpendicular to the flow.
     5636
     5637    """
     5638
     5639    from anuga.fit_interpolate.interpolate import Interpolation_function
     5640
     5641    quantity_names =['elevation',
     5642                     'stage',
     5643                     'xmomentum',
     5644                     'ymomentum']
    56345645
    56355646    # Get mesh and quantities from sww file
    56365647    X = get_mesh_and_quantities_from_file(filename,
    5637                                           quantities=['elevation',
    5638                                                       'stage',
    5639                                                       'xmomentum',
    5640                                                       'ymomentum'],
     5648                                          quantities=quantity_names,
    56415649                                          verbose=verbose)
    56425650    mesh, quantities, time = X
    56435651
    56445652
     5653    # Adjust polyline to mesh spatial origin
     5654    polyline = mesh.geo_reference.get_relative(polyline)
    56455655   
    56465656    # Find all intersections and associated triangles.
    5647    
    5648     mesh.get_intersecting_segments(polyline)
     5657    segments = mesh.get_intersecting_segments(polyline)
     5658    #print
     5659    #for x in segments:
     5660    #    print x
     5661
    56495662   
    56505663    # Then store for each triangle the length of the intersecting segment(s),
    5651     # right hand normal(s) and midpoints.
     5664    # right hand normal(s) and midpoints as interpolation_points.
     5665    interpolation_points = []
     5666    for segment in segments:
     5667        midpoint = sum(array(segment.segment))/2
     5668        interpolation_points.append(midpoint)
     5669
     5670
     5671    I = Interpolation_function(time,
     5672                               quantities,
     5673                               quantity_names=quantity_names,
     5674                               vertex_coordinates=mesh.nodes,
     5675                               triangles=mesh.triangles,
     5676                               interpolation_points=interpolation_points,
     5677                               verbose=verbose)
     5678
     5679    # Compute hydrograph
     5680    Q = []
     5681    for t in time:
     5682        total_flow=0
     5683        for i, p in enumerate(interpolation_points):
     5684            elevation, stage, uh, vh = I(t, point_id=i)
     5685            normal = segments[i].normal
     5686
     5687            # Inner product of momentum vector with segment normal [m^2/s]
     5688            normal_momentum = uh*normal[0] + vh*normal[1]
     5689
     5690            # Flow across this segment [m^3/s]
     5691            segment_flow = normal_momentum*segments[i].length
     5692
     5693            # Accumulate
     5694            total_flow += segment_flow
     5695             
     5696
     5697        # Store flow at this timestep   
     5698        Q.append(total_flow)
     5699   
     5700
    56525701
    56535702    return time, Q
  • anuga_core/source/anuga/shallow_water/test_data_manager.py

    r5276 r5288  
    74937493       
    74947494
    7495     def NOtest_get_flow_through_cross_section(self):
     7495    def test_get_flow_through_cross_section(self):
    74967496        """test_get_flow_through_cross_section(self):
    74977497
     
    75077507        w = 5 m (width of channel)
    75087508
    7509         q = u*h*w = 10 m^3/s
     7509        q = u*h*w = 10 m^3/s
     7510
     7511        #---------- First run without geo referencing       
    75107512       
    75117513        """
     
    75217523        width = 5
    75227524        length = 100
    7523         t_end = 10
     7525        t_end = 3
    75247526        points, vertices, boundary = rectangular(length, width,
    75257527                                                 length, width)
     
    75457547
    75467548
    7547         #---------- First run without geo referencing
     7549
    75487550       
    75497551        domain.set_quantity('elevation', 0.0)
     
    75717573           
    75727574
    7573 
    7574         # Check flow through the middle
    7575         cross_section = [[length/2., 0], [length/2., width]]
    7576         Q = get_flow_through_cross_section(swwfile,
    7577                                            cross_section,
    7578                                            verbose=True)
    7579 
    7580         assert allclose(Q, uh*width) 
     7575        # Check flows through the middle
     7576        for i in range(10):
     7577            x = length/2. + i*0.23674563 # Arbitrary
     7578            cross_section = [[x, 0], [x, width]]
     7579            time, Q = get_flow_through_cross_section(swwfile,
     7580                                                     cross_section,
     7581                                                     verbose=False)
     7582
     7583            assert allclose(Q, uh*width)
     7584
     7585
     7586       
     7587        # Try the same with partial lines
     7588        x = length/2.
     7589        for i in range(10):
     7590            start_point = [length/2., i*width/10.]
     7591            #print start_point
     7592                           
     7593            cross_section = [start_point, [length/2., width]]
     7594            time, Q = get_flow_through_cross_section(swwfile,
     7595                                                     cross_section,
     7596                                                     verbose=False)
     7597
     7598            #print i, Q, (width-start_point[1])
     7599            assert allclose(Q, uh*(width-start_point[1]))
     7600
     7601
     7602        # Verify no flow when line is parallel to flow
     7603        cross_section = [[length/2.-10, width/2.], [length/2.+10, width/2.]]
     7604        time, Q = get_flow_through_cross_section(swwfile,
     7605                                                 cross_section,
     7606                                                 verbose=False)
     7607
     7608        #print i, Q
     7609        assert allclose(Q, 0)       
     7610
     7611
     7612        # Try with lines on an angle (all flow still runs through here)
     7613        cross_section = [[length/2., 0], [length/2.+width, width]]
     7614        time, Q = get_flow_through_cross_section(swwfile,
     7615                                                 cross_section,
     7616                                                 verbose=False)
     7617
     7618        assert allclose(Q, uh*width)       
     7619       
     7620
     7621
    75817622                                     
    7582 
    7583 
    7584 
    7585 
    7586        
     7623    def test_get_flow_through_cross_section_with_geo(self):
     7624        """test_get_flow_through_cross_section(self):
     7625
     7626        Test that the total flow through a cross section can be
     7627        correctly obtained from an sww file.
     7628       
     7629        This test creates a flat bed with a known flow through it and tests
     7630        that the function correctly returns the expected flow.
     7631
     7632        The specifics are
     7633        u = 2 m/s
     7634        h = 1 m
     7635        w = 5 m (width of channel)
     7636
     7637        q = u*h*w = 10 m^3/s
     7638
     7639
     7640        This run tries it with georeferencing
     7641       
     7642        """
     7643
     7644        import time, os
     7645        from Numeric import array, zeros, allclose, Float, concatenate
     7646        from Scientific.IO.NetCDF import NetCDFFile
     7647
     7648        # Setup
     7649        from mesh_factory import rectangular
     7650
     7651        # Create basic mesh (100m x 5m)
     7652        width = 5
     7653        length = 100
     7654        t_end = 1
     7655        points, vertices, boundary = rectangular(length, width,
     7656                                                 length, width)
     7657
     7658        # Create shallow water domain
     7659        domain = Domain(points, vertices, boundary,
     7660                        geo_reference = Geo_reference(56,308500,6189000))
     7661
     7662        domain.default_order = 2
     7663        domain.set_minimum_storable_height(0.01)
     7664
     7665        domain.set_name('flowtest')
     7666        swwfile = domain.get_name() + '.sww'
     7667
     7668        domain.set_datadir('.')
     7669        domain.format = 'sww'
     7670        domain.smooth = True
     7671
     7672        h = 1.0
     7673        u = 2.0
     7674        uh = u*h
     7675
     7676        Br = Reflective_boundary(domain)     # Side walls
     7677        Bd = Dirichlet_boundary([h, uh, 0])  # 2 m/s across the 5 m inlet:
     7678
     7679
     7680
     7681       
     7682        domain.set_quantity('elevation', 0.0)
     7683        domain.set_quantity('stage', h)
     7684        domain.set_quantity('xmomentum', uh)
     7685        domain.set_boundary( {'left': Bd, 'right': Bd, 'top': Br, 'bottom': Br})
     7686
     7687        for t in domain.evolve(yieldstep=1, finaltime = t_end):
     7688            pass
     7689
     7690        # Check that momentum is as it should be in the interior
     7691
     7692        I = [[0, width/2.],
     7693             [length/2., width/2.],
     7694             [length, width/2.]]
     7695       
     7696        I = domain.geo_reference.get_absolute(I)
     7697        f = file_function(swwfile,
     7698                          quantities=['stage', 'xmomentum', 'ymomentum'],
     7699                          interpolation_points=I,
     7700                          verbose=False)
     7701
     7702        for t in range(t_end+1):
     7703            for i in range(3):
     7704                assert allclose(f(t, i), [1, 2, 0])
     7705           
     7706
     7707        # Check flows through the middle
     7708        for i in range(10):
     7709            x = length/2. + i*0.23674563 # Arbitrary
     7710            cross_section = [[x, 0], [x, width]]
     7711
     7712            cross_section = domain.geo_reference.get_absolute(cross_section)           
     7713            time, Q = get_flow_through_cross_section(swwfile,
     7714                                                     cross_section,
     7715                                                     verbose=False)
     7716
     7717            assert allclose(Q, uh*width)
     7718
     7719
    75877720       
    75887721    def test_get_all_swwfiles(self):
     
    76947827#    suite = unittest.makeSuite(Test_Data_Manager,'test_screen_catcher')
    76957828    suite = unittest.makeSuite(Test_Data_Manager,'test')
    7696     #suite = unittest.makeSuite(Test_Data_Manager,'test_get_flow_through_cross_section')
     7829    #suite = unittest.makeSuite(Test_Data_Manager,'test_get_flow_through_cross_section_with_geo')
    76977830    #suite = unittest.makeSuite(Test_Data_Manager,'covered_')
    76987831
  • anuga_core/source/anuga/utilities/test_polygon.py

    r5282 r5288  
    627627        assert allclose(value, [14.068965517, 7.0344827586])       
    628628
     629
     630    def test_intersection_endpoints(self):
     631        """test_intersection_endpoints(self):
     632
     633        Test that coinciding endpoints are picked up
     634        """
     635        line0 = [[0,0], [1,1]]
     636        line1 = [[1,1], [2,1]]
     637
     638        status, value = intersection(line0, line1)
     639        assert status == 1
     640        assert allclose(value, [1.0, 1.0])
     641
     642
     643        line0 = [[1,1], [2,0]]
     644        line1 = [[1,1], [2,1]]
     645
     646        status, value = intersection(line0, line1)
     647        assert status == 1
     648        assert allclose(value, [1.0, 1.0])       
     649       
    629650
    630651    def test_intersection_direction_invariance(self):
Note: See TracChangeset for help on using the changeset viewer.