Changeset 4551

Timestamp:

Jun 20, 2007, 6:34:00 PM (17 years ago)

Author:

ole

Message:

Introduced maximal_inundation_elevation (runup) for sww files

Location:

anuga_core

Files:

• documentation/user_manual/demos/runup.py (modified) (1 diff)
• source/anuga/shallow_water/data_manager.py (modified) (5 diffs)
• source/anuga/shallow_water/test_shallow_water_domain.py (modified) (3 diffs)
• source/anuga/utilities/numerical_tools.py (modified) (1 diff)
• source/anuga/utilities/test_numerical_tools.py (modified) (2 diffs)

anuga_core/documentation/user_manual/demos/runup.py

@@ -25 +25 @@
 
 domain = Domain(points, vertices, boundary) # Create domain
-domain.set_name('runup')                    # Output to bedslope.sww
+domain.set_name('runup')                    # Output to file runup.sww
 domain.set_datadir('.')                     # Use current directory for output
 #domain.set_quantities_to_be_stored(['stage',# Store all conserved quantities

anuga_core/source/anuga/shallow_water/data_manager.py

@@ -69 +69 @@
 from Numeric import concatenate, array, Float, Int, Int32, resize, sometrue, \
      searchsorted, zeros, allclose, around, reshape, transpose, sort, \
-     NewAxis, ArrayType
+     NewAxis, ArrayType, compress, take, arange, argmax
 from Scientific.IO.NetCDF import NetCDFFile
 #from shutil import copy
@@ -1860 +1860 @@
     #if verbose: bye= nsuadsfd[0] # uncomment to check catching verbose errors
     
-    #Read sww file
+    # Read sww file
     if verbose: 
         print 'Reading from %s' %swwfile
@@ -1882 +1882 @@
     if origin is None:
 
-        #Get geo_reference
-        #sww files don't have to have a geo_ref
+        # Get geo_reference
+        # sww files don't have to have a geo_ref
         try:
             geo_reference = Geo_reference(NetCDFObject=fid)
         except AttributeError, e:
-            geo_reference = Geo_reference() #Default georef object
+            geo_reference = Geo_reference() # Default georef object
 
         xllcorner = geo_reference.get_xllcorner()
@@ -1899 +1899 @@
 
 
-    #FIXME: Refactor using code from Interpolation_function.statistics (in interpolate.py)
-    #Something like print swwstats(swwname)
+    # FIXME: Refactor using code from Interpolation_function.statistics (in interpolate.py)
+    # Something like print swwstats(swwname)
     if verbose:
         print '------------------------------------------------'
@@ -5379 +5379 @@
 
 
+
+
+
+def get_maximum_inundation_elevation(sww_filename, region=None, timesteps=None, verbose=False):
+    """Compute maximum run up height from sww file.
+
+    Algorithm is as in get_maximum_inundation_elevation from shallow_water_domain
+    except that this function works with the sww file and computes the maximal
+    runup height over multiple timesteps. 
+    
+    Optional argument region restricts this to specified polygon.
+    Optional argument timesteps restricts this to a specified time step (an index)
+    or time interval (a list of indices).
+    """
+
+    # We are using nodal values here as that is what is stored in sww files.
+
+    # Water depth below which it is considered to be 0 in the model
+    # FIXME (Ole): Allow this to be specified as a keyword argument as well
+    
+    from anuga.config import minimum_allowed_height
+
+    if region is not None:
+        msg = 'region not yet implemented in get_maximum_inundation_elevation'
+        raise Exception, msg
+
+
+    root, extension = os.path.splitext(sww_filename)
+    if extension == '':
+        sww_filename += '.sww'
+    
+    # Read sww file
+    if verbose: 
+        print 'Reading from %s' %sww_filename
+    
+    from Scientific.IO.NetCDF import NetCDFFile
+    fid = NetCDFFile(sww_filename)
+
+    # Get extent and reference
+    x = fid.variables['x'][:]
+    y = fid.variables['y'][:]
+    volumes = fid.variables['volumes'][:]
+
+
+    time = fid.variables['time'][:]
+    if timesteps is not None:
+        if type(timesteps) is list or type(timesteps) is ArrayType:  
+            timesteps = ensure_numeric(timesteps, Int)
+        elif type(timesteps) is type(0):
+            timesteps = [timesteps]
+        else:
+            msg = 'timesteps must be either an integer or a sequence of integers. '
+            msg += 'I got timesteps==%s, %s' %(timesteps, type(timesteps))
+            raise Exception, msg
+    else:
+        # Take them all
+        timesteps = arange(len(time))    
+    
+
+    # Get the relevant quantities
+    elevation = fid.variables['elevation'][:] 
+    stage = fid.variables['stage'][:]
+
+
+    # Compute maximal runup for each timestep
+    maximal_runup = None
+    for i in timesteps:
+        #print 'time', time[i]
+        depth = stage[i,:] - elevation[:]
+    
+        # Get wet nodes i.e. nodes with depth>0 within given region and timesteps
+        wet_nodes = compress(depth > minimum_allowed_height, arange(len(depth)))
+
+        # Find maximum elevation among wet nodes and return
+        #runup_index = argmax(take(elevation, wet_nodes)) #FIXME Maybe get loc as well
+        runup = max(take(elevation, wet_nodes))
+
+        if runup > maximal_runup:
+            maximal_runup = runup  # This works even if maximal_runup is None
+                
+
+    return maximal_runup
+
+
 #-------------------------------------------------------------
 if __name__ == "__main__":

