Changeset 5226

Timestamp:

Apr 21, 2008, 5:21:26 PM (17 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:

• data_manager.py (modified) (2 diffs)
• test_data_manager.py (modified) (3 diffs)

anuga_core/source/anuga/shallow_water/data_manager.py

@@ -98 +98 @@
      remove_lone_verts, sww2timeseries, get_centroid_values
 from anuga.load_mesh.loadASCII import export_mesh_file
+from anuga.utilities.polygon import intersection
+
+
 # formula mappings
 
@@ -5492 +5495 @@
 
 
+# ----------------------------------------------
+# Functions to obtain diagnostics from sww files
+#-----------------------------------------------
+
+def get_mesh_and_quantities_from_sww_file(filename, quantity_names, verbose=False):
+    """Get and rebuild mesh structure and the associated quantities from sww file
+    """
+
+    #FIXME(Ole): This is work in progress
+    
+    import types
+    # FIXME (Ole): Maybe refactor filefunction using this more fundamental code.
+
+    return
+    
+    # Open NetCDF file
+    if verbose: print 'Reading', filename
+    fid = NetCDFFile(filename, 'r')
+
+    if type(quantity_names) == types.StringType:
+        quantity_names = [quantity_names]        
+
+    if quantity_names is None or len(quantity_names) < 1:
+        msg = 'No quantities are specified'
+        raise Exception, msg
+
+    # Now assert that requested quantitites (and the independent ones)
+    # are present in file 
+    missing = []
+    for quantity in ['x', 'y', 'volumes', 'time'] + quantity_names:
+        if not fid.variables.has_key(quantity):
+            missing.append(quantity)
+
+    if len(missing) > 0:
+        msg = 'Quantities %s could not be found in file %s'\
+              %(str(missing), filename)
+        fid.close()
+        raise Exception, msg
+
+    if not filename.endswith('.sww'):
+        msg = 'Filename must have extension .sww'        
+        raise Exception, msg        
+
+    # Get first timestep
+    try:
+        starttime = fid.starttime[0]
+    except ValueError:
+        msg = 'Could not read starttime from file %s' %filename
+        raise msg
+
+    # Get variables
+    time = fid.variables['time'][:]    
+
+    # Get origin
+    xllcorner = fid.xllcorner[0]
+    yllcorner = fid.yllcorner[0]
+    zone = fid.zone[0]
+    georeference = Geo_reference(zone, xllcorner, yllcorner)
+
+
+    x = fid.variables['x'][:]
+    y = fid.variables['y'][:]
+    triangles = fid.variables['volumes'][:]
+
+    x = reshape(x, (len(x),1))
+    y = reshape(y, (len(y),1))
+    vertex_coordinates = concatenate((x,y), axis=1) #m x 2 array
+
+    #if interpolation_points is not None:
+    #    # Adjust for georef
+    #    interpolation_points[:,0] -= xllcorner
+    #    interpolation_points[:,1] -= yllcorner        
+        
+
+    # Produce values for desired data points at
+    # each timestep for each quantity
+    quantities = {}
+    for name in quantity_names:
+        quantities[name] = fid.variables[name][:]
+        
+    fid.close()
+
+    # Create mesh and quad tree
+    #interpolator = Interpolate(vertex_coordinates, triangles)
+
+    #return interpolator, quantities, geo_reference, time
+
+
+def get_flow_through_cross_section(filename,
+                                   polyline,
+                                   verbose=False):
+    """Obtain flow (m^3/s) perpendicular to cross section given by the argument polyline.
+
+    Inputs:
+        filename: Name of sww file
+        polyline: Representation of desired cross section - it may contain multiple
+                  sections allowing for complex shapes.
+
+    Output:
+        Q: Hydrograph of total flow across given segments for all stored timesteps.
+
+    The normal flow is computed for each triangle intersected by the polyline and added up.
+    If multiple sections are specified normal flows may partially cancel each other.
+
+    """
+
+    # Get mesh and quantities from sww file
+    X = get_mesh_and_quantities_from_sww_file(filename, ['elevation',
+                                                         'stage',
+                                                         'xmomentum',
+                                                         'ymomentum'], verbose=verbose)
+    interpolator, quantities, geo_reference, time = X
+
+
+    
+    # Find all intersections and associated triangles.
+    
+    get_intersecting_segments(polyline)
+    
+    # Then store for each triangle the length of the intersecting segment(s),
+    # right hand normal(s) and midpoints. 
+    pass
+
 
 

anuga_core/source/anuga/shallow_water/test_data_manager.py

@@ -19 +19 @@
 from anuga.utilities.numerical_tools import ensure_numeric
 from anuga.coordinate_transforms.redfearn import degminsec2decimal_degrees
+from anuga.abstract_2d_finite_volumes.util import file_function
 
 # This is needed to run the tests of local functions
@@ -7306 +7307 @@
 
 
-        #Cleanup
+        # Cleanup
         os.remove(swwfile)
+
+
+    def NOtest_get_flow_through_cross_section(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 
+        
+        """
+
+        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
+        len = 100
+        points, vertices, boundary = rectangular(len, width, 100, 5)
+
+        # Create shallow water domain
+        domain = Domain(points, vertices, boundary)
+        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: 
+
+
+        #---------- First run without geo referencing
+        
+        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 = 50):
+            pass
+
+        # Check that momentum is as it should be in the interior
+        f = file_function(swwfile,
+                          quantities=['stage', 'xmomentum', 'ymomentum'],
+                          interpolation_points=[[0,width/2],
+                                                [len/2, width/2],
+                                                [len, width/2]],
+                          verbose=False)
+
+        for t in range(50):
+            for i in range(3):
+                assert allclose(f(t, i), [1, 2, 0])
+            
+
+
+        # Check flow through the middle
+        cross_section = [[len/2,0], [len/2,width]]
+        Q = get_flow_through_cross_section(swwfile,
+                                           cross_section,
+                                           verbose=True)
+
+        assert allclose(Q, uh*width)  
+                                      
+
+
+
+
+        
         
     def test_get_all_swwfiles(self):
@@ -7417 +7507 @@
 #    suite = unittest.makeSuite(Test_Data_Manager,'test_screen_catcher')
     suite = unittest.makeSuite(Test_Data_Manager,'test')
-    #suite = unittest.makeSuite(Test_Data_Manager,'test_urs_ungridded_holeII')
+    #suite = unittest.makeSuite(Test_Data_Manager,'test_get_flow_through_cross_section')
     #suite = unittest.makeSuite(Test_Data_Manager,'test_urs_ungridded_holeII')