Changeset 4704

Timestamp:

Sep 5, 2007, 4:42:58 PM (18 years ago)

Author:

ole

Message:

Implemented ticket:192

Location:

anuga_core/source/anuga

Files:

• abstract_2d_finite_volumes/domain.py (modified) (15 diffs)
• abstract_2d_finite_volumes/quantity.py (modified) (5 diffs)
• abstract_2d_finite_volumes/test_domain.py (modified) (2 diffs)
• abstract_2d_finite_volumes/test_quantity.py (modified) (5 diffs)
• abstract_2d_finite_volumes/util.py (modified) (1 diff)
• shallow_water/data_manager.py (modified) (25 diffs)
• shallow_water/test_data_manager.py (modified) (11 diffs)
• shallow_water/test_shallow_water_domain.py (modified) (2 diffs)

anuga_core/source/anuga/abstract_2d_finite_volumes/domain.py

@@ -26 +26 @@
 from anuga.abstract_2d_finite_volumes.region\
      import Set_region as region_set_region
+
+from anuga.utilities.polygon import inside_polygon
+from anuga.abstract_2d_finite_volumes.util import get_textual_float
 
 import types
@@ -205 +208 @@
         self.monitor_polygon = None
         self.monitor_time_interval = None                
-
+        self.monitor_indices = None
 
         # Checkpointing and storage
@@ -579 +582 @@
             self.quantities_to_be_monitored = None
             self.monitor_polygon = None
-            self.monitor_time_interval = None                    
+            self.monitor_time_interval = None
+            self.monitor_indices = None            
             return
 
@@ -596 +600 @@
                 apply_expression_to_dictionary(quantity_name, self.quantities)
 
-            # Initialise extrema
-            self.quantities_to_be_monitored[quantity_name] = [None, None]
+            # Initialise extrema information
+            info_block = {'min': None,          # Min value
+                          'max': None,          # Max value
+                          'min_location': None, # Argmin (x, y)
+                          'max_location': None, # Argmax (x, y)
+                          'min_time': None,     # Argmin (t) 
+                          'max_time': None}     # Argmax (t)
             
+            self.quantities_to_be_monitored[quantity_name] = info_block
+
+            
 
         if polygon is not None:
-            # FIXME Check input
-            pass
+            # Check input
+            if isinstance(polygon, basestring):
+
+                # Check if multiple quantities were accidentally
+                # given as separate argument rather than a list.
+                msg = 'Multiple quantities must be specified in a list. '
+                msg += 'Not as multiple arguments. '
+                msg += 'I got "%s" as a second argument' %polygon
+                
+                if polygon in self.quantities:
+                    raise Exception, msg
+                
+                try:
+                    apply_expression_to_dictionary(polygon, self.quantities)
+                except:
+                    # At least polygon wasn't an expression involving quantitites
+                    pass
+                else:
+                    raise Exception, msg
+
+                # In any case, we don't allow polygon to be a string
+                msg = 'argument "polygon" must not be a string: '
+                msg += 'I got polygon=\'%s\' ' %polygon
+                raise Exception, msg
+
+
+            # Get indices for centroids that are inside polygon
+            points = self.get_centroid_coordinates(absolute=True)
+            self.monitor_indices = inside_polygon(points, polygon)
+            
 
         if time_interval is not None:
-            # FIXME Check input
-            pass        
-
+            assert len(time_interval) == 2
 
         
@@ -758 +796 @@
         """
 
-        #Input checks
+        # Input checks
         import types, string
 
@@ -780 +818 @@
         assert type(tags) == types.ListType, msg
 
-        #Determine width of longest quantity name (for cosmetic purposes)
+        # Determine width of longest quantity name (for cosmetic purposes)
         maxwidth = 0
         for name in quantities:
@@ -787 +825 @@
                 maxwidth = w
 
-        #Output stats
+        # Output stats
         msg = 'Boundary values at time %.4f:\n' %self.time
         for tag in tags:
@@ -795 +833 @@
                 q = self.quantities[name]
 
-                #Find range of boundary values for tag and q
+                # Find range of boundary values for tag and q
                 maxval = minval = None
                 for i, ((vol_id, edge_id), B) in\
@@ -815 +853 @@
 
 
-
-    def quantity_statistics(self):
+    def update_extrema(self):
+        """Update extrema if requested by set_quantities_to_be_monitored.
+        This data is used for reporting e.g. by running
+        print domain.quantity_statistics()
+        and may also stored in output files (see data_manager in shallow_water)
+        """
+
+        if self.quantities_to_be_monitored is None:
+            return
+
+        # Observe time interval restriction if any
+        if self.monitor_time_interval is not None and\
+               (self.time < self.monitor_time_interval[0] or\
+               self.time > self.monitor_time_interval[1]):
+            return
+            
+        
+        # Update extrema for each specified quantity subject to
+        # polygon restriction (via monitor_indices).
+        for quantity_name in self.quantities_to_be_monitored:
+
+            if quantity_name in self.quantities:
+                Q = self.get_quantity(quantity_name)
+            else:
+                Q = self.create_quantity_from_expression(quantity_name)
+
+            info_block = self.quantities_to_be_monitored[quantity_name]
+
+            # Update maximum (n > None is always True)
+            maxval = Q.get_maximum_value(self.monitor_indices)
+            if maxval > info_block['max']:
+                info_block['max'] = maxval
+                maxloc = Q.get_maximum_location()
+                info_block['max_location'] = maxloc
+                info_block['max_time'] = self.time
+
+
+            # Update minimum (n < None is always False)
+            minval = Q.get_minimum_value(self.monitor_indices)
+            if info_block['min'] is None or\
+                   minval < info_block['min']:
+                info_block['min'] = minval                
+                minloc = Q.get_minimum_location()
+                info_block['min_location'] = minloc
+                info_block['min_time'] = self.time                
+        
+
+
+    def quantity_statistics(self, precision = '%.4f'):
         """Return string with statistics about quantities for printing or logging
 