anuga_core/source/anuga/shallow_water/test_shallow_water_domain.py

@@ -814 +814 @@
         """test_get_maximum_inundation_3(self)
 
-        Test real runup example
+        Test of real runup example:
         """
 
@@ -946 +946 @@
         assert allclose(wet_elevations, z)        
         
+
+
+    def test_get_maximum_inundation_from_sww(self):
+        """test_get_maximum_inundation_from_sww(self)
+
+        Test of get_maximum_inundation_elevation(sww_filename) from data_manager.py
+        
+        This is based on test_get_maximum_inundation_3(self) but works with the
+        stored results instead of with the internal data structure.
+        
+        """
+
+        from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular_cross
+        from data_manager import get_maximum_inundation_elevation
+        
+
+        initial_runup_height = -0.4
+        final_runup_height = -0.3
+
+
+        #--------------------------------------------------------------
+        # Setup computational domain
+        #--------------------------------------------------------------
+        N = 10
+        points, vertices, boundary = rectangular_cross(N, N) 
+        domain = Domain(points, vertices, boundary)
+        domain.set_name('runup_test')
+        #domain.limit2007 = 1 #FIXME: This works better with old limiters
+
+        #--------------------------------------------------------------
+        # Setup initial conditions
+        #--------------------------------------------------------------
+        def topography(x,y):
+            return -x/2                             # linear bed slope
+            
+
+        domain.set_quantity('elevation', topography)       # Use function for elevation
+        domain.set_quantity('friction', 0.)                # Zero friction 
+        domain.set_quantity('stage', initial_runup_height) # Constant negative initial stage
+
+
+        #--------------------------------------------------------------
+        # Setup boundary conditions
+        #--------------------------------------------------------------
+        Br = Reflective_boundary(domain)              # Reflective wall
+        Bd = Dirichlet_boundary([final_runup_height,  # Constant inflow
+                                 0,
+                                 0])
+
+        # All reflective to begin with (still water) 
+        domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
+
+
+        #--------------------------------------------------------------
+        # Test initial inundation height
+        #--------------------------------------------------------------
+
+        indices = domain.get_wet_elements()
+        z = domain.get_quantity('elevation').\
+            get_values(location='centroids', indices=indices)
+        assert alltrue(z<initial_runup_height)
+
+        q_ref = domain.get_maximum_inundation_elevation()
+        assert allclose(q_ref, initial_runup_height, rtol = 1.0/N) # First order accuracy 
+
+        
+        #--------------------------------------------------------------
+        # Let triangles adjust
+        #--------------------------------------------------------------
+        for t in domain.evolve(yieldstep = 0.1, finaltime = 1.0):
+            pass
+
+
+        #--------------------------------------------------------------
+        # Test inundation height again
+        #--------------------------------------------------------------
+       
+        q_ref = domain.get_maximum_inundation_elevation()
+        q = get_maximum_inundation_elevation('runup_test.sww')
+        msg = 'We got %f, should have been %f' %(q, q_ref)
+        assert allclose(q, q_ref, rtol=1.0/N), msg
+
+
+        q = get_maximum_inundation_elevation('runup_test.sww')
+        msg = 'We got %f, should have been %f' %(q, initial_runup_height)
+        assert allclose(q, initial_runup_height, rtol = 1.0/N), msg 
+
+        
+
+
+        #--------------------------------------------------------------
+        # Update boundary to allow inflow
+        #--------------------------------------------------------------
+        domain.set_boundary({'right': Bd})
+
+        
+        #--------------------------------------------------------------
+        # Evolve system through time
+        #--------------------------------------------------------------
+        q_max = None
+        for t in domain.evolve(yieldstep = 0.1, finaltime = 3.0):
+            q = domain.get_maximum_inundation_elevation()
+            if q > q_max: q_max = q
+
+    
+        #--------------------------------------------------------------
+        # Test inundation height again
+        #--------------------------------------------------------------
+
+        indices = domain.get_wet_elements()
+        z = domain.get_quantity('elevation').\
+            get_values(location='centroids', indices=indices)
+
+        assert alltrue(z<final_runup_height)
+
+        q = domain.get_maximum_inundation_elevation()
+        assert allclose(q, final_runup_height, rtol = 1.0/N) # First order accuracy
+
+        q = get_maximum_inundation_elevation('runup_test.sww', timesteps=31)
+        msg = 'We got %f, should have been %f' %(q, final_runup_height)
+        assert allclose(q, final_runup_height, rtol=1.0/N), msg        
+
+
+        q = get_maximum_inundation_elevation('runup_test.sww')
+        msg = 'We got %f, should have been %f' %(q, q_max)
+        assert allclose(q, q_max, rtol=1.0/N), msg
+
+        q = get_maximum_inundation_elevation('runup_test.sww', timesteps=range(32))
+        msg = 'We got %f, should have been %f' %(q, q_max)
+        assert allclose(q, q_max, rtol=1.0/N), msg                
+
+
 
 
@@ -4426 +4558 @@
 if __name__ == "__main__":
     suite = unittest.makeSuite(Test_Shallow_Water,'test')
