Changeset 4701

Timestamp:

Sep 4, 2007, 5:36:01 PM (17 years ago)

Author:

ole

Message:

Work towards ticket:192

Location:

anuga_core/source/anuga

Files:

• abstract_2d_finite_volumes/domain.py (modified) (8 diffs)
• shallow_water/shallow_water_domain.py (modified) (7 diffs)
• shallow_water/test_data_manager.py (modified) (2 diffs)

anuga_core/source/anuga/abstract_2d_finite_volumes/domain.py

@@ -170 +170 @@
                         
 
-        #Defaults
+        # Defaults
         from anuga.config import max_smallsteps, beta_w, beta_h, epsilon
         from anuga.config import CFL
@@ -180 +180 @@
         
 
-        #FIXME: Maybe have separate orders for h-limiter and w-limiter?
-        #Or maybe get rid of order altogether and use beta_w and beta_h
+        # FIXME: Maybe have separate orders for h-limiter and w-limiter?
+        # Or maybe get rid of order altogether and use beta_w and beta_h
         self.set_default_order(1)
         #self.default_order = 1
@@ -195 +195 @@
         
 
-        #Model time
+        # Model time
         self.time = 0.0
         self.finaltime = None
@@ -201 +201 @@
         self.starttime = 0 #Physical starttime if any (0 is 1 Jan 1970 00:00:00)
 
-        ######OBSOLETE
-        #Origin in UTM coordinates
-        #FIXME: This should be set if read by a msh file
-        #self.zone = zone
-        #self.xllcorner = xllcorner
-        #self.yllcorner = yllcorner
-
-
-        #Checkpointing and storage
+        # Monitoring
+        self.quantities_to_be_monitored = None
+        self.monitor_polygon = None
+        self.monitor_time_interval = None                
+
+
+        # Checkpointing and storage
         from anuga.config import default_datadir
         self.datadir = default_datadir
@@ -215 +213 @@
         self.checkpoint = False
 
-        #MH310505 To avoid calculating the flux across each edge twice, keep an integer (boolean) array,
-        #to be used during the flux calculation
+        # MH310505 To avoid calculating the flux across each edge twice, keep an integer (boolean) array,
+        # to be used during the flux calculation
         N = len(self) #number_of_triangles
         self.already_computed_flux = zeros((N, 3), Int)
@@ -301 +299 @@
 
         """
+        
+        # FIXME: Could we name this a bit more intuitively
+        # E.g. set_quantities_from_dictionary
         for key in quantity_dict.keys():
             self.set_quantity(key, quantity_dict[key], location='vertices')
@@ -544 +545 @@
 
 
+
+
+    def set_quantities_to_be_monitored(self, q,
+                                       polygon=None,
+                                       time_interval=None):
+        """Specify which quantities will be monitored for extrema.
+
+        q must be either:
+          - the name of a quantity
+          - a list of quantity names
+          - None
+
+        In the two first cases, the named quantities will be monitored at
+        each internal timestep
+        
+        If q is None, monitoring will be switched off altogether.
+
+        polygon (if specified) will restrict monitoring to triangles inside polygon.
+        If omitted all triangles will be included.
+
+        time_interval (if specified) will restrict monitoring to time steps in
+        that interval. If omitted all timesteps will be included.
+
+        FIXME: Derived quantities such as 'stage-elevation' will appear later
+        """
+
+        # FIXME (Ole): This is under construction. See ticket:192
+
+        if q is None:
+            self.quantities_to_be_monitored = None
+            self.monitor_polygon = None
+            self.monitor_time_interval = None                    
+            return
+
+        if isinstance(q, basestring):
+            q = [q] # Turn argument into a list
+
+        # Check correcness
+        for quantity_name in q:
+            msg = 'Quantity %s is not a valid conserved quantity'\
+                  %quantity_name
+            
+            assert quantity_name in self.conserved_quantities, msg
+
+        if polygon is not None:
+            # FIXME Check input
+            pass
+
+        if time_interval is not None:
+            # FIXME Check input
+            pass        
+
+
+        self.quantities_to_be_monitored = q
+        self.monitor_polygon = polygon
+        self.monitor_time_interval = time_interval        
+        
+
+
     #MISC
     def check_integrity(self):
@@ -742 +802 @@
 
         return msg
+
+
+
+    def quantity_statistics(self):
+        """Return string with statistics about quantities for printing or logging
+
+        Quantities reported are specified through method
+
+           set_quantities_to_be_monitored
+           
+        """
+
+        pass
+
+
 
 

anuga_core/source/anuga/shallow_water/shallow_water_domain.py

@@ -173 +173 @@
         self.forcing_terms.append(gravity)
 
-        #Stored output
+        # Stored output
         self.store = True
         self.format = 'sww'
@@ -179 +179 @@
         self.minimum_storable_height = minimum_storable_height
         self.quantities_to_be_stored = ['stage','xmomentum','ymomentum']
+        
                 
 
@@ -234 +235 @@
 
         #FIXME (Ole): rename H0 to minimum_allowed_height_in_flux_computation
-        #rename tight_slope_limiters to tight_slope_limiters.
-        #Maybe use histogram to identify isolated extreme speeds and deal with them adaptively
-        #similarly to how we used to use 1 order steps to recover.
+
+        #FIXME (Ole): Maybe use histogram to identify isolated extreme speeds
+        #and deal with them adaptively similarly to how we used to use 1 order
+        #steps to recover.
         self.minimum_allowed_height = minimum_allowed_height
         self.H0 = minimum_allowed_height   
@@ -262 +264 @@
         self.points_file_block_line_size = points_file_block_line_size
         
+        
     def set_quantities_to_be_stored(self, q):
         """Specify which quantities will be stored in the sww file.