@@ -825 +910 @@
         """
 
-        pass
-
-
+        maxlen = 128 # Max length of polygon string representation
+
+        # Output stats
+        msg = 'Monitored quantities at time %.4f:\n' %self.time
+        if self.monitor_polygon is not None:
+            p_str = str(self.monitor_polygon)
+            msg += '- Restricted by polygon: %s' %p_str[:maxlen]
+            if len(p_str) >= maxlen:
+                msg += '...\n'
+            else:
+                msg += '\n'
+
+
+        if self.monitor_time_interval is not None:
+            msg += '- Restricted by time interval: %s\n' %str(self.monitor_time_interval)
+            time_interval_start = self.monitor_time_interval[0]
+        else:
+            time_interval_start = 0.0
+
+            
+        for quantity_name, info in self.quantities_to_be_monitored.items():
+            msg += '    %s:\n' %quantity_name
+
+            msg += '      values since time = %.2f in [%s, %s]\n'\
+                   %(time_interval_start,
+                     get_textual_float(info['min'], precision),
+                     get_textual_float(info['max'], precision))                     
+                     
+            msg += '      minimum attained at time = %s, location = %s\n'\
+                   %(get_textual_float(info['min_time'], precision),
+                     get_textual_float(info['min_location'], precision))                     
+                     
+
+            msg += '      maximum attained at time = %s, location = %s\n'\
+                   %(get_textual_float(info['max_time'], precision),
+                     get_textual_float(info['max_location'], precision))                   
+
+
+        return msg
 
 
@@ -902 +1023 @@
         from anuga.config import min_timestep, max_timestep, epsilon
 
-        #FIXME: Maybe lump into a larger check prior to evolving
+        # FIXME: Maybe lump into a larger check prior to evolving
         msg = 'Boundary tags must be bound to boundary objects before '
         msg += 'evolving system, '
@@ -920 +1041 @@
 
         if finaltime is not None and duration is not None:
-            #print 'F', finaltime, duration
+            # print 'F', finaltime, duration
             msg = 'Only one of finaltime and duration may be specified'
             raise msg
@@ -932 +1053 @@
 
 
-        self.yieldtime = 0.0 #Time between 'yields'
-
-        #Initialise interval of timestep sizes (for reporting only)
+        self.yieldtime = 0.0 # Track time between 'yields'
+
+        # Initialise interval of timestep sizes (for reporting only)
         self.min_timestep = max_timestep
         self.max_timestep = min_timestep
@@ -940 +1061 @@
         self.number_of_first_order_steps = 0
 
-        #update ghosts
+        # Update ghosts
         self.update_ghosts()
 
-        #Initial update of vertex and edge values
+        # Initial update of vertex and edge values
         self.distribute_to_vertices_and_edges()
 
-        #Initial update boundary values
+        # Update extrema if necessary (for reporting)
+        self.update_extrema()
+        
+        # Initial update boundary values
         self.update_boundary()
 
-        #Or maybe restore from latest checkpoint
+        # Or maybe restore from latest checkpoint
         if self.checkpoint is True:
             self.goto_latest_checkpoint()
 
         if skip_initial_step is False:
-            yield(self.time)  #Yield initial values
+            yield(self.time)  # Yield initial values
 
         while True:
-            #Compute fluxes across each element edge
+            # Compute fluxes across each element edge
             self.compute_fluxes()
 
-            #Update timestep to fit yieldstep and finaltime
+            # Update timestep to fit yieldstep and finaltime
             self.update_timestep(yieldstep, finaltime)
 
-            #Update conserved quantities
+            # Update conserved quantities
             self.update_conserved_quantities()
 
-            #update ghosts
+            # Update ghosts
             self.update_ghosts()
 
-            #Update vertex and edge values
+            # Update vertex and edge values
             self.distribute_to_vertices_and_edges()
 
-            #Update boundary values
+            # Update extrema if necessary (for reporting)
+            self.update_extrema()            
+
+            # Update boundary values
             self.update_boundary()
 
-            #Update time
+            # Update time
             self.time += self.timestep
             self.yieldtime += self.timestep
@@ -982 +1109 @@
                 self.number_of_first_order_steps += 1
 
-            #Yield results
+            # Yield results
             if finaltime is not None and self.time >= finaltime-epsilon:
 
                 if self.time > finaltime:
-                    #FIXME (Ole, 30 April 2006): Do we need this check?
+                    # FIXME (Ole, 30 April 2006): Do we need this check?
                     print 'WARNING (domain.py): time overshot finaltime. Contact Ole.Nielsen@ga.gov.au'
                     self.time = finaltime

anuga_core/source/anuga/abstract_2d_finite_volumes/quantity.py

@@ -16 +16 @@
 
 from Numeric import array, zeros, Float, less, concatenate, NewAxis,\
-     argmax, allclose, take, reshape
+     argmax, argmin, allclose, take, reshape
 
 from anuga.utilities.numerical_tools import ensure_numeric, is_scalar
@@ -799 +799 @@
 
    
-    def get_maximum_index(self, indices=None):
-        """Return index for maximum value of quantity (on centroids)
-
-        Optional argument:
+    def get_extremum_index(self, mode=None, indices=None):
+        """Return index for maximum or minimum value of quantity (on centroids)
+
+        Optional arguments:
+            mode is either 'max'(default) or 'min'.
             indices is the set of element ids that the operation applies to.
 
         Usage:
-            i = get_maximum_index()
+            i = get_extreme_index()
 
         Notes:
-            We do not seek the maximum at vertices as each vertex can
+            We do not seek the extremum at vertices as each vertex can
             have multiple values - one for each triangle sharing it.
 
@@ -819 +820 @@
 
         # Always return absolute indices
-        i = argmax(V)
-
+        if mode is None or mode == 'max':
+            i = argmax(V)
+        elif mode == 'min':    
+            i = argmin(V)
+
+            
         if indices is None:
             return i
@@ -826 +831 @@
             return indices[i]
 
+
+    def get_maximum_index(self, indices=None):
+        """See get extreme index for details
+        """
+
+        return self.get_extremum_index(mode='max',
+                                       indices=indices)
+
+
         
     def get_maximum_value(self, indices=None):
@@ -866 +880 @@
 
         i = self.get_maximum_index(indices)
+        x, y = self.domain.get_centroid_coordinates()[i]
+
+        return x, y
+
+
+    def get_minimum_index(self, indices=None):
+        """See get extreme index for details
+        """        
+
+        return self.get_extremum_index(mode='min',
+                                       indices=indices)
+
+
+    def get_minimum_value(self, indices=None):
+        """Return minimum value of quantity (on centroids)
+
+        Optional argument:
+            indices is the set of element ids that the operation applies to.
+
+        Usage:
+            v = get_minimum_value()
+
+        See get_maximum_value for more details.    
+        """
+
+
+        i = self.get_minimum_index(indices)
+        V = self.get_values(location='centroids')
+        
+        return V[i]
+        
+
+    def get_minimum_location(self, indices=None):
+        """Return location of minimum value of quantity (on centroids)
+
+        Optional argument:
+            indices is the set of element ids that the operation applies to.
+
+        Usage:
+            x, y = get_minimum_location()
+
+
+        Notes:
+            We do not seek the maximum at vertices as each vertex can
+            have multiple values - one for each triangle sharing it.
+
+            If there are multiple cells with same maximum value, the
+            first cell encountered in the triangle array is returned.       
+        """
+
+        i = self.get_minimum_index(indices)
         x, y = self.domain.get_centroid_coordinates()[i]
 

anuga_core/source/anuga/abstract_2d_finite_volumes/test_domain.py

@@ -219 +219 @@
         assert domain.quantities_to_be_monitored.has_key('stage')
         assert domain.quantities_to_be_monitored.has_key('stage-elevation')
-        assert domain.quantities_to_be_monitored['stage'][0] == None
-        assert domain.quantities_to_be_monitored['stage'][1] == None                
-        
-
-        # Check that invalid request are dealt with
+        for key in domain.quantities_to_be_monitored['stage'].keys():
+            assert domain.quantities_to_be_monitored['stage'][key] is None
+
+
+        # Check that invalid requests are dealt with
         try:
             domain.set_quantities_to_be_monitored(['yyyyy'])        
@@ -238 +238 @@
         else:
             msg = 'Should have caught illegal quantity'
+            raise Exception, msg
+
+        try:
+            domain.set_quantities_to_be_monitored('stage', 'stage-elevation')
+        except:
+            pass
+        else:
+            msg = 'Should have caught too many arguments'
+            raise Exception, msg
+
+        try:
+            domain.set_quantities_to_be_monitored('stage', 'blablabla')
+        except:
+            pass
+        else:
+            msg = 'Should have caught polygon as a string'
             raise Exception, msg        
+
+
+
+        # Now try with a polygon restriction
+        domain.set_quantities_to_be_monitored('xmomentum',
+                                              polygon=[[1,1], [1,3], [3,3], [3,1]],
+                                              time_interval = [0,3])
+        assert domain.monitor_indices[0] == 1
+        assert domain.monitor_time_interval[0] == 0
+        assert domain.monitor_time_interval[1] == 3        
+        
 
     def test_set_quantity_from_expression(self):

anuga_core/source/anuga/abstract_2d_finite_volumes/test_quantity.py

@@ -122 +122 @@
                                                [3.0, -1.5, -1.5]])
 
-    def test_get_maximum_1(self):
+    def test_get_extrema_1(self):
         quantity = Conserved_quantity(self.mesh4,
                                       [[1,2,3], [5,5,5], [0,0,9], [-6, 3, 3]])
@@ -130 +130 @@
         assert v == 5
 
+        v = quantity.get_minimum_value()
+        assert v == 0        
+
         i = quantity.get_maximum_index()
         assert i == 1
+
+        i = quantity.get_minimum_index()
+        assert i == 3        
         
         x,y = quantity.get_maximum_location()
@@ -140 +146 @@
         v = quantity.get_values(interpolation_points = [[x,y]])
         assert allclose(v, 5)
+
+
+        x,y = quantity.get_minimum_location()
+        v = quantity.get_values(interpolation_points = [[x,y]])
+        assert allclose(v, 0)
+
 
     def test_get_maximum_2(self):
@@ -162 +174 @@
         assert v == 6
 
+        v = quantity.get_minimum_value()
+        assert v == 2        
+
         i = quantity.get_maximum_index()
         assert i == 3
+
+        i = quantity.get_minimum_index()
+        assert i == 0        
         
         x,y = quantity.get_maximum_location()
@@ -173 +191 @@
         assert allclose(v, 6)
 
+        x,y = quantity.get_minimum_location()        
+        v = quantity.get_values(interpolation_points = [[x,y]])
+        assert allclose(v, 2)
 
         #Multiple locations for maximum -

anuga_core/source/anuga/abstract_2d_finite_volumes/util.py

@@ -431 +431 @@
         raise ValueError, msg
     
+
+def get_textual_float(value, format = '%.2f'):
+    """Get textual representation of floating point numbers
+    and accept None as valid entry
+
+    format is a string - default = '%.2f'
+    """
+
+    if value is None:
+        return 'None'
+    else:
+        try:
+            float(value)
+        except:
+            # May this is a vector
+            if len(value) > 1:
+                s = '('
+                for v in value:
+                    s += get_textual_float(v, format) + ', '
+                    
+                s = s[:-2] + ')' # Strip trailing comma and close
+                return s
+            else:
+                raise 'Illegal input to get_textual_float:', value
+        else:
+            return format %float(value)
 
 

anuga_core/source/anuga/shallow_water/data_manager.py

@@ -321 +321 @@
 
         if hasattr(domain, 'minimum_storable_height'):
-            self.minimum_storable_height =  domain.minimum_storable_height
+            self.minimum_storable_height = domain.minimum_storable_height
         else:
             self.minimum_storable_height = default_minimum_storable_height
@@ -331 +331 @@
             description = 'Output from anuga.abstract_2d_finite_volumes suitable for plotting'
             self.writer = Write_sww()
-            self.writer.header(fid, domain.starttime,
-                                         self.number_of_volumes,
-                                         self.domain.number_of_full_nodes,
-                                         description=description,
-                                         smoothing=domain.smooth,
-                                         order=domain.default_order)
+            self.writer.store_header(fid,
+                                     domain.starttime,
+                                     self.number_of_volumes,
+                                     self.domain.number_of_full_nodes,
+                                     description=description,
+                                     smoothing=domain.smooth,
+                                     order=domain.default_order)
+
+            # Extra optional information
             if hasattr(domain, 'texture'):
-                fid.texture = domain.texture               
-    #        if domain.geo_reference is not None:
-    #            domain.geo_reference.write_NetCDF(fid)
-            
-#             fid.institution = 'Geoscience Australia'
-#             fid.description = 'Output from anuga.abstract_2d_finite_volumes suitable for plotting'
-
-#             if domain.smooth:
-#                 fid.smoothing = 'Yes'
-#             else:
-#                 fid.smoothing = 'No'
-
-#             fid.order = domain.default_order
-
-#             #Reference point
-#             #Start time in seconds since the epoch (midnight 1/1/1970)
-#             #FIXME: Use Georef
-#             fid.starttime = domain.starttime
-
-#             # dimension definitions
-#             fid.createDimension('number_of_volumes', self.number_of_volumes)
-#             fid.createDimension('number_of_vertices', 3)
-
-#             if domain.smooth is True:
-#                 #fid.createDimension('number_of_points', len(domain.vertexlist))
-#                 fid.createDimension('number_of_points', self.domain.number_of_full_nodes)
-
-#                 # FIXME(Ole): This will cause sww files for paralle domains to
-#                 # have ghost nodes stored (but not used by triangles).
-#                 # To clean this up, we have to change get_vertex_values and friends in
-#                 # quantity.py (but I can't be bothered right now)
-#             else:
-#                 fid.createDimension('number_of_points', 3*self.number_of_volumes)
-
-#             fid.createDimension('number_of_timesteps', None) #extensible
-
-#             # variable definitions
-#             fid.createVariable('x', self.precision, ('number_of_points',))
-#             fid.createVariable('y', self.precision, ('number_of_points',))
-#             fid.createVariable('elevation', self.precision, ('number_of_points',))
-#             if domain.geo_reference is not None:
-#                 domain.geo_reference.write_NetCDF(fid)
-
-#             #FIXME: Backwards compatibility
-#             fid.createVariable('z', self.precision, ('number_of_points',))
-#             #################################
-
-#             fid.createVariable('volumes', Int, ('number_of_volumes',
-#                                                 'number_of_vertices'))
-
-#             fid.createVariable('time', Float,  # Always use full precision lest two timesteps 
-#                                              # close to each other may appear as the same step 
-#                                ('number_of_timesteps',))
-
-#             fid.createVariable('stage', self.precision,
-#                                ('number_of_timesteps',
-#                                 'number_of_points'))
-
-#             fid.createVariable('xmomentum', self.precision,
-#                                ('number_of_timesteps',
-#                                 'number_of_points'))
-
-#             fid.createVariable('ymomentum', self.precision,
-#                                ('number_of_timesteps',
-#                                 'number_of_points'))
-
-        #Close
+                fid.texture = domain.texture
+
+            if domain.quantities_to_be_monitored is not None:
+                fid.createDimension('singleton', 1)
+                for q in domain.quantities_to_be_monitored:
+                    #print 'doing', q
+                    fid.createVariable(q+':extrema', self.precision,
+                                       ('numbers_in_range',))
+                    fid.createVariable(q+':min_location', self.precision,
+                                       ('numbers_in_range',))
+                    fid.createVariable(q+':max_location', self.precision,
+                                       ('numbers_in_range',))
+                    fid.createVariable(q+':min_time', self.precision,
+                                       ('singleton',))
+                    fid.createVariable(q+':max_time', self.precision,
+                                       ('singleton',))
+
+                    
         fid.close()
 
@@ -438 +392 @@
         #
         points = concatenate( (X[:,NewAxis],Y[:,NewAxis]), axis=1 )
-        self.writer.triangulation(fid, points,
-                                            V.astype(volumes.typecode()),
-                                            Z,
-                                            points_georeference= \
-                                            domain.geo_reference)
-        #
-                                
-#         volumes[:] = V.astype(volumes.typecode())
-#         x[:] = X
-#         y[:] = Y
-#         z[:] = Z
-
-#         #FIXME: Backwards compatibility
-#         z = fid.variables['z']
-#         z[:] = Z
-        ################################
-
-
-
-        #Close
+        self.writer.store_triangulation(fid,
+                                        points,
+                                        V.astype(volumes.typecode()),
+                                        Z,
+                                        points_georeference= \
+                                        domain.geo_reference)
+
+        # Close
         fid.close()
+
 
     def store_timestep(self, names):
@@ -470 +413 @@
         from Numeric import choose
         
-        #Get NetCDF       
+        # Get NetCDF       
         retries = 0
         file_open = False
         while not file_open and retries < 10:
             try:
-                fid = NetCDFFile(self.filename, 'a') #Open existing file
+                fid = NetCDFFile(self.filename, 'a') # Open existing file
             except IOError:
-                #This could happen if someone was reading the file.
-                #In that case, wait a while and try again
+                # This could happen if someone was reading the file.
+                # In that case, wait a while and try again
                 msg = 'Warning (store_timestep): File %s could not be opened'\
                       %self.filename
@@ -494 +437 @@
 
 
-        #Check to see if the file is already too big:
+        # Check to see if the file is already too big:
         time = fid.variables['time']
         i = len(time)+1
@@ -500 +443 @@
         file_size_increase =  file_size/i
         if file_size + file_size_increase > self.max_size*(2**self.recursion):
-            #in order to get the file name and start time correct,
-            #I change the domian.filename and domain.starttime.
-            #This is the only way to do this without changing
-            #other modules (I think).
-
-            #write a filename addon that won't break swollens reader
-            #(10.sww is bad)
+            # In order to get the file name and start time correct,
+            # I change the domain.filename and domain.starttime.
+            # This is the only way to do this without changing
+            # other modules (I think).
+
+            # Write a filename addon that won't break swollens reader
+            # (10.sww is bad)
             filename_ext = '_time_%s'%self.domain.time
             filename_ext = filename_ext.replace('.', '_')
-            #remember the old filename, then give domain a
-            #name with the extension
+            
+            # Remember the old filename, then give domain a
+            # name with the extension
             old_domain_filename = self.domain.get_name()
             if not self.recursion:
@@ -516 +460 @@
 
 
-            #change the domain starttime to the current time
+            # Change the domain starttime to the current time
             old_domain_starttime = self.domain.starttime
             self.domain.starttime = self.domain.time
 
-            #build a new data_structure.
+            # Build a new data_structure.
             next_data_structure=\
                 Data_format_sww(self.domain, mode=self.mode,\
@@ -554 +498 @@
 
             if 'stage' in names and 'xmomentum' in names and \
-                   'ymomentum' in names:
+               'ymomentum' in names:
 
                 # Get stage
@@ -573 +517 @@
                 ymomentum, _ = Q.get_vertex_values(xy = False,
                                           precision = self.precision)
-                self.writer.quantities(fid, 
-                                                 time=self.domain.time,
-                                                 precision=self.precision,
-                                                 stage=stage,
-                                                 xmomentum=xmomentum,
-                                                 ymomentum=ymomentum)
+                
+                # Write quantities to NetCDF
+                self.writer.store_quantities(fid, 
+                                             time=self.domain.time,
+                                             precision=self.precision,
+                                             stage=stage,
+                                             xmomentum=xmomentum,
+                                             ymomentum=ymomentum)
             else:
                 # This is producing a sww that is not standard.
-                #Store time
+                # Store time
                 time[i] = self.domain.time
                 
@@ -588 +534 @@
                     Q = domain.quantities[name]
                     A,V = Q.get_vertex_values(xy = False,
-                                          precision = self.precision)
-
-                    #FIXME: Make this general (see below)
+                                              precision = self.precision)
+
+                    # FIXME: Make this general (see below)
                     if name == 'stage':
                         z = fid.variables['elevation']
@@ -603 +549 @@
                    #As in....
                    #eval( name + '[i,:] = A.astype(self.precision)' )
-                   #FIXME: But we need a UNIT test for that before
+                   #FIXME (Ole): But we need a UNIT test for that before
                    # refactoring
 
 
+
+            # Update extrema if requested
+            domain = self.domain
+            if domain.quantities_to_be_monitored is not None:
+                for q, info in domain.quantities_to_be_monitored.items():
+
+                    if info['min'] is not None:
+                        fid.variables[q + ':extrema'][0] = info['min']
+                        fid.variables[q + ':min_location'][:] =\
+                                        info['min_location']
+                        fid.variables[q + ':min_time'][0] = info['min_time']
+                        
+                    if info['max'] is not None:
+                        fid.variables[q + ':extrema'][1] = info['max']
+                        fid.variables[q + ':max_location'][:] =\
+                                        info['max_location']
+                        fid.variables[q + ':max_time'][0] = info['max_time']
+
+            
 
             #Flush and close
@@ -614 +579 @@
 
 
-#Class for handling checkpoints data
+# Class for handling checkpoints data
 class Data_format_cpt(Data_format):
     """Interface to native NetCDF format (.cpt)
