Changeset 9018

Timestamp:

Nov 8, 2013, 8:45:50 PM (11 years ago)

Author:

steve

Message:

Added in the possibility to store centroid values into sww files

Location:

trunk/anuga_core/source/anuga

Files:

• file/sww.py (modified) (18 diffs)
• file/test_read_sww.py (modified) (2 diffs)
• file/test_sww.py (modified) (1 diff)
• shallow_water/shallow_water_domain.py (modified) (4 diffs)

trunk/anuga_core/source/anuga/file/sww.py

@@ -71 +71 @@
         self.recursion = recursion
         self.mode = mode
+        
         if hasattr(domain, 'max_size'):
             self.max_size = domain.max_size # File size max is 2Gig
         else:
             self.max_size = max_size
+            
+        if hasattr(domain, 'store_centroids'):
+            self.store_centroids = domain.store_centroids
+        else:
+            self.store_centroids = False
+            
         if hasattr(domain, 'minimum_storable_height'):
             self.minimum_storable_height = domain.minimum_storable_height
@@ -105 +112 @@
                           'suitable for plotting'
                           
-            self.writer = Write_sww(static_quantities, dynamic_quantities)
+            self.writer = Write_sww(static_quantities, dynamic_quantities, store_centroids=self.store_centroids)
             self.writer.store_header(fid,
                                      domain.starttime,
@@ -194 +201 @@
         # Get names of static quantities
         static_quantities = {}
+        static_quantities_centroid = {}
         for name in self.writer.static_quantities:
             Q = domain.quantities[name]
@@ -199 +207 @@
                                        precision=self.precision)
             static_quantities[name] = A
+            
+            if self.store_centroids:
+                static_quantities_centroid[name] = Q.centroid_values
         
         # Store static quantities        
         self.writer.store_static_quantities(fid, **static_quantities)
-                                            
+        
+        if  self.store_centroids:
+            self.writer.store_static_quantities_centroid(fid, **static_quantities_centroid) 
+                                        
         fid.close()
 
@@ -319 +333 @@
             # Now store dynamic quantities
             dynamic_quantities = {}
+            dynamic_quantities_centroid = {}
+            
             for name in self.writer.dynamic_quantities:
                 #netcdf_array = fid.variables[name]
@@ -341 +357 @@
                 dynamic_quantities[name] = A
                 
+                if self.store_centroids:
+                    dynamic_quantities_centroid[name] = Q.centroid_values
+                
                                         
             # Store dynamic quantities