@@ -286 +289 @@
             q = [q] # Turn argument into a list
 
-        #Check correcness
+        # Check correcness
         for quantity_name in q:
             msg = 'Quantity %s is not a valid conserved quantity'\
@@ -930 +933 @@
         #MH090605 if second order,
         #perform the extrapolation and limiting on
-        #all of the conserved quantitie
+        #all of the conserved quantities
 
         if (domain._order_ == 1):
@@ -944 +947 @@
         for name in domain.conserved_quantities:
             Q = domain.quantities[name]
+
             if domain._order_ == 1:
                 Q.extrapolate_first_order()
             elif domain._order_ == 2:
+
+                # Experiment
+                #if name == 'stage':
+                #    #print name, 'second'
+                #    Q.extrapolate_second_order()
+                #    Q.limit()
+                #else:
+                #    #print name, 'first'                
+                #    Q.extrapolate_first_order()
+                #    #Q.extrapolate_second_order()
+                #    #Q.limit()                
+                
                 Q.extrapolate_second_order()
                 Q.limit()

anuga_core/source/anuga/shallow_water/test_data_manager.py

@@ -59 +59 @@
         ######################
         #Initial condition - with jumps
-
-
         bed = domain.quantities['elevation'].vertex_values
         stage = zeros(bed.shape, Float)
@@ -233 +231 @@
         self.domain.format = 'sww'
         self.domain.smooth = True
-
-        sww = get_dataobject(self.domain)
-        sww.store_connectivity()
+        sww = get_dataobject(self.domain)        
+
         for t in self.domain.evolve(yieldstep = 1, finaltime = 1):
             pass
             
-        #Get NetCDF
-        fid = NetCDFFile(sww.filename, 'r')  #Open existing file for append
-
-        from anuga.shallow_water.shallow_water_domain import Domain 
-        sww_quantities = Domain.conserved_quantities
+        # Get NetCDF
+        fid = NetCDFFile(sww.filename, 'r') # Open existing file for append
+
         # Get the variables
         range = fid.variables['stage_range'][:]
         assert allclose(range,[-0.93519, 0.15]) or\
                allclose(range,[-0.9352743, 0.15]) # Old slope limiters
+        
         range = fid.variables['xmomentum_range'][:]
-        
-        #assert allclose(range,[0,0.46950444])
         assert allclose(range,[0,0.4695096]) or\
-               allclose(range,[0,0.47790655]) # Old slope limiters             
+               allclose(range,[0,0.47790655]) # Old slope limiters
+        
         range = fid.variables['ymomentum_range'][:]
-        
         #assert allclose(range,[0,0.02174380])
         assert allclose(range,[0,0.02174439]) or\
-               allclose(range,[0,0.02283983]) # Old slope limiters                     
+               allclose(range,[0,0.02283983]) # Old slope limiters 
+        
+        fid.close()
+        os.remove(sww.filename)
+
+    def NOtest_sww_extrema(self):
+        """Test that extrema of quantities can be retrieved at every vertex
+        Extrema are updated at every *internal* timestep
+        """
+
+        domain = self.domain
+        
+        domain.set_name('datatest' + str(id(self)))
+        domain.format = 'sww'
+        domain.smooth = True
+
+        assert domain.quantities_to_be_monitored is None
+        assert domain.monitor_polygon is None
+        assert domain.monitor_time_interval is None        
+        
+        domain.set_quantities_to_be_monitored(['stage', 'ymomentum'])
+        
+        assert len(domain.quantities_to_be_monitored) == 2
+        assert domain.quantities_to_be_monitored[0] == 'stage'
+        assert domain.quantities_to_be_monitored[1] == 'ymomentum'        
+        assert domain.monitor_polygon is None
+        assert domain.monitor_time_interval is None        
+        
+        sww = get_dataobject(domain)
+
+        for t in domain.evolve(yieldstep = 1, finaltime = 1):
+            print domain.timestepping_statistics()
+            print domain.quantity_statistics()
+
+            
+        # Get NetCDF
+        fid = NetCDFFile(sww.filename, 'r') # Open existing file for append
+
+        # Get the variables
+        extrema = fid.variables['stage_extrema'][:]
+        assert allclose(range, [])
+
+        extrema = fid.variables['xmomentum_extrema'][:]
+        assert allclose(range,[])
+        
+        extrema = fid.variables['ymomentum_extrema'][:]
+        assert allclose(range,[])
+
         
         fid.close()