Changeset 5226


Ignore:
Timestamp:
Apr 21, 2008, 5:21:26 PM (15 years ago)
Author:
ole
Message:

Work done during Water Down Under 2008.
Started algorithm for flow through a cross section as per ticket:175 (test disabled).

Location:
anuga_core/source/anuga/shallow_water
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • anuga_core/source/anuga/shallow_water/data_manager.py

    r5189 r5226  
    9898     remove_lone_verts, sww2timeseries, get_centroid_values
    9999from anuga.load_mesh.loadASCII import export_mesh_file
     100from anuga.utilities.polygon import intersection
     101
     102
    100103# formula mappings
    101104
     
    54925495
    54935496
     5497# ----------------------------------------------
     5498# Functions to obtain diagnostics from sww files
     5499#-----------------------------------------------
     5500
     5501def get_mesh_and_quantities_from_sww_file(filename, quantity_names, verbose=False):
     5502    """Get and rebuild mesh structure and the associated quantities from sww file
     5503    """
     5504
     5505    #FIXME(Ole): This is work in progress
     5506   
     5507    import types
     5508    # FIXME (Ole): Maybe refactor filefunction using this more fundamental code.
     5509
     5510    return
     5511   
     5512    # Open NetCDF file
     5513    if verbose: print 'Reading', filename
     5514    fid = NetCDFFile(filename, 'r')
     5515
     5516    if type(quantity_names) == types.StringType:
     5517        quantity_names = [quantity_names]       
     5518
     5519    if quantity_names is None or len(quantity_names) < 1:
     5520        msg = 'No quantities are specified'
     5521        raise Exception, msg
     5522
     5523    # Now assert that requested quantitites (and the independent ones)
     5524    # are present in file
     5525    missing = []
     5526    for quantity in ['x', 'y', 'volumes', 'time'] + quantity_names:
     5527        if not fid.variables.has_key(quantity):
     5528            missing.append(quantity)
     5529
     5530    if len(missing) > 0:
     5531        msg = 'Quantities %s could not be found in file %s'\
     5532              %(str(missing), filename)
     5533        fid.close()
     5534        raise Exception, msg
     5535
     5536    if not filename.endswith('.sww'):
     5537        msg = 'Filename must have extension .sww'       
     5538        raise Exception, msg       
     5539
     5540    # Get first timestep
     5541    try:
     5542        starttime = fid.starttime[0]
     5543    except ValueError:
     5544        msg = 'Could not read starttime from file %s' %filename
     5545        raise msg
     5546
     5547    # Get variables
     5548    time = fid.variables['time'][:]   
     5549
     5550    # Get origin
     5551    xllcorner = fid.xllcorner[0]
     5552    yllcorner = fid.yllcorner[0]
     5553    zone = fid.zone[0]
     5554    georeference = Geo_reference(zone, xllcorner, yllcorner)
     5555
     5556
     5557    x = fid.variables['x'][:]
     5558    y = fid.variables['y'][:]
     5559    triangles = fid.variables['volumes'][:]
     5560
     5561    x = reshape(x, (len(x),1))
     5562    y = reshape(y, (len(y),1))
     5563    vertex_coordinates = concatenate((x,y), axis=1) #m x 2 array
     5564
     5565    #if interpolation_points is not None:
     5566    #    # Adjust for georef
     5567    #    interpolation_points[:,0] -= xllcorner
     5568    #    interpolation_points[:,1] -= yllcorner       
     5569       
     5570
     5571    # Produce values for desired data points at
     5572    # each timestep for each quantity
     5573    quantities = {}
     5574    for name in quantity_names:
     5575        quantities[name] = fid.variables[name][:]
     5576       
     5577    fid.close()
     5578
     5579    # Create mesh and quad tree
     5580    #interpolator = Interpolate(vertex_coordinates, triangles)
     5581
     5582    #return interpolator, quantities, geo_reference, time
     5583
     5584
     5585def get_flow_through_cross_section(filename,
     5586                                   polyline,
     5587                                   verbose=False):
     5588    """Obtain flow (m^3/s) perpendicular to cross section given by the argument polyline.
     5589
     5590    Inputs:
     5591        filename: Name of sww file
     5592        polyline: Representation of desired cross section - it may contain multiple
     5593                  sections allowing for complex shapes.
     5594
     5595    Output:
     5596        Q: Hydrograph of total flow across given segments for all stored timesteps.
     5597
     5598    The normal flow is computed for each triangle intersected by the polyline and added up.
     5599    If multiple sections are specified normal flows may partially cancel each other.
     5600
     5601    """
     5602
     5603    # Get mesh and quantities from sww file
     5604    X = get_mesh_and_quantities_from_sww_file(filename, ['elevation',
     5605                                                         'stage',
     5606                                                         'xmomentum',
     5607                                                         'ymomentum'], verbose=verbose)
     5608    interpolator, quantities, geo_reference, time = X
     5609
     5610
     5611   
     5612    # Find all intersections and associated triangles.
     5613   
     5614    get_intersecting_segments(polyline)
     5615   
     5616    # Then store for each triangle the length of the intersecting segment(s),
     5617    # right hand normal(s) and midpoints.
     5618    pass
     5619
    54945620
    54955621
  • anuga_core/source/anuga/shallow_water/test_data_manager.py

    r5189 r5226  
    1919from anuga.utilities.numerical_tools import ensure_numeric
    2020from anuga.coordinate_transforms.redfearn import degminsec2decimal_degrees
     21from anuga.abstract_2d_finite_volumes.util import file_function
    2122
    2223# This is needed to run the tests of local functions
     
    73067307
    73077308
    7308         #Cleanup
     7309        # Cleanup
    73097310        os.remove(swwfile)
     7311
     7312
     7313    def NOtest_get_flow_through_cross_section(self):
     7314        """test_get_flow_through_cross_section(self):
     7315
     7316        Test that the total flow through a cross section can be
     7317        correctly obtained from an sww file.
     7318       
     7319        This test creates a flat bed with a known flow through it and tests
     7320        that the function correctly returns the expected flow.
     7321
     7322        The specifics are
     7323        u = 2 m/s
     7324        h = 1 m
     7325        w = 5 m (width of channel)
     7326
     7327        q = u*h*w = 10 m^3/s
     7328       
     7329        """
     7330
     7331        import time, os
     7332        from Numeric import array, zeros, allclose, Float, concatenate
     7333        from Scientific.IO.NetCDF import NetCDFFile
     7334
     7335        # Setup
     7336        from mesh_factory import rectangular
     7337
     7338        # Create basic mesh (100m x 5m)
     7339        width = 5
     7340        len = 100
     7341        points, vertices, boundary = rectangular(len, width, 100, 5)
     7342
     7343        # Create shallow water domain
     7344        domain = Domain(points, vertices, boundary)
     7345        domain.default_order = 2
     7346        domain.set_minimum_storable_height(0.01)
     7347
     7348        domain.set_name('flowtest')
     7349        swwfile = domain.get_name() + '.sww'
     7350
     7351        domain.set_datadir('.')
     7352        domain.format = 'sww'
     7353        domain.smooth = True
     7354
     7355        h = 1.0
     7356        u = 2.0
     7357        uh = u*h
     7358
     7359        Br = Reflective_boundary(domain)     # Side walls
     7360        Bd = Dirichlet_boundary([h, uh, 0])  # 2 m/s across the 5 m inlet:
     7361
     7362
     7363        #---------- First run without geo referencing
     7364       
     7365        domain.set_quantity('elevation', 0.0)
     7366        domain.set_quantity('stage', h)
     7367        domain.set_quantity('xmomentum', uh)
     7368        domain.set_boundary( {'left': Bd, 'right': Bd, 'top': Br, 'bottom': Br})
     7369
     7370        for t in domain.evolve(yieldstep=1, finaltime = 50):
     7371            pass
     7372
     7373        # Check that momentum is as it should be in the interior
     7374        f = file_function(swwfile,
     7375                          quantities=['stage', 'xmomentum', 'ymomentum'],
     7376                          interpolation_points=[[0,width/2],
     7377                                                [len/2, width/2],
     7378                                                [len, width/2]],
     7379                          verbose=False)
     7380
     7381        for t in range(50):
     7382            for i in range(3):
     7383                assert allclose(f(t, i), [1, 2, 0])
     7384           
     7385
     7386
     7387        # Check flow through the middle
     7388        cross_section = [[len/2,0], [len/2,width]]
     7389        Q = get_flow_through_cross_section(swwfile,
     7390                                           cross_section,
     7391                                           verbose=True)
     7392
     7393        assert allclose(Q, uh*width) 
     7394                                     
     7395
     7396
     7397
     7398
     7399       
    73107400       
    73117401    def test_get_all_swwfiles(self):
     
    74177507#    suite = unittest.makeSuite(Test_Data_Manager,'test_screen_catcher')
    74187508    suite = unittest.makeSuite(Test_Data_Manager,'test')
    7419     #suite = unittest.makeSuite(Test_Data_Manager,'test_urs_ungridded_holeII')
     7509    #suite = unittest.makeSuite(Test_Data_Manager,'test_get_flow_through_cross_section')
    74207510    #suite = unittest.makeSuite(Test_Data_Manager,'test_urs_ungridded_holeII')
    74217511
Note: See TracChangeset for help on using the changeset viewer.