@@ -752 +717 @@
         Q = domain.quantities[name]
         A,V = Q.get_vertex_values(xy=False,
-                  precision = self.precision)
+                                  precision = self.precision)
 
         stage[i,:] = A.astype(self.precision)
@@ -2863 +2828 @@
 
 
-    #print number_of_latitudes, number_of_longitudes
+    # print number_of_latitudes, number_of_longitudes
     number_of_points = number_of_latitudes*number_of_longitudes
     number_of_volumes = (number_of_latitudes-1)*(number_of_longitudes-1)*2
@@ -2883 +2848 @@
                     basename_in + '_e.nc')
     
-    #Create new file
+    # Create new file
     starttime = times[0]
     sww = Write_sww()
-    sww.header(outfile, times, number_of_volumes,
-                         number_of_points, description=description,
-                         verbose=verbose)
-
-    #Store
+    sww.store_header(outfile, times, number_of_volumes,
+                     number_of_points, description=description,
+                     verbose=verbose)
+
+    # Store
     from anuga.coordinate_transforms.redfearn import redfearn
     x = zeros(number_of_points, Float)  #Easting
@@ -2897 +2862 @@
 
     if verbose: print 'Making triangular grid'
-    #Check zone boundaries
+
+    # Check zone boundaries
     refzone, _, _ = redfearn(latitudes[0],longitudes[0])
 