-            self.writer.store_quantities(fid,
+            slice_index = self.writer.store_quantities(fid,
                                          time=self.domain.time,
                                          sww_precision=self.precision,
                                          **dynamic_quantities)
+            
+            # Store dynamic quantities
+            if self.store_centroids:
+                self.writer.store_quantities_centroid(fid,
+                                                      slice_index= slice_index,
+                                                      sww_precision=self.precision,
+                                                      **dynamic_quantities_centroid)            
 
 
@@ -430 +456 @@
         
         for q in filter(lambda n:n != 'x' and n != 'y' and n != 'time' and n != 'volumes' and \
-                        '_range' not in n, \
+                        '_range' not in n and '_c' not in n , \
                         fin.variables.keys()):
+            #print q
             if len(fin.variables[q].shape) == 1: # Not a time-varying quantity
                 self.quantities[q] = num.ravel(num.array(fin.variables[q], num.float)).reshape(M,3)
@@ -466 +493 @@
     """
 
-    def __init__(self, static_quantities, dynamic_quantities):
+    def __init__(self, static_quantities, dynamic_quantities, store_centroids=False):
         """Initialise Write_sww with two list af quantity names: 
         
@@ -479 +506 @@
         self.static_quantities = static_quantities   
         self.dynamic_quantities = dynamic_quantities
+        self.store_centroids = store_centroids
 
 
@@ -577 +605 @@
             outfile.createVariable(q + Write_sww.RANGE, sww_precision,
                                    ('numbers_in_range',))
+            
+            if self.store_centroids:
+                outfile.createVariable(q + '_c', sww_precision,
+                                   ('number_of_volumes',))
                                    
             # Initialise ranges with small and large sentinels.
@@ -702 +734 @@
 
 
+    def write_dynamic_quantities(self, outfile, quantities,
+                    times, precis = netcdf_float32, verbose = False):   
+        """
+            Write out given quantities to file.
+        """
+        
+
+        for q in quantities:
+            outfile.createVariable(q, precis, ('number_of_timesteps',
+                                                      'number_of_points'))
+            outfile.createVariable(q + Write_sts.RANGE, precis,
+                                   ('numbers_in_range',))
+            
+            if self.store_centroids:
+                outfile.createVariable(q + '_c', precis,
+                                   ('number_of_timesteps','number_of_volumes'))
+
+            # Initialise ranges with small and large sentinels.
+            # If this was in pure Python we could have used None sensibly
+            outfile.variables[q+Write_sts.RANGE][0] = max_float  # Min
+            outfile.variables[q+Write_sts.RANGE][1] = -max_float # Max
+
+        # Doing sts_precision instead of Float gives cast errors.
+        outfile.createVariable('time', netcdf_float, ('number_of_timesteps',))
+
+        if isinstance(times, (list, num.ndarray)):
+            outfile.variables['time'][:] = times    # Store time relative
+
+        if verbose:
+            log.critical('------------------------------------------------')
+            log.critical('Statistics:')
+            log.critical('    t in [%f, %f], len(t) == %d'
+                         % (num.min(times), num.max(times), len(times.flat)))
 
     def store_parallel_data(self,
@@ -795 +860 @@
 
                     
-        
+
+    def store_static_quantities_centroid(self, 
+                                            outfile, 
+                                            sww_precision=num.float32,
+                                            verbose=False, 
+                                            **quant):
+        """
+        Write the static centroid quantity info.
+
+        **quant is extra keyword arguments passed in. These must be
+          the numpy arrays to be stored in the sww file at each timestep.
+
+        The argument sww_precision allows for storing as either 
+        * single precision (default): num.float32
+        * double precision: num.float64 or num.float 
+
+        Precondition:
+            store_triangulation and
+            store_header have been called.
+        """
+
+        # The dictionary quant must contain numpy arrays for each name.
+        # These will typically be quantities from Domain such as friction 
+        #
+        # Arrays not listed in static_quantitiues will be ignored, silently.
+        #
+        # This method will also write the ranges for each quantity, 
+        # e.g. stage_range, xmomentum_range and ymomentum_range
+        for q in self.static_quantities:
+            if not quant.has_key(q):
+                msg = 'Values for quantity %s was not specified in ' % q
+                msg += 'store_quantities so they cannot be stored.'
+                raise NewQuantity, msg
+            else:
+                q_values = ensure_numeric(quant[q])
+                
+                x = q_values.astype(sww_precision)
+                outfile.variables[q+'_c'][:] = x
+        
+
+                    
+       
 
     def store_quantities(self, 
@@ -862 +968 @@
                 if q_values_max > q_range[1]:
                     outfile.variables[q + Write_sww.RANGE][1] = q_values_max
-
+                    
+        return slice_index
+                    
+                    
+
+    def store_quantities_centroid(self, 
+                                  outfile, 
+                                  sww_precision=num.float32,
+                                  slice_index=None,
+                                  verbose=False, 
+                                  **quant):
+        """
+        Write the quantity centroid info at each timestep.
+
+        **quant is extra keyword arguments passed in. These must be
+          the numpy arrays to be stored in the sww file at each timestep.
+
+        if the time array is already been built, use the slice_index
+        to specify the index.
+
+        Otherwise, use time to increase the time dimension
+
+        Maybe make this general, but the viewer assumes these quantities,
+        so maybe we don't want it general - unless the viewer is general
+        
+        The argument sww_precision allows for storing as either 
+        * single precision (default): num.float32
+        * double precision: num.float64 or num.float 
+
+        Precondition:
+            store_triangulation and
+            store_header have been called.
+        """
+
+        assert slice_index is not None, 'slice_index should be set in store_quantities'
+          
+
+        # Write the named dynamic quantities
+        # The dictionary quant must contain numpy arrays for each name.
+        # These will typically be the conserved quantities from Domain 
+        # (Typically stage,  xmomentum, ymomentum).
+        #
+        # Arrays not listed in dynamic_quantitiues will be ignored, silently.
+        #
+        # This method will also write the ranges for each quantity, 
+        # e.g. stage_range, xmomentum_range and ymomentum_range
+        for q in self.dynamic_quantities:
+            if not quant.has_key(q):
+                msg = 'Values for quantity %s was not specified in ' % q
+                msg += 'store_quantities so they cannot be stored.'
+                raise NewQuantity, msg
+            else:
+                q_values = ensure_numeric(quant[q])
+                
+                q_retyped = q_values.astype(sww_precision)
+                outfile.variables[q+'_c'][slice_index] = q_retyped
 
 
@@ -947 +1108 @@
     # or              concatenate((x[:,num.newaxis], x[:,num.newaxis]), axis=1)
 
-    conserved_quantities = []
+    dynamic_quantities = []
     interpolated_quantities = {}
-    other_quantities = []
+    static_quantities = []
 
     # get geo_reference
@@ -962 +1123 @@
         dimensions = fid.variables[quantity].dimensions
         if 'number_of_timesteps' in dimensions:
-            conserved_quantities.append(quantity)
+            dynamic_quantities.append(quantity)
             interpolated_quantities[quantity] = \
                   interpolated_quantity(fid.variables[quantity][:], time_interp)
         else:
-            other_quantities.append(quantity)
-
-    other_quantities.remove('x')
-    other_quantities.remove('y')
-    #other_quantities.remove('z')
-    other_quantities.remove('volumes')
+            static_quantities.append(quantity)
+
+    #print static_quantities
+    #print dynamic_quantities
+    
     try:
-        other_quantities.remove('stage_range')
-        other_quantities.remove('xmomentum_range')
-        other_quantities.remove('ymomentum_range')
-        other_quantities.remove('elevation_range')
+        dynamic_quantities.remove('stage_c')
+        dynamic_quantities.remove('xmomentum_c')
+        dynamic_quantities.remove('ymomentum_c')
+        dynamic_quantities.remove('elevation_c')
     except:
         pass
-
-    conserved_quantities.remove('time')
+    
+    try:
+        static_quantities.remove('elevation_c')
+    except:
+        pass
+        
+    
+    static_quantities.remove('x')
+    static_quantities.remove('y')
+    #other_quantities.remove('z')
+    static_quantities.remove('volumes')
+    try:
+        static_quantities.remove('stage_range')
+        static_quantities.remove('xmomentum_range')
+        static_quantities.remove('ymomentum_range')
+        static_quantities.remove('elevation_range')
+    except:
+        pass
+
+    dynamic_quantities.remove('time')
 
     if verbose: log.critical('    building domain')
@@ -999 +1177 @@
         unique = True
     if unique:
-        coordinates, volumes, boundary = weed(coordinates, volumes,boundary)
+        coordinates, volumes, boundary = weed(coordinates, volumes, boundary)
 
       
@@ -1016 +1194 @@
     domain.geo_reference = geo_reference
 
-    for quantity in other_quantities:
+    for quantity in static_quantities:
         try:
             NaN = fid.variables[quantity].missing_value
@@ -1041 +1219 @@
         domain.set_quantity(quantity, X)
     #
-    for quantity in conserved_quantities:
+    for quantity in dynamic_quantities:
         try:
             NaN = fid.variables[quantity].missing_value

trunk/anuga_core/source/anuga/file/test_read_sww.py

@@ -112 +112 @@
         assert num.allclose(sww_file.y, domain.get_vertex_coordinates()[:,1])
 
-
+        
         assert num.allclose(sww_file.time, [0.0, 1.0, 2.0, 3.0, 4.0])
         
@@ -179 +179 @@
         os.remove(source)
         
+    def test_read_sww_with_centroids(self):
+        """
+        Save to an sww file and then read back the info.
+        Here we store the info "uniquely"
+        """
+
+        #---------------------------------------------------------------------
+        # Import necessary modules
+        #---------------------------------------------------------------------
+        from anuga.abstract_2d_finite_volumes.mesh_factory import \
+            rectangular_cross
+        from anuga.shallow_water.shallow_water_domain import Domain
+        from boundaries import Reflective_boundary
+        from anuga.abstract_2d_finite_volumes.generic_boundary_conditions\
+                            import Dirichlet_boundary, Time_boundary
+
+        #---------------------------------------------------------------------
+        # Setup computational domain
+        #---------------------------------------------------------------------
+        length = 8.
+        width = 4.
+        dx = dy = 2   # Resolution: Length of subdivisions on both axes
+        
+        inc = 0.05 # Elevation increment
+
+        points, vertices, boundary = rectangular_cross(int(length/dx), 
+                                                       int(width/dy),
+                                                       len1=length, 
+                                                       len2=width)
+        domain = Domain(points, vertices, boundary)
+        domain.set_name('read_sww_test_c'+str(domain.processor))  # Output name
+        domain.set_quantities_to_be_stored({'elevation': 2,
+                                            'stage': 2,
+                                            'xmomentum': 2,
+                                            'ymomentum': 2,
+                                            'friction': 1})
+
+        domain.set_store_vertices_uniquely(True)
+        domain.set_store_centroids(True)
+        
+        #---------------------------------------------------------------------
+        # Setup initial conditions
+        #---------------------------------------------------------------------
+        domain.set_quantity('elevation', 0.0)    # Flat bed initially
+        domain.set_quantity('friction', 0.01)    # Constant friction
+        domain.set_quantity('stage', 0.0)        # Dry initial condition
+
+        #------------------------------------------------------------------
+        # Setup boundary conditions
+        #------------------------------------------------------------------
+        Bi = Dirichlet_boundary([0.4, 0, 0])          # Inflow
+        Br = Reflective_boundary(domain)              # Solid reflective wall
+        Bo = Dirichlet_boundary([-5, 0, 0])           # Outflow
+
+        domain.set_boundary({'left': Bi, 'right': Bo, 'top': Br, 'bottom': Br})
+
+        #-------------------------------------------------------------------
+        # Evolve system through time
+        #-------------------------------------------------------------------
+
+        for t in domain.evolve(yieldstep=1, finaltime=4.0):
+            pass
+            
+            
+        # Check that quantities have been stored correctly    
+        source = domain.get_name() + '.sww'
+
+
+        #x = fid.variables['x'][:]
+        #y = fid.variables['y'][:]
+        #stage = fid.variables['stage'][:]
+        #elevation = fid.variables['elevation'][:]
+        #fid.close()
+                   
+        #assert len(stage.shape) == 2
+        #assert len(elevation.shape) == 2        
+        
+        #M, N = stage.shape
+                
+        sww_file = sww.Read_sww(source)
+
+        #print 'last frame number',sww_file.get_last_frame_number()
+
+        assert num.allclose(sww_file.x, domain.get_vertex_coordinates()[:,0])
+        assert num.allclose(sww_file.y, domain.get_vertex_coordinates()[:,1])
+
+        
+        assert num.allclose(sww_file.time, [0.0, 1.0, 2.0, 3.0, 4.0])
+        
+        M = domain.get_number_of_triangles()
+        
+        assert num.allclose(num.reshape(num.arange(3*M), (M,3)), sww_file.vertices)
+
+        last_frame_number = sww_file.get_last_frame_number() 
+        assert last_frame_number == 4
+
+        assert num.allclose(sww_file.get_bounds(), [0.0, length, 0.0, width])
+
+        assert 'stage'     in sww_file.quantities.keys()
+        assert 'friction'  in sww_file.quantities.keys()
+        assert 'elevation' in sww_file.quantities.keys()
+        assert 'xmomentum' in sww_file.quantities.keys()
+        assert 'ymomentum' in sww_file.quantities.keys()
+
+
+        for qname, q in sww_file.read_quantities(last_frame_number).items():
+            assert num.allclose(domain.get_quantity(qname).get_values(), q)
+
+        #-----------------------------------------
+        # Start the evolution off again at frame 3
+        #-----------------------------------------
+        sww_file.read_quantities(last_frame_number-1)
+
+        points, vertices, boundary = rectangular_cross(int(length/dx), 
+                                                       int(width/dy),
+                                                       len1=length, 
+                                                       len2=width)
+        new_domain = Domain(points, vertices, boundary)
+        new_domain.set_quantities_to_be_stored(None)
+
+        new_domain.set_store_vertices_uniquely(True)
+
+        for qname, q in sww_file.read_quantities(last_frame_number-1).items():
+            new_domain.set_quantity(qname, q)    
+
+        #------------------------------------------------------------------
+        # Setup boundary conditions
+        #------------------------------------------------------------------
+        Bi = Dirichlet_boundary([0.4, 0, 0])          # Inflow
+        Br = Reflective_boundary(new_domain)          # Solid reflective wall
+        Bo = Dirichlet_boundary([-5, 0, 0])           # Outflow
+
+        new_domain.set_boundary({'left': Bi, 'right': Bo, 'top': Br, 'bottom': Br})
+
+        #-------------------------------------------------------------------
+        # Evolve system through time
+        #-------------------------------------------------------------------
+
+        for t in new_domain.evolve(yieldstep=1.0, finaltime=1.0):
+             pass
+
+        # Compare  new_domain and domain quantities 
+        for quantity in domain.get_quantity_names():
+            dv = domain.get_quantity(quantity).get_values()
+            ndv = new_domain.get_quantity(quantity).get_values()
+
+            #print dv-ndv
+
+            assert num.allclose( dv, ndv, rtol=5.e-2, atol=5.e-2)
+
+        # Clean up
+        #os.remove(source)
+        
 
 if __name__ == "__main__":

trunk/anuga_core/source/anuga/file/test_sww.py

@@ -460 +460 @@
 if __name__ == "__main__":
     suite = unittest.makeSuite(Test_sww, 'test')
-    runner = unittest.TextTestRunner(verbosity=1)
+    runner = unittest.TextTestRunner(verbosity=2)
     runner.run(suite)

trunk/anuga_core/source/anuga/shallow_water/shallow_water_domain.py

@@ -208 +208 @@
         #-------------------------------
         self.set_store(True)
+        self.set_store_centroids(False)
         self.set_store_vertices_uniquely(False)
         self.quantities_to_be_stored = {'elevation': 1, 
@@ -213 +214 @@
                                         'xmomentum': 2, 
                                         'ymomentum': 2}
-
-        
-
 
         #-------------------------------
@@ -231 +229 @@
     def set_defaults(self):
         """Set the default values in this routine. That way we can inherit class
-        and just over redefine the defaults for the new class
+        and just redefine the defaults for the new class
         """
 
@@ -451 +449 @@
         return self.store
 
+
+    def set_store_centroids(self, flag=True):
+        """Set whether centroid data is saved to sww file.
+        """
+
+        self.store_centroids = flag
+        
+    def get_store_centroids(self):
+        """Get whether data saved to sww file.
+        """
+        
+        return self.store_centroids   
+    
+        
         
     def set_sloped_mannings_function(self, flag=True):