Changeset 5729


Ignore:
Timestamp:
Sep 4, 2008, 2:45:27 PM (14 years ago)
Author:
ole
Message:

Refactored get_flow_through_cross_section and added a runtime version to
shallow_water_domain. It still needs some work and testing.

Location:
anuga_core/source/anuga
Files:
11 edited

Legend:

Unmodified
Added
Removed
  • anuga_core/source/anuga/abstract_2d_finite_volumes/domain.py

    r5666 r5729  
    369369
    370370    def get_quantity(self, name, location='vertices', indices = None):
    371         """Get quantity object.
     371        """Get pointer to quantity object.
    372372
    373373        name: Name of quantity
  • anuga_core/source/anuga/abstract_2d_finite_volumes/neighbour_mesh.py

    r5321 r5729  
    11131113       
    11141114
     1115       
    11151116
    11161117class Triangle_intersection:
     
    11471148   
    11481149        return s
     1150
     1151def segment_midpoints(segments):
     1152    """Calculate midpoints of all segments
     1153   
     1154    Inputs:
     1155       segments: List of instances of class Segment
     1156       
     1157    Ouputs:
     1158       midpoints: List of points
     1159    """
     1160   
     1161    midpoints = []
     1162    msg = 'Elements of input list to segment_midpoints must be of class Triangle_intersection'
     1163    for segment in segments:
     1164        assert isinstance(segment, Triangle_intersection), msg
     1165       
     1166        midpoint = sum(array(segment.segment))/2
     1167        midpoints.append(midpoint)
     1168
     1169    return midpoints
     1170   
     1171   
     1172   
  • anuga_core/source/anuga/abstract_2d_finite_volumes/quantity.py

    r5522 r5729  
    753753                                        verbose=False,
    754754                                        use_cache=False):
    755         # FIXME: Use this function for the time being. Later move code in here
     755        """ Set values based on geo referenced geospatial data object.
     756        """
    756757
    757758        points = geospatial_data.get_data_points(absolute=False)
    758759        values = geospatial_data.get_attributes()
    759760        data_georef = geospatial_data.get_geo_reference()
    760 
    761761
    762762
     
    10071007
    10081008    def get_interpolated_values(self, interpolation_points):
    1009 
    1010         # Interpolation object based on internal (discontinuous triangles)
     1009        """ Get values at interpolation points
     1010       
     1011        If interpolation points have been given previously, the
     1012        associated matrices will be reused to save time.
     1013       
     1014        The argument interpolation points must be given as either a list of absolute UTM coordinates or
     1015        a geospatial data object.
     1016        """
     1017
     1018        # FIXME (Ole): Could do with an input check (should be generalised and used widely)
     1019        # That will check that interpolation points is either a list of points, Nx2 array, or geospatial
     1020       
     1021        # Ensure points are converted to coordinates relative to mesh origin
     1022        # FIXME (Ole): This could all be refactored using the 'change_points_geo_ref' method
     1023        # of Class geo_reference. The purpose is to make interpolation points relative
     1024        # to the mesh origin.
     1025        #
     1026        # Speed is also a consideration here.
     1027       
     1028        if isinstance(interpolation_points, Geospatial_data):       
     1029            # Ensure interpolation points are in absolute UTM coordinates
     1030            interpolation_points = interpolation_points.get_data_points(absolute=True)
     1031           
     1032        # Reconcile interpolation points with georeference of domain
     1033        interpolation_points = self.domain.geo_reference.get_relative(interpolation_points)
     1034        interpolation_points = ensure_numeric(interpolation_points)
     1035       
     1036        # Get internal (discontinuous) triangles for use with the interpolation object.
    10111037        x, y, vertex_values, triangles = self.get_vertex_values(xy=True,
    10121038                                                                smooth=False)
    1013         # FIXME: This concat should roll into get_vertex_values
     1039        # FIXME (Ole): This concat should roll into get_vertex_values
    10141040        vertex_coordinates = concatenate((x[:, NewAxis], y[:, NewAxis]),
    10151041                                         axis=1)
     
    10201046            I = self.interpolation_object
    10211047
    1022             if allclose(interpolation_points, I._point_coordinates):
    1023                 can_reuse = True
     1048            if allclose(interpolation_points.shape,
     1049                        I._point_coordinates.shape):
     1050                if allclose(interpolation_points, I._point_coordinates):
     1051                    can_reuse = True
    10241052               
    10251053
     
    10471075
    10481076        return X, Compatible list, Numeric array (see below)
    1049         interpolation_points: List of x, y coordinates where value is
    1050         sought (using interpolation). If points are given, values of
    1051         location and indices are ignored
    1052        
    1053         location: Where values are to be stored.
    1054                   Permissible options are: vertices, edges, centroids
    1055                   and unique vertices. Default is 'vertices'
     1077       
     1078        Inputs:
     1079           interpolation_points: List of x, y coordinates where value is
     1080                                 sought (using interpolation). If points
     1081                                 are given, values of location and indices
     1082                                 are ignored. Assume either absolute UTM
     1083                                 coordinates or geospatial data object.
     1084       
     1085           location: Where values are to be stored.
     1086                     Permissible options are: vertices, edges, centroids
     1087                     and unique vertices. Default is 'vertices'
    10561088
    10571089
     
    10751107        internal ordering.
    10761108        """
     1109       
    10771110        from Numeric import take
    10781111
     
    10861119       
    10871120        # FIXME (Ole): Consider deprecating 'edges' - but not if it is used
    1088         # elsewhere in ANUGA.
     1121        # elsewhere in ANUGA.
     1122        # Edges have already been deprecated in set_values, see changeset:5521,
     1123        # but *might* be useful in get_values. Any thoughts anyone?
     1124       
    10891125        if location not in ['vertices', 'centroids', 'edges',
    10901126                            'unique vertices']:
  • anuga_core/source/anuga/abstract_2d_finite_volumes/test_quantity.py

    r5521 r5729  
    23542354
    23552355
     2356    def test_get_interpolated_values_with_georef(self):
     2357   
     2358        zone = 56
     2359        xllcorner = 308500
     2360        yllcorner = 6189000
     2361        a = [0.0, 0.0]
     2362        b = [0.0, 2.0]
     2363        c = [2.0,0.0]
     2364        d = [0.0, 4.0]
     2365        e = [2.0, 2.0]
     2366        f = [4.0,0.0]
     2367
     2368        points = [a, b, c, d, e, f]
     2369        #bac, bce, ecf, dbe
     2370        vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
     2371
     2372        domain = Domain(points, vertices,
     2373                        geo_reference=Geo_reference(zone,xllcorner,yllcorner))
     2374
     2375        quantity = Quantity(domain)
     2376        quantity.set_values(lambda x, y: x+2*y) #2 4 4 6
     2377
     2378        #First pick one point (and turn it into absolute coordinates)
     2379        x, y = 2.0/3, 8.0/3
     2380        v = quantity.get_values(interpolation_points = [[x+xllcorner,y+yllcorner]])
     2381        assert allclose(v, 6)
     2382       
     2383
     2384        # Then another to test that algorithm won't blindly
     2385        # reuse interpolation matrix
     2386        x, y = 4.0/3, 4.0/3
     2387        v = quantity.get_values(interpolation_points = [[x+xllcorner,y+yllcorner]])
     2388        assert allclose(v, 4)       
     2389       
     2390        # Try two points
     2391        pts = [[2.0/3 + xllcorner, 8.0/3 + yllcorner],
     2392               [4.0/3 + xllcorner, 4.0/3 + yllcorner]]         
     2393        v = quantity.get_values(interpolation_points=pts)
     2394        assert allclose(v, [6, 4])               
     2395       
     2396        # Test it using the geospatial data format with absolute input points and default georef
     2397        pts = Geospatial_data(data_points=pts)
     2398        v = quantity.get_values(interpolation_points=pts)
     2399        assert allclose(v, [6, 4])                               
     2400       
     2401       
     2402        # Test it using the geospatial data format with relative input points
     2403        pts = Geospatial_data(data_points=[[2.0/3, 8.0/3], [4.0/3, 4.0/3]],
     2404                              geo_reference=Geo_reference(zone,xllcorner,yllcorner))
     2405        v = quantity.get_values(interpolation_points=pts)
     2406        assert allclose(v, [6, 4])                       
     2407       
     2408       
     2409       
    23562410
    23572411    def test_getting_some_vertex_values(self):
  • anuga_core/source/anuga/geospatial_data/geospatial_data.py

    r5421 r5729  
    15431543                    alpha = alpha)
    15441544
    1545         points_geo=domain.geo_reference.change_points_geo_ref(points)
     1545       
     1546        # Convert points to geospatial data for use with get_values below
     1547        points_geo = Geospatial_data(points, domain.geo_reference)
     1548       
    15461549        #returns the predicted elevation of the points that were "split" out
    15471550        #of the original data set for one particular alpha
  • anuga_core/source/anuga/geospatial_data/test_geospatial_data.py

    r5203 r5729  
    18921892            os.remove(filename)
    18931893           
    1894     #        print value, alpha
     1894            # print value, alpha
    18951895            assert (alpha==0.01)
    18961896
     
    19761976
    19771977    #suite = unittest.makeSuite(Test_Geospatial_data, 'test_write_csv_attributes_lat_long')
    1978 #    suite = unittest.makeSuite(Test_Geospatial_data, 'test_find_optimal_smoothing_parameter')
    1979 #    suite = unittest.makeSuite(Test_Geospatial_data, 'test_split1')
     1978    #suite = unittest.makeSuite(Test_Geospatial_data, 'test_find_optimal_smoothing_parameter')
     1979    #suite = unittest.makeSuite(Test_Geospatial_data, 'test_split1')
    19801980    suite = unittest.makeSuite(Test_Geospatial_data, 'test')
    19811981    runner = unittest.TextTestRunner() #verbosity=2)
  • anuga_core/source/anuga/shallow_water/data_manager.py

    r5726 r5729  
    100100from anuga.abstract_2d_finite_volumes.util import get_revision_number, \
    101101     remove_lone_verts, sww2timeseries, get_centroid_values
     102     
     103from anuga.abstract_2d_finite_volumes.neighbour_mesh import segment_midpoints
    102104from anuga.load_mesh.loadASCII import export_mesh_file
    103105from anuga.utilities.polygon import intersection
     
    63706372    segments = mesh.get_intersecting_segments(polyline, verbose=verbose)
    63716373   
    6372        
    6373     # Then store for each triangle the length of the intersecting segment(s),
    6374     # right hand normal(s) and midpoints as interpolation_points.
    6375     interpolation_points = []
    6376     for segment in segments:
    6377         midpoint = sum(array(segment.segment))/2
    6378         interpolation_points.append(midpoint)
    6379 
    6380 
     6374    # Get midpoints
     6375    interpolation_points = segment_midpoints(segments)       
     6376
     6377    # Interpolate
    63816378    if verbose:
    63826379        print 'Interpolating - ',       
    63836380        print 'total number of interpolation points = %d'\
    63846381              %len(interpolation_points)
    6385              
    6386              
    63876382   
    63886383    I = Interpolation_function(time,
     
    64456440    for t in time:
    64466441        total_flow=0
    6447         for i, p in enumerate(interpolation_points):
     6442        for i in range(len(interpolation_points)):
    64486443            elevation, stage, uh, vh = interpolation_function(t, point_id=i)
    64496444            normal = segments[i].normal
  • anuga_core/source/anuga/shallow_water/shallow_water_domain.py

    r5672 r5729  
    8383from Numeric import compress, arange
    8484
    85 
     85from anuga.abstract_2d_finite_volumes.neighbour_mesh import segment_midpoints
    8686from anuga.abstract_2d_finite_volumes.domain import Domain as Generic_Domain
    8787from anuga.abstract_2d_finite_volumes.generic_boundary_conditions\
     
    374374        return self.get_quantity('elevation').\
    375375               get_maximum_location(indices=wet_elements)   
     376               
     377               
     378               
     379               
     380    def get_flow_through_cross_section(self, polyline,
     381                                       verbose=False):               
     382        """Get the total flow through an arbitrary poly line.       
     383       
     384        This is a run-time equivalent of the function with same name in data_manager.py
     385       
     386        Input:
     387            polyline: Representation of desired cross section - it may contain
     388                      multiple sections allowing for complex shapes. Assume
     389                      absolute UTM coordinates.
     390                      Format [[x0, y0], [x1, y1], ...]       
     391                 
     392        Output:       
     393            Q: Total flow [m^3/s] across given segments.
     394       
     395         
     396        """       
     397       
     398        # Adjust polyline to mesh spatial origin
     399        polyline = self.geo_reference.get_relative(polyline)
     400       
     401        # Find all intersections and associated triangles.
     402        segments = self.get_intersecting_segments(polyline, verbose=verbose)
     403
     404        msg = 'No segments found'
     405        assert len(segments) > 0, msg
     406       
     407        # Get midpoints
     408        midpoints = segment_midpoints(segments)       
     409       
     410        # FIXME (Ole): HACK - need to make midpoints Geospatial instances
     411        midpoints = self.geo_reference.get_absolute(midpoints)       
     412       
     413        # Compute flow       
     414        if verbose: print 'Computing flow through specified cross section'
     415       
     416        # Get interpolated values
     417        xmomentum = self.get_quantity('xmomentum')
     418        ymomentum = self.get_quantity('ymomentum')       
     419       
     420        uh = xmomentum.get_values(interpolation_points=midpoints)
     421        vh = ymomentum.get_values(interpolation_points=midpoints)       
     422       
     423        # Compute and sum flows across each segment
     424        total_flow=0
     425        for i in range(len(uh)):
     426           
     427            # Inner product of momentum vector with segment normal [m^2/s]
     428            normal = segments[i].normal
     429            normal_momentum = uh[i]*normal[0] + vh[i]*normal[1]
     430               
     431            # Flow across this segment [m^3/s]
     432            segment_flow = normal_momentum*segments[i].length
     433
     434            # Accumulate
     435            total_flow += segment_flow
     436
     437           
     438        return total_flow
     439       
    376440
    377441    def check_integrity(self):
  • anuga_core/source/anuga/shallow_water/test_data_manager.py

    r5726 r5729  
    1007310073        u = 2 m/s
    1007410074        h = 1 m
    10075         w = 5 m (width of channel)
    10076 
    10077         q = u*h*w = 10 m^3/s
     10075        w = 3 m (width of channel)
     10076
     10077        q = u*h*w = 6 m^3/s
    1007810078
    1007910079        #---------- First run without geo referencing       
     
    1011210112
    1011310113        Br = Reflective_boundary(domain)     # Side walls
    10114         Bd = Dirichlet_boundary([h, uh, 0])  # 2 m/s across the 5 m inlet:
     10114        Bd = Dirichlet_boundary([h, uh, 0])  # 2 m/s across the 3 m inlet:
    1011510115
    1011610116
     
    1019910199        u = 2 m/s
    1020010200        h = 2 m
    10201         w = 5 m (width of channel)
    10202 
    10203         q = u*h*w = 20 m^3/s
     10201        w = 3 m (width of channel)
     10202
     10203        q = u*h*w = 12 m^3/s
    1020410204
    1020510205
     
    1024310243
    1024410244        Br = Reflective_boundary(domain)     # Side walls
    10245         Bd = Dirichlet_boundary([w, uh, 0])  # 2 m/s across the 5 m inlet:
     10245        Bd = Dirichlet_boundary([w, uh, 0])  # 2 m/s across the 3 m inlet:
    1024610246
    1024710247
     
    1030010300        u = 2 m/s
    1030110301        h = 1 m
    10302         w = 5 m (width of channel)
    10303 
    10304         q = u*h*w = 10 m^3/s
     10302        w = 3 m (width of channel)
     10303
     10304        q = u*h*w = 6 m^3/s
    1030510305        Es = h + 0.5*v*v/g  # Specific energy head [m]
    1030610306        Et = w + 0.5*v*v/g  # Total energy head [m]       
     
    1034610346
    1034710347        Br = Reflective_boundary(domain)     # Side walls
    10348         Bd = Dirichlet_boundary([w, uh, 0])  # 2 m/s across the 5 m inlet:
     10348        Bd = Dirichlet_boundary([w, uh, 0])  # 2 m/s across the 3 m inlet:
    1034910349
    1035010350       
  • anuga_core/source/anuga/shallow_water/test_shallow_water_domain.py

    r5450 r5729  
    14061406           
    14071407       
     1408       
     1409    def test_get_flow_through_cross_section_with_geo(self):
     1410        """test_get_flow_through_cross_section(self):
     1411
     1412        Test that the total flow through a cross section can be
     1413        correctly obtained at run-time from the ANUGA domain.
     1414       
     1415        This test creates a flat bed with a known flow through it and tests
     1416        that the function correctly returns the expected flow.
     1417
     1418        The specifics are
     1419        e = -1 m
     1420        u = 2 m/s
     1421        h = 2 m
     1422        w = 3 m (width of channel)
     1423
     1424        q = u*h*w = 12 m^3/s
     1425
     1426
     1427        This run tries it with georeferencing and with elevation = -1
     1428       
     1429        """
     1430
     1431        import time, os
     1432        from Numeric import array, zeros, allclose, Float, concatenate
     1433        from Scientific.IO.NetCDF import NetCDFFile
     1434
     1435        # Setup
     1436        from mesh_factory import rectangular
     1437
     1438        # Create basic mesh (20m x 3m)
     1439        width = 3
     1440        length = 20
     1441        t_end = 1
     1442        points, vertices, boundary = rectangular(length, width,
     1443                                                 length, width)
     1444
     1445        # Create shallow water domain
     1446        domain = Domain(points, vertices, boundary,
     1447                        geo_reference=Geo_reference(56,308500,6189000))
     1448
     1449        domain.default_order = 2
     1450        domain.set_quantities_to_be_stored(None)
     1451
     1452
     1453        e = -1.0
     1454        w = 1.0
     1455        h = w-e
     1456        u = 2.0
     1457        uh = u*h
     1458
     1459        Br = Reflective_boundary(domain)     # Side walls
     1460        Bd = Dirichlet_boundary([w, uh, 0])  # 2 m/s across the 3 m inlet:
     1461
     1462
     1463        # Initial conditions
     1464        domain.set_quantity('elevation', e)
     1465        domain.set_quantity('stage', w)
     1466        domain.set_quantity('xmomentum', uh)
     1467        domain.set_boundary( {'left': Bd, 'right': Bd, 'top': Br, 'bottom': Br})
     1468       
     1469       
     1470        # Interpolation points down the middle
     1471        I = [[0, width/2.],
     1472             [length/2., width/2.],
     1473             [length, width/2.]]
     1474        interpolation_points = domain.geo_reference.get_absolute(I)       
     1475       
     1476        # Shortcuts to quantites
     1477        stage = domain.get_quantity('stage')       
     1478        xmomentum = domain.get_quantity('xmomentum')       
     1479        ymomentum = domain.get_quantity('ymomentum')               
     1480
     1481        for t in domain.evolve(yieldstep=0.1, finaltime = t_end):
     1482            # Check that quantities are they should be in the interior
     1483
     1484            w_t = stage.get_values(interpolation_points)           
     1485            uh_t = xmomentum.get_values(interpolation_points)
     1486            vh_t = ymomentum.get_values(interpolation_points)           
     1487           
     1488            assert allclose(w_t, w)
     1489            assert allclose(uh_t, uh)           
     1490            assert allclose(vh_t, 0.0)                       
     1491           
     1492           
     1493            # Check flows through the middle
     1494            for i in range(5):
     1495                x = length/2. + i*0.23674563 # Arbitrary
     1496                cross_section = [[x, 0], [x, width]]
     1497
     1498                cross_section = domain.geo_reference.get_absolute(cross_section)           
     1499                Q = domain.get_flow_through_cross_section(cross_section,
     1500                                                          verbose=False)
     1501
     1502                assert allclose(Q, uh*width)
     1503
     1504
     1505           
     1506       
     1507       
    14081508
    14091509    def test_another_runup_example(self):
     
    14311531        domain = Domain(points, vertices, boundary) # Create domain
    14321532        domain.set_quantities_to_be_stored(None)
    1433         domain.set_maximum_allowed_speed(100) #
     1533        domain.set_maximum_allowed_speed(100) #FIXME (Ole): try to remove this
    14341534       
    14351535        # FIXME (Ole): Need tests where this is commented out
    1436         domain.tight_slope_limiters = 0 # Backwards compatibility (14/4/7)                 
     1536        domain.tight_slope_limiters = 0 # Backwards compatibility (14/4/7)
    14371537        domain.H0 = 0 # Backwards compatibility (6/2/7)
    14381538        domain.use_centroid_velocities = 0 # Backwards compatibility (7/5/8)
     
    58055905if __name__ == "__main__":
    58065906
    5807     suite = unittest.makeSuite(Test_Shallow_Water,'test')
    5808 
    5809     #suite = unittest.makeSuite(Test_Shallow_Water,'test_bedslope_problem_first_order_moresteps')   
     5907    #suite = unittest.makeSuite(Test_Shallow_Water,'test')
     5908
     5909    suite = unittest.makeSuite(Test_Shallow_Water,'test_get_flow_through_cross_section_with_geo')   
    58105910    #suite = unittest.makeSuite(Test_Shallow_Water,'test_fitting_using_shallow_water_domain')   
    58115911    #suite = unittest.makeSuite(Test_Shallow_Water,'test_tight_slope_limiters')
  • anuga_core/source/anuga/utilities/numerical_tools.py

    r5561 r5729  
    313313
    314314    return bins
     315   
     316
    315317
    316318def get_machine_precision():
Note: See TracChangeset for help on using the changeset viewer.