@@ -4727 +4693 @@
     # work out sww_times and the index range this covers
     sww = Write_sww()
-    sww.header(outfile, times, len(volumes), len(points_utm),
-                         verbose=verbose)
+    sww.store_header(outfile, times, len(volumes), len(points_utm),
+                     verbose=verbose)
     outfile.mean_stage = mean_stage
     outfile.zscale = zscale
 
-    sww.triangulation(outfile, points_utm, volumes,
+    sww.store_triangulation(outfile, points_utm, volumes,
                             elevation, zone,  new_origin=origin,
                             verbose=verbose)
@@ -4745 +4711 @@
             xmomentum = ua*h
             ymomentum = -1*va*h # -1 since in mux files south is positive.
-            sww.quantities(outfile, 
-                                     slice_index=j - mux_times_start_i,
-                                     verbose=verbose,
-                                     stage=stage,
-                                     xmomentum=xmomentum,
-                                     ymomentum=ymomentum)
+            sww.store_quantities(outfile, 
+                                 slice_index=j - mux_times_start_i,
+                                 verbose=verbose,
+                                 stage=stage,
+                                 xmomentum=xmomentum,
+                                 ymomentum=ymomentum)
         j += 1
     if verbose: sww.verbose_quantities(outfile)
@@ -4778 +4744 @@
 
 class Write_sww:
-    from anuga.shallow_water.shallow_water_domain import Domain 
+    from anuga.shallow_water.shallow_water_domain import Domain
+
+    # FIXME (Ole): Hardwiring the conserved quantities like
+    # this could be a problem. I would prefer taking them from
+    # the instantiation of Domain.
     sww_quantities = Domain.conserved_quantities
+
+
     RANGE = '_range'
+    EXTREMA = ':extrema'
 
     def __init__(self):
         pass
     
-    def header(self,outfile, times, number_of_volumes,
-                         number_of_points,
-                         description='Converted from XXX',
-                         smoothing=True,
-                         order=1, verbose=False):
+    def store_header(self,
+                     outfile,
+                     times,
+                     number_of_volumes,
+                     number_of_points,
+                     description='Converted from XXX',
+                     smoothing=True,
+                     order=1, verbose=False):
         """
         outfile - the name of the file that will be written
@@ -4863 +4839 @@
                                ('numbers_in_range',))
 
+
         # Initialise ranges with small and large sentinels.
         # If this was in pure Python we could have used None sensibly
@@ -4868 +4845 @@
         outfile.variables[q+Write_sww.RANGE][1] = -max_float # Max
 
-        #FIXME: Backwards compatibility
+        # FIXME: Backwards compatibility
         outfile.createVariable('z', precision, ('number_of_points',))
         #################################
@@ -4900 +4877 @@
 
         
-    def triangulation(self,outfile, points_utm, volumes,
-                                elevation, zone=None, new_origin=None, 
-                                points_georeference=None, verbose=False):
+    def store_triangulation(self,
+                            outfile,
+                            points_utm,
+                            volumes,
+                            elevation, zone=None, new_origin=None, 
+                            points_georeference=None, verbose=False):
         """
         
@@ -4994 +4974 @@
         outfile.variables[q+Write_sww.RANGE][1] = max(elevation)
 
-    def quantities(self, outfile, precision=Float,
-                             slice_index=None, time=None,
-                             verbose=False, **quant):
+
+    def store_quantities(self, outfile, precision=Float,
+                         slice_index=None, time=None,
+                         verbose=False, **quant):
         """
         Write the quantity info.
@@ -5021 +5002 @@
             file_time[slice_index] = time    
 
-        # write the conserved quantities from Doamin.
+        # write the conserved quantities from Domain.
         # Typically stage,  xmomentum, ymomentum
         # other quantities will be ignored, silently.
@@ -5049 +5030 @@
                                        outfile.variables[q+Write_sww.RANGE][1])
         print '------------------------------------------------'
+
 
         

anuga_core/source/anuga/shallow_water/test_data_manager.py

@@ -256 +256 @@
         os.remove(sww.filename)
 
-    def NOtest_sww_extrema(self):
+    def test_sww_extrema(self):
         """Test that extrema of quantities can be retrieved at every vertex
         Extrema are updated at every *internal* timestep
@@ -263 +263 @@
         domain = self.domain
         
-        domain.set_name('datatest' + str(id(self)))
+        domain.set_name('extrema_test' + str(id(self)))
         domain.format = 'sww'
         domain.smooth = True
@@ -271 +271 @@
         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'        
+        domain.set_quantities_to_be_monitored(['xmomentum',
+                                               'ymomentum',
+                                               'stage-elevation'])
+        
+        assert len(domain.quantities_to_be_monitored) == 3
+        assert domain.quantities_to_be_monitored.has_key('stage-elevation')
+        assert domain.quantities_to_be_monitored.has_key('xmomentum')                
+        assert domain.quantities_to_be_monitored.has_key('ymomentum')        
         assert domain.monitor_polygon is None
         assert domain.monitor_time_interval is None        
@@ -282 +285 @@
 
         for t in domain.evolve(yieldstep = 1, finaltime = 1):
-            print domain.timestepping_statistics()
-            print domain.quantity_statistics()
+            pass
+            #print domain.timestepping_statistics()
+            #print domain.quantity_statistics(precision = '%.8f')
 
             
@@ -290 +294 @@
 
         # 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,[])
+        extrema = fid.variables['stage-elevation:extrema'][:]
+        assert allclose(extrema, [0.00, 0.30])
+
+        loc = fid.variables['stage-elevation:min_location'][:]
+        assert allclose(loc, [0.16666667, 0.33333333])
+
+        loc = fid.variables['stage-elevation:max_location'][:]        
+        assert allclose(loc, [0.8333333, 0.16666667])        
+
+        time = fid.variables['stage-elevation:max_time'][:]
+        assert allclose(time, 0.0)                
+
+        extrema = fid.variables['xmomentum:extrema'][:]
+        assert allclose(extrema,[-0.06062178, 0.47886313])
+        
+        extrema = fid.variables['ymomentum:extrema'][:]
+        assert allclose(extrema,[0.00, 0.06241221])        
 
         
@@ -6723 +6736 @@
         number_of_points = len(points_utm)
         sww = Write_sww()
-        sww.header(outfile, times, number_of_volumes,
+        sww.store_header(outfile, times, number_of_volumes,
                          number_of_points, description='fully sick testing',
-                             verbose=self.verbose)
-        sww.triangulation(outfile, points_utm, volumes,
-                                    elevation,  new_origin=new_origin,
-                                    verbose=self.verbose)       
+                         verbose=self.verbose)
+        sww.store_triangulation(outfile, points_utm, volumes,
+                                elevation,  new_origin=new_origin,
+                                verbose=self.verbose)       
         outfile.close()
         fid = NetCDFFile(filename)
@@ -6755 +6768 @@
         number_of_points = len(points_utm)
         sww = Write_sww()
-        sww.header(outfile, times, number_of_volumes,
+        sww.store_header(outfile, times, number_of_volumes,
                          number_of_points, description='fully sick testing',
                          verbose=self.verbose)
-        sww.triangulation(outfile, points_utm, volumes,
-                                    elevation,  new_origin=new_origin,
-                                    verbose=self.verbose)       
+        sww.store_triangulation(outfile, points_utm, volumes,
+                                elevation,  new_origin=new_origin,
+                                verbose=self.verbose)       
         outfile.close()
         fid = NetCDFFile(filename)
@@ -6791 +6804 @@
         number_of_points = len(points_utm)
         sww = Write_sww()
-        sww.header(outfile, times, number_of_volumes,
+        sww.store_header(outfile, times, number_of_volumes,
                          number_of_points, description='fully sick testing',
                          verbose=self.verbose)
-        sww.triangulation(outfile, points_utm, volumes,
-                                    elevation,  new_origin=new_origin,
-                                    verbose=self.verbose)       
+        sww.store_triangulation(outfile, points_utm, volumes,
+                                elevation,  new_origin=new_origin,
+                                verbose=self.verbose)       
         outfile.close()
         fid = NetCDFFile(filename)
@@ -6830 +6843 @@
         number_of_points = len(points_utm)
         sww = Write_sww()
-        sww.header(outfile, times, number_of_volumes,
+        sww.store_header(outfile, times, number_of_volumes,
                          number_of_points, description='fully sick testing',
                          verbose=self.verbose)
-        sww.triangulation(outfile, points_utm, volumes,
-                                    elevation,  new_origin=new_origin,
-                                    points_georeference=points_georeference,
-                                    verbose=self.verbose)       
+        sww.store_triangulation(outfile, points_utm, volumes,
+                                elevation,  new_origin=new_origin,
+                                points_georeference=points_georeference,
+                                verbose=self.verbose)       
         outfile.close()
         fid = NetCDFFile(filename)
@@ -6866 +6879 @@
         number_of_points = len(points_utm)
         sww = Write_sww()
-        sww.header(outfile, times, number_of_volumes,
+        sww.store_header(outfile, times, number_of_volumes,
                          number_of_points, description='fully sick testing',
                          verbose=self.verbose)
-        sww.triangulation(outfile, points_utm, volumes,
-                                    elevation,  new_origin=new_origin,
-                                    points_georeference=points_georeference,
-                                    verbose=self.verbose)       
+        sww.store_triangulation(outfile, points_utm, volumes,
+                                elevation,  new_origin=new_origin,
+                                points_georeference=points_georeference,
+                                verbose=self.verbose)       
         outfile.close()
         fid = NetCDFFile(filename)
@@ -7281 +7294 @@
     #suite = unittest.makeSuite(Test_Data_Manager,'test_sww_header')
     #suite = unittest.makeSuite(Test_Data_Manager,'test_export_grid_parallel')
-    #suite = unittest.makeSuite(Test_Data_Manager,'t')
     suite = unittest.makeSuite(Test_Data_Manager,'test')
+    #suite = unittest.makeSuite(Test_Data_Manager,'test_sww_extrema')
 
     

anuga_core/source/anuga/shallow_water/test_shallow_water_domain.py

@@ -4635 +4635 @@
 
 
+
+    def test_extrema(self):
+        """Test that extrema of quantities are computed correctly
+        Extrema are updated at every *internal* timestep
+        """
+
+        from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular_cross
+
+        initial_runup_height = -0.4
+        final_runup_height = -0.3
+
+
+        #--------------------------------------------------------------
+        # Setup computational domain
+        #--------------------------------------------------------------
+        N = 5
+        points, vertices, boundary = rectangular_cross(N, N) 
+        domain = Domain(points, vertices, boundary)
+        domain.set_name('extrema_test')
+
+        #--------------------------------------------------------------
+        # 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
+        domain.set_quantities_to_be_monitored(['stage', 'stage-elevation'],
+                                              time_interval = [0.5, 2.7],
+                                              polygon = [[0,0], [0,1], [1,1], [1,0]])
+        
+        assert len(domain.quantities_to_be_monitored) == 2
+        assert domain.quantities_to_be_monitored.has_key('stage')
+        assert domain.quantities_to_be_monitored.has_key('stage-elevation')
+        for key in domain.quantities_to_be_monitored['stage'].keys():
+            assert domain.quantities_to_be_monitored['stage'][key] is None        
+
+
+        #--------------------------------------------------------------
+        # 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})
+
+
+        #--------------------------------------------------------------
+        # Let triangles adjust and check extrema 
+        #--------------------------------------------------------------
+        for t in domain.evolve(yieldstep = 0.1, finaltime = 1.0):
+            domain.quantity_statistics() # Run it silently
+
+
+
+        #--------------------------------------------------------------
+        # Test extrema
+        #--------------------------------------------------------------
+
+        stage = domain.quantities_to_be_monitored['stage']
+        assert stage['min'] <= stage['max']
+
+        #print stage['min'], stage['max'] 
+        assert allclose(stage['min'], initial_runup_height,
+                        rtol = 1.0/N) # First order accuracy
+
+
+        depth = domain.quantities_to_be_monitored['stage-elevation']
+        assert depth['min'] <= depth['max'] 
+        assert depth['min'] >= 0.0
+        assert depth['max'] >= 0.0        
+        ##assert depth[1] <= ?? initial_runup_height        
+
+
+        #--------------------------------------------------------------
+        # Update boundary to allow inflow
+        #--------------------------------------------------------------
+        domain.set_boundary({'right': Bd})
+
+        
+        #--------------------------------------------------------------
+        # Evolve system through time
+        #--------------------------------------------------------------
+        for t in domain.evolve(yieldstep = 0.1, finaltime = 3.0):
+            #domain.write_time()
+            domain.quantity_statistics() # Run it silently           
+            
+    
+        #--------------------------------------------------------------
+        # Test extrema again
+        #--------------------------------------------------------------
+
+        stage = domain.quantities_to_be_monitored['stage']
+        assert stage['min'] <= stage['max']
+
+        assert allclose(stage['min'], initial_runup_height,
+                        rtol = 1.0/N) # First order accuracy        
+
+        depth = domain.quantities_to_be_monitored['stage-elevation']
+        assert depth['min'] <= depth['max'] 
+        assert depth['min'] >= 0.0
+        assert depth['max'] >= 0.0        
+        
+
+
     def test_tight_slope_limiters(self):
         """Test that new slope limiters (Feb 2007) don't induce extremely
@@ -4872 +4983 @@
 if __name__ == "__main__":
 
-    suite = unittest.makeSuite(Test_Shallow_Water,'test')    
+    #suite = unittest.makeSuite(Test_Shallow_Water,'test')    
+    suite = unittest.makeSuite(Test_Shallow_Water,'test_extrema')    
     #suite = unittest.makeSuite(Test_Shallow_Water,'test_tight_slope_limiters')
     #suite = unittest.makeSuite(Test_Shallow_Water,'test_get_maximum_inundation_from_sww')