Changeset 4704


Ignore:
Timestamp:
Sep 5, 2007, 4:42:58 PM (17 years ago)
Author:
ole
Message:

Implemented ticket:192

Location:
anuga_core/source/anuga
Files:
8 edited

Legend:

Unmodified
Added
Removed
  • anuga_core/source/anuga/abstract_2d_finite_volumes/domain.py

    r4702 r4704  
    2626from anuga.abstract_2d_finite_volumes.region\
    2727     import Set_region as region_set_region
     28
     29from anuga.utilities.polygon import inside_polygon
     30from anuga.abstract_2d_finite_volumes.util import get_textual_float
    2831
    2932import types
     
    205208        self.monitor_polygon = None
    206209        self.monitor_time_interval = None               
    207 
     210        self.monitor_indices = None
    208211
    209212        # Checkpointing and storage
     
    579582            self.quantities_to_be_monitored = None
    580583            self.monitor_polygon = None
    581             self.monitor_time_interval = None                   
     584            self.monitor_time_interval = None
     585            self.monitor_indices = None           
    582586            return
    583587
     
    596600                apply_expression_to_dictionary(quantity_name, self.quantities)
    597601
    598             # Initialise extrema
    599             self.quantities_to_be_monitored[quantity_name] = [None, None]
     602            # Initialise extrema information
     603            info_block = {'min': None,          # Min value
     604                          'max': None,          # Max value
     605                          'min_location': None, # Argmin (x, y)
     606                          'max_location': None, # Argmax (x, y)
     607                          'min_time': None,     # Argmin (t)
     608                          'max_time': None}     # Argmax (t)
    600609           
     610            self.quantities_to_be_monitored[quantity_name] = info_block
     611
     612           
    601613
    602614        if polygon is not None:
    603             # FIXME Check input
    604             pass
     615            # Check input
     616            if isinstance(polygon, basestring):
     617
     618                # Check if multiple quantities were accidentally
     619                # given as separate argument rather than a list.
     620                msg = 'Multiple quantities must be specified in a list. '
     621                msg += 'Not as multiple arguments. '
     622                msg += 'I got "%s" as a second argument' %polygon
     623               
     624                if polygon in self.quantities:
     625                    raise Exception, msg
     626               
     627                try:
     628                    apply_expression_to_dictionary(polygon, self.quantities)
     629                except:
     630                    # At least polygon wasn't an expression involving quantitites
     631                    pass
     632                else:
     633                    raise Exception, msg
     634
     635                # In any case, we don't allow polygon to be a string
     636                msg = 'argument "polygon" must not be a string: '
     637                msg += 'I got polygon=\'%s\' ' %polygon
     638                raise Exception, msg
     639
     640
     641            # Get indices for centroids that are inside polygon
     642            points = self.get_centroid_coordinates(absolute=True)
     643            self.monitor_indices = inside_polygon(points, polygon)
     644           
    605645
    606646        if time_interval is not None:
    607             # FIXME Check input
    608             pass       
    609 
     647            assert len(time_interval) == 2
    610648
    611649       
     
    758796        """
    759797
    760         #Input checks
     798        # Input checks
    761799        import types, string
    762800
     
    780818        assert type(tags) == types.ListType, msg
    781819
    782         #Determine width of longest quantity name (for cosmetic purposes)
     820        # Determine width of longest quantity name (for cosmetic purposes)
    783821        maxwidth = 0
    784822        for name in quantities:
     
    787825                maxwidth = w
    788826
    789         #Output stats
     827        # Output stats
    790828        msg = 'Boundary values at time %.4f:\n' %self.time
    791829        for tag in tags:
     
    795833                q = self.quantities[name]
    796834
    797                 #Find range of boundary values for tag and q
     835                # Find range of boundary values for tag and q
    798836                maxval = minval = None
    799837                for i, ((vol_id, edge_id), B) in\
     
    815853
    816854
    817 
    818     def quantity_statistics(self):
     855    def update_extrema(self):
     856        """Update extrema if requested by set_quantities_to_be_monitored.
     857        This data is used for reporting e.g. by running
     858        print domain.quantity_statistics()
     859        and may also stored in output files (see data_manager in shallow_water)
     860        """
     861
     862        if self.quantities_to_be_monitored is None:
     863            return
     864
     865        # Observe time interval restriction if any
     866        if self.monitor_time_interval is not None and\
     867               (self.time < self.monitor_time_interval[0] or\
     868               self.time > self.monitor_time_interval[1]):
     869            return
     870           
     871       
     872        # Update extrema for each specified quantity subject to
     873        # polygon restriction (via monitor_indices).
     874        for quantity_name in self.quantities_to_be_monitored:
     875
     876            if quantity_name in self.quantities:
     877                Q = self.get_quantity(quantity_name)
     878            else:
     879                Q = self.create_quantity_from_expression(quantity_name)
     880
     881            info_block = self.quantities_to_be_monitored[quantity_name]
     882
     883            # Update maximum (n > None is always True)
     884            maxval = Q.get_maximum_value(self.monitor_indices)
     885            if maxval > info_block['max']:
     886                info_block['max'] = maxval
     887                maxloc = Q.get_maximum_location()
     888                info_block['max_location'] = maxloc
     889                info_block['max_time'] = self.time
     890
     891
     892            # Update minimum (n < None is always False)
     893            minval = Q.get_minimum_value(self.monitor_indices)
     894            if info_block['min'] is None or\
     895                   minval < info_block['min']:
     896                info_block['min'] = minval               
     897                minloc = Q.get_minimum_location()
     898                info_block['min_location'] = minloc
     899                info_block['min_time'] = self.time               
     900       
     901
     902
     903    def quantity_statistics(self, precision = '%.4f'):
    819904        """Return string with statistics about quantities for printing or logging
    820905
     
    825910        """
    826911
    827         pass
    828 
    829 
     912        maxlen = 128 # Max length of polygon string representation
     913
     914        # Output stats
     915        msg = 'Monitored quantities at time %.4f:\n' %self.time
     916        if self.monitor_polygon is not None:
     917            p_str = str(self.monitor_polygon)
     918            msg += '- Restricted by polygon: %s' %p_str[:maxlen]
     919            if len(p_str) >= maxlen:
     920                msg += '...\n'
     921            else:
     922                msg += '\n'
     923
     924
     925        if self.monitor_time_interval is not None:
     926            msg += '- Restricted by time interval: %s\n' %str(self.monitor_time_interval)
     927            time_interval_start = self.monitor_time_interval[0]
     928        else:
     929            time_interval_start = 0.0
     930
     931           
     932        for quantity_name, info in self.quantities_to_be_monitored.items():
     933            msg += '    %s:\n' %quantity_name
     934
     935            msg += '      values since time = %.2f in [%s, %s]\n'\
     936                   %(time_interval_start,
     937                     get_textual_float(info['min'], precision),
     938                     get_textual_float(info['max'], precision))                     
     939                     
     940            msg += '      minimum attained at time = %s, location = %s\n'\
     941                   %(get_textual_float(info['min_time'], precision),
     942                     get_textual_float(info['min_location'], precision))                     
     943                     
     944
     945            msg += '      maximum attained at time = %s, location = %s\n'\
     946                   %(get_textual_float(info['max_time'], precision),
     947                     get_textual_float(info['max_location'], precision))                   
     948
     949
     950        return msg
    830951
    831952
     
    9021023        from anuga.config import min_timestep, max_timestep, epsilon
    9031024
    904         #FIXME: Maybe lump into a larger check prior to evolving
     1025        # FIXME: Maybe lump into a larger check prior to evolving
    9051026        msg = 'Boundary tags must be bound to boundary objects before '
    9061027        msg += 'evolving system, '
     
    9201041
    9211042        if finaltime is not None and duration is not None:
    922             #print 'F', finaltime, duration
     1043            # print 'F', finaltime, duration
    9231044            msg = 'Only one of finaltime and duration may be specified'
    9241045            raise msg
     
    9321053
    9331054
    934         self.yieldtime = 0.0 #Time between 'yields'
    935 
    936         #Initialise interval of timestep sizes (for reporting only)
     1055        self.yieldtime = 0.0 # Track time between 'yields'
     1056
     1057        # Initialise interval of timestep sizes (for reporting only)
    9371058        self.min_timestep = max_timestep
    9381059        self.max_timestep = min_timestep
     
    9401061        self.number_of_first_order_steps = 0
    9411062
    942         #update ghosts
     1063        # Update ghosts
    9431064        self.update_ghosts()
    9441065
    945         #Initial update of vertex and edge values
     1066        # Initial update of vertex and edge values
    9461067        self.distribute_to_vertices_and_edges()
    9471068
    948         #Initial update boundary values
     1069        # Update extrema if necessary (for reporting)
     1070        self.update_extrema()
     1071       
     1072        # Initial update boundary values
    9491073        self.update_boundary()
    9501074
    951         #Or maybe restore from latest checkpoint
     1075        # Or maybe restore from latest checkpoint
    9521076        if self.checkpoint is True:
    9531077            self.goto_latest_checkpoint()
    9541078
    9551079        if skip_initial_step is False:
    956             yield(self.time)  #Yield initial values
     1080            yield(self.time)  # Yield initial values
    9571081
    9581082        while True:
    959             #Compute fluxes across each element edge
     1083            # Compute fluxes across each element edge
    9601084            self.compute_fluxes()
    9611085
    962             #Update timestep to fit yieldstep and finaltime
     1086            # Update timestep to fit yieldstep and finaltime
    9631087            self.update_timestep(yieldstep, finaltime)
    9641088
    965             #Update conserved quantities
     1089            # Update conserved quantities
    9661090            self.update_conserved_quantities()
    9671091
    968             #update ghosts
     1092            # Update ghosts
    9691093            self.update_ghosts()
    9701094
    971             #Update vertex and edge values
     1095            # Update vertex and edge values
    9721096            self.distribute_to_vertices_and_edges()
    9731097
    974             #Update boundary values
     1098            # Update extrema if necessary (for reporting)
     1099            self.update_extrema()           
     1100
     1101            # Update boundary values
    9751102            self.update_boundary()
    9761103
    977             #Update time
     1104            # Update time
    9781105            self.time += self.timestep
    9791106            self.yieldtime += self.timestep
     
    9821109                self.number_of_first_order_steps += 1
    9831110
    984             #Yield results
     1111            # Yield results
    9851112            if finaltime is not None and self.time >= finaltime-epsilon:
    9861113
    9871114                if self.time > finaltime:
    988                     #FIXME (Ole, 30 April 2006): Do we need this check?
     1115                    # FIXME (Ole, 30 April 2006): Do we need this check?
    9891116                    print 'WARNING (domain.py): time overshot finaltime. Contact Ole.Nielsen@ga.gov.au'
    9901117                    self.time = finaltime
  • anuga_core/source/anuga/abstract_2d_finite_volumes/quantity.py

    r4679 r4704  
    1616
    1717from Numeric import array, zeros, Float, less, concatenate, NewAxis,\
    18      argmax, allclose, take, reshape
     18     argmax, argmin, allclose, take, reshape
    1919
    2020from anuga.utilities.numerical_tools import ensure_numeric, is_scalar
     
    799799
    800800   
    801     def get_maximum_index(self, indices=None):
    802         """Return index for maximum value of quantity (on centroids)
    803 
    804         Optional argument:
     801    def get_extremum_index(self, mode=None, indices=None):
     802        """Return index for maximum or minimum value of quantity (on centroids)
     803
     804        Optional arguments:
     805            mode is either 'max'(default) or 'min'.
    805806            indices is the set of element ids that the operation applies to.
    806807
    807808        Usage:
    808             i = get_maximum_index()
     809            i = get_extreme_index()
    809810
    810811        Notes:
    811             We do not seek the maximum at vertices as each vertex can
     812            We do not seek the extremum at vertices as each vertex can
    812813            have multiple values - one for each triangle sharing it.
    813814
     
    819820
    820821        # Always return absolute indices
    821         i = argmax(V)
    822 
     822        if mode is None or mode == 'max':
     823            i = argmax(V)
     824        elif mode == 'min':   
     825            i = argmin(V)
     826
     827           
    823828        if indices is None:
    824829            return i
     
    826831            return indices[i]
    827832
     833
     834    def get_maximum_index(self, indices=None):
     835        """See get extreme index for details
     836        """
     837
     838        return self.get_extremum_index(mode='max',
     839                                       indices=indices)
     840
     841
    828842       
    829843    def get_maximum_value(self, indices=None):
     
    866880
    867881        i = self.get_maximum_index(indices)
     882        x, y = self.domain.get_centroid_coordinates()[i]
     883
     884        return x, y
     885
     886
     887    def get_minimum_index(self, indices=None):
     888        """See get extreme index for details
     889        """       
     890
     891        return self.get_extremum_index(mode='min',
     892                                       indices=indices)
     893
     894
     895    def get_minimum_value(self, indices=None):
     896        """Return minimum value of quantity (on centroids)
     897
     898        Optional argument:
     899            indices is the set of element ids that the operation applies to.
     900
     901        Usage:
     902            v = get_minimum_value()
     903
     904        See get_maximum_value for more details.   
     905        """
     906
     907
     908        i = self.get_minimum_index(indices)
     909        V = self.get_values(location='centroids')
     910       
     911        return V[i]
     912       
     913
     914    def get_minimum_location(self, indices=None):
     915        """Return location of minimum value of quantity (on centroids)
     916
     917        Optional argument:
     918            indices is the set of element ids that the operation applies to.
     919
     920        Usage:
     921            x, y = get_minimum_location()
     922
     923
     924        Notes:
     925            We do not seek the maximum at vertices as each vertex can
     926            have multiple values - one for each triangle sharing it.
     927
     928            If there are multiple cells with same maximum value, the
     929            first cell encountered in the triangle array is returned.       
     930        """
     931
     932        i = self.get_minimum_index(indices)
    868933        x, y = self.domain.get_centroid_coordinates()[i]
    869934
  • anuga_core/source/anuga/abstract_2d_finite_volumes/test_domain.py

    r4702 r4704  
    219219        assert domain.quantities_to_be_monitored.has_key('stage')
    220220        assert domain.quantities_to_be_monitored.has_key('stage-elevation')
    221         assert domain.quantities_to_be_monitored['stage'][0] == None
    222         assert domain.quantities_to_be_monitored['stage'][1] == None               
    223        
    224 
    225         # Check that invalid request are dealt with
     221        for key in domain.quantities_to_be_monitored['stage'].keys():
     222            assert domain.quantities_to_be_monitored['stage'][key] is None
     223
     224
     225        # Check that invalid requests are dealt with
    226226        try:
    227227            domain.set_quantities_to_be_monitored(['yyyyy'])       
     
    238238        else:
    239239            msg = 'Should have caught illegal quantity'
     240            raise Exception, msg
     241
     242        try:
     243            domain.set_quantities_to_be_monitored('stage', 'stage-elevation')
     244        except:
     245            pass
     246        else:
     247            msg = 'Should have caught too many arguments'
     248            raise Exception, msg
     249
     250        try:
     251            domain.set_quantities_to_be_monitored('stage', 'blablabla')
     252        except:
     253            pass
     254        else:
     255            msg = 'Should have caught polygon as a string'
    240256            raise Exception, msg       
     257
     258
     259
     260        # Now try with a polygon restriction
     261        domain.set_quantities_to_be_monitored('xmomentum',
     262                                              polygon=[[1,1], [1,3], [3,3], [3,1]],
     263                                              time_interval = [0,3])
     264        assert domain.monitor_indices[0] == 1
     265        assert domain.monitor_time_interval[0] == 0
     266        assert domain.monitor_time_interval[1] == 3       
     267       
    241268
    242269    def test_set_quantity_from_expression(self):
  • anuga_core/source/anuga/abstract_2d_finite_volumes/test_quantity.py

    r4592 r4704  
    122122                                               [3.0, -1.5, -1.5]])
    123123
    124     def test_get_maximum_1(self):
     124    def test_get_extrema_1(self):
    125125        quantity = Conserved_quantity(self.mesh4,
    126126                                      [[1,2,3], [5,5,5], [0,0,9], [-6, 3, 3]])
     
    130130        assert v == 5
    131131
     132        v = quantity.get_minimum_value()
     133        assert v == 0       
     134
    132135        i = quantity.get_maximum_index()
    133136        assert i == 1
     137
     138        i = quantity.get_minimum_index()
     139        assert i == 3       
    134140       
    135141        x,y = quantity.get_maximum_location()
     
    140146        v = quantity.get_values(interpolation_points = [[x,y]])
    141147        assert allclose(v, 5)
     148
     149
     150        x,y = quantity.get_minimum_location()
     151        v = quantity.get_values(interpolation_points = [[x,y]])
     152        assert allclose(v, 0)
     153
    142154
    143155    def test_get_maximum_2(self):
     
    162174        assert v == 6
    163175
     176        v = quantity.get_minimum_value()
     177        assert v == 2       
     178
    164179        i = quantity.get_maximum_index()
    165180        assert i == 3
     181
     182        i = quantity.get_minimum_index()
     183        assert i == 0       
    166184       
    167185        x,y = quantity.get_maximum_location()
     
    173191        assert allclose(v, 6)
    174192
     193        x,y = quantity.get_minimum_location()       
     194        v = quantity.get_values(interpolation_points = [[x,y]])
     195        assert allclose(v, 2)
    175196
    176197        #Multiple locations for maximum -
  • anuga_core/source/anuga/abstract_2d_finite_volumes/util.py

    r4663 r4704  
    431431        raise ValueError, msg
    432432   
     433
     434def get_textual_float(value, format = '%.2f'):
     435    """Get textual representation of floating point numbers
     436    and accept None as valid entry
     437
     438    format is a string - default = '%.2f'
     439    """
     440
     441    if value is None:
     442        return 'None'
     443    else:
     444        try:
     445            float(value)
     446        except:
     447            # May this is a vector
     448            if len(value) > 1:
     449                s = '('
     450                for v in value:
     451                    s += get_textual_float(v, format) + ', '
     452                   
     453                s = s[:-2] + ')' # Strip trailing comma and close
     454                return s
     455            else:
     456                raise 'Illegal input to get_textual_float:', value
     457        else:
     458            return format %float(value)
    433459
    434460
  • anuga_core/source/anuga/shallow_water/data_manager.py

    r4699 r4704  
    321321
    322322        if hasattr(domain, 'minimum_storable_height'):
    323             self.minimum_storable_height =  domain.minimum_storable_height
     323            self.minimum_storable_height = domain.minimum_storable_height
    324324        else:
    325325            self.minimum_storable_height = default_minimum_storable_height
     
    331331            description = 'Output from anuga.abstract_2d_finite_volumes suitable for plotting'
    332332            self.writer = Write_sww()
    333             self.writer.header(fid, domain.starttime,
    334                                          self.number_of_volumes,
    335                                          self.domain.number_of_full_nodes,
    336                                          description=description,
    337                                          smoothing=domain.smooth,
    338                                          order=domain.default_order)
     333            self.writer.store_header(fid,
     334                                     domain.starttime,
     335                                     self.number_of_volumes,
     336                                     self.domain.number_of_full_nodes,
     337                                     description=description,
     338                                     smoothing=domain.smooth,
     339                                     order=domain.default_order)
     340
     341            # Extra optional information
    339342            if hasattr(domain, 'texture'):
    340                 fid.texture = domain.texture               
    341     #        if domain.geo_reference is not None:
    342     #            domain.geo_reference.write_NetCDF(fid)
    343            
    344 #             fid.institution = 'Geoscience Australia'
    345 #             fid.description = 'Output from anuga.abstract_2d_finite_volumes suitable for plotting'
    346 
    347 #             if domain.smooth:
    348 #                 fid.smoothing = 'Yes'
    349 #             else:
    350 #                 fid.smoothing = 'No'
    351 
    352 #             fid.order = domain.default_order
    353 
    354 #             #Reference point
    355 #             #Start time in seconds since the epoch (midnight 1/1/1970)
    356 #             #FIXME: Use Georef
    357 #             fid.starttime = domain.starttime
    358 
    359 #             # dimension definitions
    360 #             fid.createDimension('number_of_volumes', self.number_of_volumes)
    361 #             fid.createDimension('number_of_vertices', 3)
    362 
    363 #             if domain.smooth is True:
    364 #                 #fid.createDimension('number_of_points', len(domain.vertexlist))
    365 #                 fid.createDimension('number_of_points', self.domain.number_of_full_nodes)
    366 
    367 #                 # FIXME(Ole): This will cause sww files for paralle domains to
    368 #                 # have ghost nodes stored (but not used by triangles).
    369 #                 # To clean this up, we have to change get_vertex_values and friends in
    370 #                 # quantity.py (but I can't be bothered right now)
    371 #             else:
    372 #                 fid.createDimension('number_of_points', 3*self.number_of_volumes)
    373 
    374 #             fid.createDimension('number_of_timesteps', None) #extensible
    375 
    376 #             # variable definitions
    377 #             fid.createVariable('x', self.precision, ('number_of_points',))
    378 #             fid.createVariable('y', self.precision, ('number_of_points',))
    379 #             fid.createVariable('elevation', self.precision, ('number_of_points',))
    380 #             if domain.geo_reference is not None:
    381 #                 domain.geo_reference.write_NetCDF(fid)
    382 
    383 #             #FIXME: Backwards compatibility
    384 #             fid.createVariable('z', self.precision, ('number_of_points',))
    385 #             #################################
    386 
    387 #             fid.createVariable('volumes', Int, ('number_of_volumes',
    388 #                                                 'number_of_vertices'))
    389 
    390 #             fid.createVariable('time', Float,  # Always use full precision lest two timesteps
    391 #                                              # close to each other may appear as the same step
    392 #                                ('number_of_timesteps',))
    393 
    394 #             fid.createVariable('stage', self.precision,
    395 #                                ('number_of_timesteps',
    396 #                                 'number_of_points'))
    397 
    398 #             fid.createVariable('xmomentum', self.precision,
    399 #                                ('number_of_timesteps',
    400 #                                 'number_of_points'))
    401 
    402 #             fid.createVariable('ymomentum', self.precision,
    403 #                                ('number_of_timesteps',
    404 #                                 'number_of_points'))
    405 
    406         #Close
     343                fid.texture = domain.texture
     344
     345            if domain.quantities_to_be_monitored is not None:
     346                fid.createDimension('singleton', 1)
     347                for q in domain.quantities_to_be_monitored:
     348                    #print 'doing', q
     349                    fid.createVariable(q+':extrema', self.precision,
     350                                       ('numbers_in_range',))
     351                    fid.createVariable(q+':min_location', self.precision,
     352                                       ('numbers_in_range',))
     353                    fid.createVariable(q+':max_location', self.precision,
     354                                       ('numbers_in_range',))
     355                    fid.createVariable(q+':min_time', self.precision,
     356                                       ('singleton',))
     357                    fid.createVariable(q+':max_time', self.precision,
     358                                       ('singleton',))
     359
     360                   
    407361        fid.close()
    408362
     
    438392        #
    439393        points = concatenate( (X[:,NewAxis],Y[:,NewAxis]), axis=1 )
    440         self.writer.triangulation(fid, points,
    441                                             V.astype(volumes.typecode()),
    442                                             Z,
    443                                             points_georeference= \
    444                                             domain.geo_reference)
    445         #
    446                                
    447 #         volumes[:] = V.astype(volumes.typecode())
    448 #         x[:] = X
    449 #         y[:] = Y
    450 #         z[:] = Z
    451 
    452 #         #FIXME: Backwards compatibility
    453 #         z = fid.variables['z']
    454 #         z[:] = Z
    455         ################################
    456 
    457 
    458 
    459         #Close
     394        self.writer.store_triangulation(fid,
     395                                        points,
     396                                        V.astype(volumes.typecode()),
     397                                        Z,
     398                                        points_georeference= \
     399                                        domain.geo_reference)
     400
     401        # Close
    460402        fid.close()
     403
    461404
    462405    def store_timestep(self, names):
     
    470413        from Numeric import choose
    471414       
    472         #Get NetCDF       
     415        # Get NetCDF       
    473416        retries = 0
    474417        file_open = False
    475418        while not file_open and retries < 10:
    476419            try:
    477                 fid = NetCDFFile(self.filename, 'a') #Open existing file
     420                fid = NetCDFFile(self.filename, 'a') # Open existing file
    478421            except IOError:
    479                 #This could happen if someone was reading the file.
    480                 #In that case, wait a while and try again
     422                # This could happen if someone was reading the file.
     423                # In that case, wait a while and try again
    481424                msg = 'Warning (store_timestep): File %s could not be opened'\
    482425                      %self.filename
     
    494437
    495438
    496         #Check to see if the file is already too big:
     439        # Check to see if the file is already too big:
    497440        time = fid.variables['time']
    498441        i = len(time)+1
     
    500443        file_size_increase =  file_size/i
    501444        if file_size + file_size_increase > self.max_size*(2**self.recursion):
    502             #in order to get the file name and start time correct,
    503             #I change the domian.filename and domain.starttime.
    504             #This is the only way to do this without changing
    505             #other modules (I think).
    506 
    507             #write a filename addon that won't break swollens reader
    508             #(10.sww is bad)
     445            # In order to get the file name and start time correct,
     446            # I change the domain.filename and domain.starttime.
     447            # This is the only way to do this without changing
     448            # other modules (I think).
     449
     450            # Write a filename addon that won't break swollens reader
     451            # (10.sww is bad)
    509452            filename_ext = '_time_%s'%self.domain.time
    510453            filename_ext = filename_ext.replace('.', '_')
    511             #remember the old filename, then give domain a
    512             #name with the extension
     454           
     455            # Remember the old filename, then give domain a
     456            # name with the extension
    513457            old_domain_filename = self.domain.get_name()
    514458            if not self.recursion:
     
    516460
    517461
    518             #change the domain starttime to the current time
     462            # Change the domain starttime to the current time
    519463            old_domain_starttime = self.domain.starttime
    520464            self.domain.starttime = self.domain.time
    521465
    522             #build a new data_structure.
     466            # Build a new data_structure.
    523467            next_data_structure=\
    524468                Data_format_sww(self.domain, mode=self.mode,\
     
    554498
    555499            if 'stage' in names and 'xmomentum' in names and \
    556                    'ymomentum' in names:
     500               'ymomentum' in names:
    557501
    558502                # Get stage
     
    573517                ymomentum, _ = Q.get_vertex_values(xy = False,
    574518                                          precision = self.precision)
    575                 self.writer.quantities(fid,
    576                                                  time=self.domain.time,
    577                                                  precision=self.precision,
    578                                                  stage=stage,
    579                                                  xmomentum=xmomentum,
    580                                                  ymomentum=ymomentum)
     519               
     520                # Write quantities to NetCDF
     521                self.writer.store_quantities(fid,
     522                                             time=self.domain.time,
     523                                             precision=self.precision,
     524                                             stage=stage,
     525                                             xmomentum=xmomentum,
     526                                             ymomentum=ymomentum)
    581527            else:
    582528                # This is producing a sww that is not standard.
    583                 #Store time
     529                # Store time
    584530                time[i] = self.domain.time
    585531               
     
    588534                    Q = domain.quantities[name]
    589535                    A,V = Q.get_vertex_values(xy = False,
    590                                           precision = self.precision)
    591 
    592                     #FIXME: Make this general (see below)
     536                                              precision = self.precision)
     537
     538                    # FIXME: Make this general (see below)
    593539                    if name == 'stage':
    594540                        z = fid.variables['elevation']
     
    603549                   #As in....
    604550                   #eval( name + '[i,:] = A.astype(self.precision)' )
    605                    #FIXME: But we need a UNIT test for that before
     551                   #FIXME (Ole): But we need a UNIT test for that before
    606552                   # refactoring
    607553
    608554
     555
     556            # Update extrema if requested
     557            domain = self.domain
     558            if domain.quantities_to_be_monitored is not None:
     559                for q, info in domain.quantities_to_be_monitored.items():
     560
     561                    if info['min'] is not None:
     562                        fid.variables[q + ':extrema'][0] = info['min']
     563                        fid.variables[q + ':min_location'][:] =\
     564                                        info['min_location']
     565                        fid.variables[q + ':min_time'][0] = info['min_time']
     566                       
     567                    if info['max'] is not None:
     568                        fid.variables[q + ':extrema'][1] = info['max']
     569                        fid.variables[q + ':max_location'][:] =\
     570                                        info['max_location']
     571                        fid.variables[q + ':max_time'][0] = info['max_time']
     572
     573           
    609574
    610575            #Flush and close
     
    614579
    615580
    616 #Class for handling checkpoints data
     581# Class for handling checkpoints data
    617582class Data_format_cpt(Data_format):
    618583    """Interface to native NetCDF format (.cpt)
     
    752717        Q = domain.quantities[name]
    753718        A,V = Q.get_vertex_values(xy=False,
    754                   precision = self.precision)
     719                                  precision = self.precision)
    755720
    756721        stage[i,:] = A.astype(self.precision)
     
    28632828
    28642829
    2865     #print number_of_latitudes, number_of_longitudes
     2830    # print number_of_latitudes, number_of_longitudes
    28662831    number_of_points = number_of_latitudes*number_of_longitudes
    28672832    number_of_volumes = (number_of_latitudes-1)*(number_of_longitudes-1)*2
     
    28832848                    basename_in + '_e.nc')
    28842849   
    2885     #Create new file
     2850    # Create new file
    28862851    starttime = times[0]
    28872852    sww = Write_sww()
    2888     sww.header(outfile, times, number_of_volumes,
    2889                          number_of_points, description=description,
    2890                          verbose=verbose)
    2891 
    2892     #Store
     2853    sww.store_header(outfile, times, number_of_volumes,
     2854                     number_of_points, description=description,
     2855                     verbose=verbose)
     2856
     2857    # Store
    28932858    from anuga.coordinate_transforms.redfearn import redfearn
    28942859    x = zeros(number_of_points, Float)  #Easting
     
    28972862
    28982863    if verbose: print 'Making triangular grid'
    2899     #Check zone boundaries
     2864
     2865    # Check zone boundaries
    29002866    refzone, _, _ = redfearn(latitudes[0],longitudes[0])
    29012867
     
    47274693    # work out sww_times and the index range this covers
    47284694    sww = Write_sww()
    4729     sww.header(outfile, times, len(volumes), len(points_utm),
    4730                          verbose=verbose)
     4695    sww.store_header(outfile, times, len(volumes), len(points_utm),
     4696                     verbose=verbose)
    47314697    outfile.mean_stage = mean_stage
    47324698    outfile.zscale = zscale
    47334699
    4734     sww.triangulation(outfile, points_utm, volumes,
     4700    sww.store_triangulation(outfile, points_utm, volumes,
    47354701                            elevation, zone,  new_origin=origin,
    47364702                            verbose=verbose)
     
    47454711            xmomentum = ua*h
    47464712            ymomentum = -1*va*h # -1 since in mux files south is positive.
    4747             sww.quantities(outfile,
    4748                                      slice_index=j - mux_times_start_i,
    4749                                      verbose=verbose,
    4750                                      stage=stage,
    4751                                      xmomentum=xmomentum,
    4752                                      ymomentum=ymomentum)
     4713            sww.store_quantities(outfile,
     4714                                 slice_index=j - mux_times_start_i,
     4715                                 verbose=verbose,
     4716                                 stage=stage,
     4717                                 xmomentum=xmomentum,
     4718                                 ymomentum=ymomentum)
    47534719        j += 1
    47544720    if verbose: sww.verbose_quantities(outfile)
     
    47784744
    47794745class Write_sww:
    4780     from anuga.shallow_water.shallow_water_domain import Domain
     4746    from anuga.shallow_water.shallow_water_domain import Domain
     4747
     4748    # FIXME (Ole): Hardwiring the conserved quantities like
     4749    # this could be a problem. I would prefer taking them from
     4750    # the instantiation of Domain.
    47814751    sww_quantities = Domain.conserved_quantities
     4752
     4753
    47824754    RANGE = '_range'
     4755    EXTREMA = ':extrema'
    47834756
    47844757    def __init__(self):
    47854758        pass
    47864759   
    4787     def header(self,outfile, times, number_of_volumes,
    4788                          number_of_points,
    4789                          description='Converted from XXX',
    4790                          smoothing=True,
    4791                          order=1, verbose=False):
     4760    def store_header(self,
     4761                     outfile,
     4762                     times,
     4763                     number_of_volumes,
     4764                     number_of_points,
     4765                     description='Converted from XXX',
     4766                     smoothing=True,
     4767                     order=1, verbose=False):
    47924768        """
    47934769        outfile - the name of the file that will be written
     
    48634839                               ('numbers_in_range',))
    48644840
     4841
    48654842        # Initialise ranges with small and large sentinels.
    48664843        # If this was in pure Python we could have used None sensibly
     
    48684845        outfile.variables[q+Write_sww.RANGE][1] = -max_float # Max
    48694846
    4870         #FIXME: Backwards compatibility
     4847        # FIXME: Backwards compatibility
    48714848        outfile.createVariable('z', precision, ('number_of_points',))
    48724849        #################################
     
    49004877
    49014878       
    4902     def triangulation(self,outfile, points_utm, volumes,
    4903                                 elevation, zone=None, new_origin=None,
    4904                                 points_georeference=None, verbose=False):
     4879    def store_triangulation(self,
     4880                            outfile,
     4881                            points_utm,
     4882                            volumes,
     4883                            elevation, zone=None, new_origin=None,
     4884                            points_georeference=None, verbose=False):
    49054885        """
    49064886       
     
    49944974        outfile.variables[q+Write_sww.RANGE][1] = max(elevation)
    49954975
    4996     def quantities(self, outfile, precision=Float,
    4997                              slice_index=None, time=None,
    4998                              verbose=False, **quant):
     4976
     4977    def store_quantities(self, outfile, precision=Float,
     4978                         slice_index=None, time=None,
     4979                         verbose=False, **quant):
    49994980        """
    50004981        Write the quantity info.
     
    50215002            file_time[slice_index] = time   
    50225003
    5023         # write the conserved quantities from Doamin.
     5004        # write the conserved quantities from Domain.
    50245005        # Typically stage,  xmomentum, ymomentum
    50255006        # other quantities will be ignored, silently.
     
    50495030                                       outfile.variables[q+Write_sww.RANGE][1])
    50505031        print '------------------------------------------------'
     5032
    50515033
    50525034       
  • anuga_core/source/anuga/shallow_water/test_data_manager.py

    r4701 r4704  
    256256        os.remove(sww.filename)
    257257
    258     def NOtest_sww_extrema(self):
     258    def test_sww_extrema(self):
    259259        """Test that extrema of quantities can be retrieved at every vertex
    260260        Extrema are updated at every *internal* timestep
     
    263263        domain = self.domain
    264264       
    265         domain.set_name('datatest' + str(id(self)))
     265        domain.set_name('extrema_test' + str(id(self)))
    266266        domain.format = 'sww'
    267267        domain.smooth = True
     
    271271        assert domain.monitor_time_interval is None       
    272272       
    273         domain.set_quantities_to_be_monitored(['stage', 'ymomentum'])
    274        
    275         assert len(domain.quantities_to_be_monitored) == 2
    276         assert domain.quantities_to_be_monitored[0] == 'stage'
    277         assert domain.quantities_to_be_monitored[1] == 'ymomentum'       
     273        domain.set_quantities_to_be_monitored(['xmomentum',
     274                                               'ymomentum',
     275                                               'stage-elevation'])
     276       
     277        assert len(domain.quantities_to_be_monitored) == 3
     278        assert domain.quantities_to_be_monitored.has_key('stage-elevation')
     279        assert domain.quantities_to_be_monitored.has_key('xmomentum')               
     280        assert domain.quantities_to_be_monitored.has_key('ymomentum')       
    278281        assert domain.monitor_polygon is None
    279282        assert domain.monitor_time_interval is None       
     
    282285
    283286        for t in domain.evolve(yieldstep = 1, finaltime = 1):
    284             print domain.timestepping_statistics()
    285             print domain.quantity_statistics()
     287            pass
     288            #print domain.timestepping_statistics()
     289            #print domain.quantity_statistics(precision = '%.8f')
    286290
    287291           
     
    290294
    291295        # Get the variables
    292         extrema = fid.variables['stage_extrema'][:]
    293         assert allclose(range, [])
    294 
    295         extrema = fid.variables['xmomentum_extrema'][:]
    296         assert allclose(range,[])
    297        
    298         extrema = fid.variables['ymomentum_extrema'][:]
    299         assert allclose(range,[])
     296        extrema = fid.variables['stage-elevation:extrema'][:]
     297        assert allclose(extrema, [0.00, 0.30])
     298
     299        loc = fid.variables['stage-elevation:min_location'][:]
     300        assert allclose(loc, [0.16666667, 0.33333333])
     301
     302        loc = fid.variables['stage-elevation:max_location'][:]       
     303        assert allclose(loc, [0.8333333, 0.16666667])       
     304
     305        time = fid.variables['stage-elevation:max_time'][:]
     306        assert allclose(time, 0.0)               
     307
     308        extrema = fid.variables['xmomentum:extrema'][:]
     309        assert allclose(extrema,[-0.06062178, 0.47886313])
     310       
     311        extrema = fid.variables['ymomentum:extrema'][:]
     312        assert allclose(extrema,[0.00, 0.06241221])       
    300313
    301314       
     
    67236736        number_of_points = len(points_utm)
    67246737        sww = Write_sww()
    6725         sww.header(outfile, times, number_of_volumes,
     6738        sww.store_header(outfile, times, number_of_volumes,
    67266739                         number_of_points, description='fully sick testing',
    6727                              verbose=self.verbose)
    6728         sww.triangulation(outfile, points_utm, volumes,
    6729                                     elevation,  new_origin=new_origin,
    6730                                     verbose=self.verbose)       
     6740                         verbose=self.verbose)
     6741        sww.store_triangulation(outfile, points_utm, volumes,
     6742                                elevation,  new_origin=new_origin,
     6743                                verbose=self.verbose)       
    67316744        outfile.close()
    67326745        fid = NetCDFFile(filename)
     
    67556768        number_of_points = len(points_utm)
    67566769        sww = Write_sww()
    6757         sww.header(outfile, times, number_of_volumes,
     6770        sww.store_header(outfile, times, number_of_volumes,
    67586771                         number_of_points, description='fully sick testing',
    67596772                         verbose=self.verbose)
    6760         sww.triangulation(outfile, points_utm, volumes,
    6761                                     elevation,  new_origin=new_origin,
    6762                                     verbose=self.verbose)       
     6773        sww.store_triangulation(outfile, points_utm, volumes,
     6774                                elevation,  new_origin=new_origin,
     6775                                verbose=self.verbose)       
    67636776        outfile.close()
    67646777        fid = NetCDFFile(filename)
     
    67916804        number_of_points = len(points_utm)
    67926805        sww = Write_sww()
    6793         sww.header(outfile, times, number_of_volumes,
     6806        sww.store_header(outfile, times, number_of_volumes,
    67946807                         number_of_points, description='fully sick testing',
    67956808                         verbose=self.verbose)
    6796         sww.triangulation(outfile, points_utm, volumes,
    6797                                     elevation,  new_origin=new_origin,
    6798                                     verbose=self.verbose)       
     6809        sww.store_triangulation(outfile, points_utm, volumes,
     6810                                elevation,  new_origin=new_origin,
     6811                                verbose=self.verbose)       
    67996812        outfile.close()
    68006813        fid = NetCDFFile(filename)
     
    68306843        number_of_points = len(points_utm)
    68316844        sww = Write_sww()
    6832         sww.header(outfile, times, number_of_volumes,
     6845        sww.store_header(outfile, times, number_of_volumes,
    68336846                         number_of_points, description='fully sick testing',
    68346847                         verbose=self.verbose)
    6835         sww.triangulation(outfile, points_utm, volumes,
    6836                                     elevation,  new_origin=new_origin,
    6837                                     points_georeference=points_georeference,
    6838                                     verbose=self.verbose)       
     6848        sww.store_triangulation(outfile, points_utm, volumes,
     6849                                elevation,  new_origin=new_origin,
     6850                                points_georeference=points_georeference,
     6851                                verbose=self.verbose)       
    68396852        outfile.close()
    68406853        fid = NetCDFFile(filename)
     
    68666879        number_of_points = len(points_utm)
    68676880        sww = Write_sww()
    6868         sww.header(outfile, times, number_of_volumes,
     6881        sww.store_header(outfile, times, number_of_volumes,
    68696882                         number_of_points, description='fully sick testing',
    68706883                         verbose=self.verbose)
    6871         sww.triangulation(outfile, points_utm, volumes,
    6872                                     elevation,  new_origin=new_origin,
    6873                                     points_georeference=points_georeference,
    6874                                     verbose=self.verbose)       
     6884        sww.store_triangulation(outfile, points_utm, volumes,
     6885                                elevation,  new_origin=new_origin,
     6886                                points_georeference=points_georeference,
     6887                                verbose=self.verbose)       
    68756888        outfile.close()
    68766889        fid = NetCDFFile(filename)
     
    72817294    #suite = unittest.makeSuite(Test_Data_Manager,'test_sww_header')
    72827295    #suite = unittest.makeSuite(Test_Data_Manager,'test_export_grid_parallel')
    7283     #suite = unittest.makeSuite(Test_Data_Manager,'t')
    72847296    suite = unittest.makeSuite(Test_Data_Manager,'test')
     7297    #suite = unittest.makeSuite(Test_Data_Manager,'test_sww_extrema')
    72857298
    72867299   
  • anuga_core/source/anuga/shallow_water/test_shallow_water_domain.py

    r4686 r4704  
    46354635
    46364636
     4637
     4638    def test_extrema(self):
     4639        """Test that extrema of quantities are computed correctly
     4640        Extrema are updated at every *internal* timestep
     4641        """
     4642
     4643        from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular_cross
     4644
     4645        initial_runup_height = -0.4
     4646        final_runup_height = -0.3
     4647
     4648
     4649        #--------------------------------------------------------------
     4650        # Setup computational domain
     4651        #--------------------------------------------------------------
     4652        N = 5
     4653        points, vertices, boundary = rectangular_cross(N, N)
     4654        domain = Domain(points, vertices, boundary)
     4655        domain.set_name('extrema_test')
     4656
     4657        #--------------------------------------------------------------
     4658        # Setup initial conditions
     4659        #--------------------------------------------------------------
     4660        def topography(x,y):
     4661            return -x/2                             # linear bed slope
     4662           
     4663
     4664        domain.set_quantity('elevation', topography)       # Use function for elevation
     4665        domain.set_quantity('friction', 0.)                # Zero friction
     4666        domain.set_quantity('stage', initial_runup_height) # Constant negative initial stage
     4667        domain.set_quantities_to_be_monitored(['stage', 'stage-elevation'],
     4668                                              time_interval = [0.5, 2.7],
     4669                                              polygon = [[0,0], [0,1], [1,1], [1,0]])
     4670       
     4671        assert len(domain.quantities_to_be_monitored) == 2
     4672        assert domain.quantities_to_be_monitored.has_key('stage')
     4673        assert domain.quantities_to_be_monitored.has_key('stage-elevation')
     4674        for key in domain.quantities_to_be_monitored['stage'].keys():
     4675            assert domain.quantities_to_be_monitored['stage'][key] is None       
     4676
     4677
     4678        #--------------------------------------------------------------
     4679        # Setup boundary conditions
     4680        #--------------------------------------------------------------
     4681        Br = Reflective_boundary(domain)              # Reflective wall
     4682        Bd = Dirichlet_boundary([final_runup_height,  # Constant inflow
     4683                                 0,
     4684                                 0])
     4685
     4686        # All reflective to begin with (still water)
     4687        domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
     4688
     4689
     4690        #--------------------------------------------------------------
     4691        # Let triangles adjust and check extrema
     4692        #--------------------------------------------------------------
     4693        for t in domain.evolve(yieldstep = 0.1, finaltime = 1.0):
     4694            domain.quantity_statistics() # Run it silently
     4695
     4696
     4697
     4698        #--------------------------------------------------------------
     4699        # Test extrema
     4700        #--------------------------------------------------------------
     4701
     4702        stage = domain.quantities_to_be_monitored['stage']
     4703        assert stage['min'] <= stage['max']
     4704
     4705        #print stage['min'], stage['max']
     4706        assert allclose(stage['min'], initial_runup_height,
     4707                        rtol = 1.0/N) # First order accuracy
     4708
     4709
     4710        depth = domain.quantities_to_be_monitored['stage-elevation']
     4711        assert depth['min'] <= depth['max']
     4712        assert depth['min'] >= 0.0
     4713        assert depth['max'] >= 0.0       
     4714        ##assert depth[1] <= ?? initial_runup_height       
     4715
     4716
     4717        #--------------------------------------------------------------
     4718        # Update boundary to allow inflow
     4719        #--------------------------------------------------------------
     4720        domain.set_boundary({'right': Bd})
     4721
     4722       
     4723        #--------------------------------------------------------------
     4724        # Evolve system through time
     4725        #--------------------------------------------------------------
     4726        for t in domain.evolve(yieldstep = 0.1, finaltime = 3.0):
     4727            #domain.write_time()
     4728            domain.quantity_statistics() # Run it silently           
     4729           
     4730   
     4731        #--------------------------------------------------------------
     4732        # Test extrema again
     4733        #--------------------------------------------------------------
     4734
     4735        stage = domain.quantities_to_be_monitored['stage']
     4736        assert stage['min'] <= stage['max']
     4737
     4738        assert allclose(stage['min'], initial_runup_height,
     4739                        rtol = 1.0/N) # First order accuracy       
     4740
     4741        depth = domain.quantities_to_be_monitored['stage-elevation']
     4742        assert depth['min'] <= depth['max']
     4743        assert depth['min'] >= 0.0
     4744        assert depth['max'] >= 0.0       
     4745       
     4746
     4747
    46374748    def test_tight_slope_limiters(self):
    46384749        """Test that new slope limiters (Feb 2007) don't induce extremely
     
    48724983if __name__ == "__main__":
    48734984
    4874     suite = unittest.makeSuite(Test_Shallow_Water,'test')   
     4985    #suite = unittest.makeSuite(Test_Shallow_Water,'test')   
     4986    suite = unittest.makeSuite(Test_Shallow_Water,'test_extrema')   
    48754987    #suite = unittest.makeSuite(Test_Shallow_Water,'test_tight_slope_limiters')
    48764988    #suite = unittest.makeSuite(Test_Shallow_Water,'test_get_maximum_inundation_from_sww')
Note: See TracChangeset for help on using the changeset viewer.