-    #suite = unittest.makeSuite(Test_Shallow_Water,'test_get_maximum_inundation_3')
+    #suite = unittest.makeSuite(Test_Shallow_Water,'test_get_maximum_inundation_from_sww')
     #suite = unittest.makeSuite(Test_Shallow_Water,'test_temp')    
     

anuga_core/source/anuga/utilities/numerical_tools.py

@@ -246 +246 @@
                      If not, an attempt is made to convert it to a Numeric
                      array
+        A: Scalar.   Return 0-dimensional array of length 1, containing that value
+        A: String.   Array of ASCII values
         typecode: Numeric type. If specified, use this in the conversion.
                                 If not, let Numeric decide

anuga_core/source/anuga/utilities/test_numerical_tools.py

@@ -77 +77 @@
     def test_ensure_numeric(self):
         from numerical_tools import ensure_numeric
-        from Numeric import ArrayType, Float, array
+        from Numeric import ArrayType, Float, Int, array
 
         A = [1,2,3,4]
@@ -122 +122 @@
         assert A is not B   #Not the same object
 
+        # Check scalars
+        A = 1
+        B = ensure_numeric(A, Float)
+        #print A, B[0], len(B), type(B) 
+        #print B.shape
+        assert A == B
+
+        B = ensure_numeric(A, Int)        
+        #print A, B
+        #print B.shape
+        assert A == B
+
+        # Error situation
+
+        B = ensure_numeric('hello', Int)                
+        assert allclose(B, [104, 101, 108, 108, 111])
 
     def test_gradient(self):