Changeset 5726


Ignore:
Timestamp:
Sep 3, 2008, 4:44:23 PM (14 years ago)
Author:
ole
Message:

Energy through cross section and simple tests done for sww file input.
This is work towards ticket:295

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

Legend:

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

    r5723 r5726  
    64666466   
    64676467
     6468def get_energy_through_cross_section(filename,
     6469                                     polyline,
     6470                                     kind = 'total',
     6471                                     verbose=False):
     6472    """Obtain flow (m^3/s) perpendicular to specified cross section.
     6473
     6474    Inputs:
     6475        filename: Name of sww file
     6476        polyline: Representation of desired cross section - it may contain
     6477                  multiple sections allowing for complex shapes. Assume
     6478                  absolute UTM coordinates.
     6479                  Format [[x0, y0], [x1, y1], ...]
     6480        kind:     Select which energy to compute.
     6481                  Options are 'specific' and 'total' (default)
     6482
     6483    Output:
     6484        time: All stored times in sww file
     6485        Q: Average energy [m] across given segments for all stored times.
     6486
     6487    The average velocity is computed for each triangle intersected by the polyline
     6488    and averaged weigted by segment lengths.
     6489
     6490    The typical usage of this function would be to get average energy of flow in a channel,
     6491    and the polyline would then be a cross section perpendicular to the flow.
     6492
     6493   
     6494    #FIXME (Ole) - need name for this energy reflecting that its dimension is [m].
     6495    """
     6496
     6497    from anuga.config import g, epsilon, velocity_protection as h0       
     6498           
     6499    quantity_names =['elevation',
     6500                     'stage',
     6501                     'xmomentum',
     6502                     'ymomentum']
     6503
     6504
     6505    # Get values for quantities at each midpoint of poly line from sww file
     6506    X = get_interpolated_quantities_at_polyline_midpoints(filename,
     6507                                                          quantity_names=quantity_names,
     6508                                                          polyline=polyline,
     6509                                                          verbose=verbose)   
     6510    segments, interpolation_function = X
     6511
     6512   
     6513    # Get vectors for time and interpolation_points
     6514    time = interpolation_function.time
     6515    interpolation_points = interpolation_function.interpolation_points   
     6516
     6517    if verbose: print 'Computing %s energy' %kind
     6518   
     6519    # Compute total length of polyline for use with weighted averages
     6520    total_line_length = 0.0
     6521    for segment in segments:
     6522        total_line_length += segment.length
     6523       
     6524    # Compute energy
     6525    E = []
     6526    for t in time:
     6527        average_energy=0.0
     6528        for i, p in enumerate(interpolation_points):
     6529            elevation, stage, uh, vh = interpolation_function(t, point_id=i)
     6530           
     6531            # Depth
     6532            h = depth = stage-elevation
     6533           
     6534            # Average velocity across this segment
     6535            if h > epsilon:
     6536                # Use protection against degenerate velocities
     6537                u = uh/(h + h0/h)
     6538                v = vh/(h + h0/h)
     6539            else:
     6540                u = v = 0.0
     6541               
     6542            speed_squared = u*u + v*v   
     6543            kinetic_energy = 0.5*speed_squared/g
     6544           
     6545            if kind == 'specific':
     6546                segment_energy = depth + kinetic_energy
     6547            elif kind == 'total':
     6548                segment_energy = stage + kinetic_energy               
     6549            else:
     6550                msg = 'Energy kind must be either "specific" or "total".'
     6551                msg += ' I got %s' %kind
     6552               
     6553
     6554            # Add to weighted average
     6555            weigth = segments[i].length/total_line_length
     6556            average_energy += segment_energy*weigth
     6557             
     6558
     6559        # Store energy at this timestep   
     6560        E.append(average_energy)
     6561
     6562
     6563    return time, E
     6564   
     6565
    64686566
    64696567
  • anuga_core/source/anuga/shallow_water/test_data_manager.py

    r5672 r5726  
    1717from anuga.shallow_water import *
    1818from anuga.shallow_water.data_manager import *
    19 from anuga.config import epsilon
     19from anuga.config import epsilon, g
    2020from anuga.utilities.anuga_exceptions import ANUGAError
    2121from anuga.utilities.numerical_tools import ensure_numeric
     
    1019810198        The specifics are
    1019910199        u = 2 m/s
    10200         h = 1 m
     10200        h = 2 m
    1020110201        w = 5 m (width of channel)
    1020210202
    10203         q = u*h*w = 10 m^3/s
    10204 
    10205 
    10206         This run tries it with georeferencing
     10203        q = u*h*w = 20 m^3/s
     10204
     10205
     10206        This run tries it with georeferencing and with elevation = -1
    1020710207       
    1020810208        """
     
    1023610236        domain.smooth = True
    1023710237
    10238         h = 1.0
     10238        e = -1.0
     10239        w = 1.0
     10240        h = w-e
    1023910241        u = 2.0
    1024010242        uh = u*h
    1024110243
    1024210244        Br = Reflective_boundary(domain)     # Side walls
    10243         Bd = Dirichlet_boundary([h, uh, 0])  # 2 m/s across the 5 m inlet:
    10244 
    10245 
    10246 
    10247        
    10248         domain.set_quantity('elevation', 0.0)
    10249         domain.set_quantity('stage', h)
     10245        Bd = Dirichlet_boundary([w, uh, 0])  # 2 m/s across the 5 m inlet:
     10246
     10247
     10248
     10249       
     10250        domain.set_quantity('elevation', e)
     10251        domain.set_quantity('stage', w)
    1025010252        domain.set_quantity('xmomentum', uh)
    1025110253        domain.set_boundary( {'left': Bd, 'right': Bd, 'top': Br, 'bottom': Br})
     
    1026810270        for t in range(t_end+1):
    1026910271            for i in range(3):
    10270                 assert allclose(f(t, i), [1, 2, 0])
     10272                #print i, t, f(t, i)           
     10273                assert allclose(f(t, i), [w, uh, 0])
    1027110274           
    1027210275
     
    1028410287
    1028510288
     10289           
     10290    def test_get_energy_through_cross_section(self):
     10291        """test_get_energy_through_cross_section(self):
     10292
     10293        Test that the specific and total energy through a cross section can be
     10294        correctly obtained from an sww file.
     10295       
     10296        This test creates a flat bed with a known flow through it and tests
     10297        that the function correctly returns the expected energies.
     10298
     10299        The specifics are
     10300        u = 2 m/s
     10301        h = 1 m
     10302        w = 5 m (width of channel)
     10303
     10304        q = u*h*w = 10 m^3/s
     10305        Es = h + 0.5*v*v/g  # Specific energy head [m]
     10306        Et = w + 0.5*v*v/g  # Total energy head [m]       
     10307
     10308
     10309        This test uses georeferencing
     10310       
     10311        """
     10312
     10313        import time, os
     10314        from Numeric import array, zeros, allclose, Float, concatenate
     10315        from Scientific.IO.NetCDF import NetCDFFile
     10316
     10317        # Setup
     10318        from mesh_factory import rectangular
     10319
     10320        # Create basic mesh (20m x 3m)
     10321        width = 3
     10322        length = 20
     10323        t_end = 1
     10324        points, vertices, boundary = rectangular(length, width,
     10325                                                 length, width)
     10326
     10327        # Create shallow water domain
     10328        domain = Domain(points, vertices, boundary,
     10329                        geo_reference = Geo_reference(56,308500,6189000))
     10330
     10331        domain.default_order = 2
     10332        domain.set_minimum_storable_height(0.01)
     10333
     10334        domain.set_name('flowtest')
     10335        swwfile = domain.get_name() + '.sww'
     10336
     10337        domain.set_datadir('.')
     10338        domain.format = 'sww'
     10339        domain.smooth = True
     10340
     10341        e = -1.0
     10342        w = 1.0
     10343        h = w-e
     10344        u = 2.0
     10345        uh = u*h
     10346
     10347        Br = Reflective_boundary(domain)     # Side walls
     10348        Bd = Dirichlet_boundary([w, uh, 0])  # 2 m/s across the 5 m inlet:
     10349
     10350       
     10351        domain.set_quantity('elevation', e)
     10352        domain.set_quantity('stage', w)
     10353        domain.set_quantity('xmomentum', uh)
     10354        domain.set_boundary( {'left': Bd, 'right': Bd, 'top': Br, 'bottom': Br})
     10355
     10356        for t in domain.evolve(yieldstep=1, finaltime = t_end):
     10357            pass
     10358
     10359        # Check that momentum is as it should be in the interior
     10360
     10361        I = [[0, width/2.],
     10362             [length/2., width/2.],
     10363             [length, width/2.]]
     10364       
     10365        I = domain.geo_reference.get_absolute(I)
     10366        f = file_function(swwfile,
     10367                          quantities=['stage', 'xmomentum', 'ymomentum'],
     10368                          interpolation_points=I,
     10369                          verbose=False)
     10370
     10371        for t in range(t_end+1):
     10372            for i in range(3):
     10373                #print i, t, f(t, i)
     10374                assert allclose(f(t, i), [w, uh, 0])
     10375           
     10376
     10377        # Check energies through the middle
     10378        for i in range(5):
     10379            x = length/2. + i*0.23674563 # Arbitrary
     10380            cross_section = [[x, 0], [x, width]]
     10381
     10382            cross_section = domain.geo_reference.get_absolute(cross_section)           
     10383           
     10384            time, Es = get_energy_through_cross_section(swwfile,
     10385                                                       cross_section,
     10386                                                       kind='specific',
     10387                                                       verbose=False)
     10388            assert allclose(Es, h + 0.5*u*u/g)
     10389           
     10390            time, Et = get_energy_through_cross_section(swwfile,
     10391                                                        cross_section,
     10392                                                        kind='total',
     10393                                                        verbose=False)
     10394            assert allclose(Et, w + 0.5*u*u/g)           
     10395
     10396           
    1028610397       
    1028710398    def test_get_all_swwfiles(self):
Note: See TracChangeset for help on using the changeset viewer.