Changeset 7342


Ignore:
Timestamp:
Aug 9, 2009, 1:20:42 AM (16 years ago)
Author:
ole
Message:

Refactored stofrage of quantities to use new concept of static and dynamic quantities.
This is in preparation for flexibly allowing quantities such as elevation or friction
to be time dependent.

All tests and the okushiri validation pass.

Location:
anuga_core
Files:
9 edited

Legend:

Unmodified
Added
Removed

  • anuga_core/documentation/user_manual/anuga_user_manual.tex

    r7313 r7342  
    386386In addition, the following statement could be used to state that
    387387quantities \code{stage}, \code{xmomentum} and \code{ymomentum} are
    388 to be stored:
    389 
    390 \begin{verbatim}
    391 domain.set_quantities_to_be_stored(['stage', 'xmomentum', 'ymomentum'])
    392 \end{verbatim}
    393 
    394 However, this is not necessary, as the above quantities are always stored by default.
     388to be stored at every timestep and \code{elevation} only once at
     389the beginning of the simulation:
     390
     391\begin{verbatim}
     392domain.set_quantities_to_be_stored({'stage': 2, 'xmomentum': 2, 'ymomentum': 2, 'elevation': 1})
     393\end{verbatim}
     394
     395However, this is not necessary, as the above is the default behaviour.
    395396
    396397\subsection{Initial Conditions}
    … …  
    19591960data from the SWW files with that of the model internally.
    19601961\end{methoddesc}
     1962
     1963
     1964\begin{methoddesc}{\emph{<domain>}.set_quantities_to_be_stored}{quantity_dictionary}
     1965Module: \module{shallow_water.shallow_water_domain}
     1966
     1967Selects quantities that is to be stored in the sww files.
     1968The argument can be None, in which case nothing is stored.
     1969
     1970Otherwise, the argument must be a dictionary where the keys are names of quantities
     1971already defined within ANUGA and the values are either 1 or 2. If the value is 1, the quantity
     1972will be stored once at the beginning of the simulation, if the value is 2 it will be stored
     1973at each timestep. The ANUGA default is equivalent to the call
     1974\begin{verbatim}
     1975domain.set_quantities_to_be_stored({'elevation': 1,
     1976                                    'stage': 2,
     1977                                    'xmomentum': 2,
     1978                                    'ymomentum': 2})
     1979\end{verbatim}
     1980\end{methoddesc}
     1981
    19611982
    19621983% Structural methods

  • anuga_core/source/anuga/abstract_2d_finite_volumes/domain.py

    r7317 r7342  
    162162        self.quantities = {}
    163163
    164         # FIXME: remove later - maybe OK, though....
    165164        for name in self.conserved_quantities:
    166165            self.quantities[name] = Quantity(self)

  • anuga_core/source/anuga/abstract_2d_finite_volumes/test_util.py

    r7340 r7342  
    141141        domain1.set_datadir('.')
    142142        domain1.set_name('spatio_temporal_boundary_source_%d' %(id(self)))
    143         domain1.quantities_to_be_stored = ['stage', 'xmomentum', 'ymomentum']
    144143
    145144        #Bed-slope, friction and IC at vertices (and interpolated elsewhere)
    … …  
    359358        domain1.set_datadir('.')
    360359        domain1.set_name('spatio_temporal_boundary_source_%d' %(id(self)))
    361         domain1.quantities_to_be_stored = ['stage', 'xmomentum', 'ymomentum']
    362360
    363361        #Bed-slope, friction and IC at vertices (and interpolated elsewhere)
    … …  
    820818            #Store and advance time
    821819            domain.time = t
    822             domain.store_timestep(domain.conserved_quantities)
     820            domain.store_timestep()
    823821            t += dt
    824822
    … …  
    16181616        sww = SWW_file(domain)
    16191617        sww.store_connectivity()
    1620         sww.store_timestep(['stage', 'xmomentum', 'ymomentum','elevation'])
     1618        sww.store_timestep()
    16211619        domain.set_quantity('stage', 10.0) # This is automatically limited
    16221620        # so it will not be less than the elevation
    16231621        domain.time = 2.
    1624         sww.store_timestep(['stage','elevation', 'xmomentum', 'ymomentum'])
     1622        sww.store_timestep()
    16251623
    16261624        # test the function
    … …  
    17381736        sww = SWW_file(domain)
    17391737        sww.store_connectivity()
    1740         sww.store_timestep(['stage', 'xmomentum', 'ymomentum'])
     1738        sww.store_timestep()
    17411739        domain.set_quantity('stage', 10.0) # This is automatically limited
    17421740        # so it will not be less than the elevation
    17431741        domain.time = 2.
    1744         sww.store_timestep(['stage', 'xmomentum', 'ymomentum'])
     1742        sww.store_timestep()
    17451743
    17461744        # test the function
    … …  
    18571855        sww = SWW_file(domain)
    18581856        sww.store_connectivity()
    1859         sww.store_timestep(['stage', 'xmomentum', 'ymomentum','elevation'])
     1857        sww.store_timestep()
    18601858        domain.set_quantity('stage', 10.0) # This is automatically limited
    18611859        # so it will not be less than the elevation
    18621860        domain.time = 2.
    1863         sww.store_timestep(['stage','elevation', 'xmomentum', 'ymomentum'])
     1861        sww.store_timestep()
    18641862
    18651863        # test the function

  • anuga_core/source/anuga/damage_modelling/test_inundation_damage.py

    r7340 r7342  
    114114        sww = SWW_file(self.domain)
    115115        sww.store_connectivity()
    116         sww.store_timestep(['stage', 'xmomentum', 'ymomentum'])
     116        sww.store_timestep()
    117117        self.domain.time = 2.
    118         sww.store_timestep(['stage', 'xmomentum', 'ymomentum'])
     118        sww.store_timestep()
    119119        self.sww = sww # so it can be deleted
    120120       
    … …  
    184184        sww = SWW_file(domain)
    185185        sww.store_connectivity()
    186         sww.store_timestep(['stage', 'xmomentum', 'ymomentum'])
     186        sww.store_timestep()
    187187        domain.time = 2.
    188         sww.store_timestep(['stage', 'xmomentum', 'ymomentum'])
     188        sww.store_timestep()
    189189        self.swwII = sww # so it can be deleted
    190190
    … …  
    284284        sww = SWW_file(domain)
    285285        sww.store_connectivity()
    286         sww.store_timestep(['stage', 'xmomentum', 'ymomentum'])
     286        sww.store_timestep()
    287287        domain.set_quantity('stage', -0.3)
    288288        domain.time = 2.
    289         sww.store_timestep(['stage', 'xmomentum', 'ymomentum'])
     289        sww.store_timestep()
    290290       
    291291        #Create a csv file
    … …  
    356356        sww = SWW_file(domain)
    357357        sww.store_connectivity()
    358         sww.store_timestep(['stage', 'xmomentum', 'ymomentum'])
     358        sww.store_timestep()
    359359        domain.set_quantity('stage', -0.3)
    360360        domain.time = 2.
    361         sww.store_timestep(['stage', 'xmomentum', 'ymomentum'])
     361        sww.store_timestep()
    362362       
    363363        #Create a csv file

  • anuga_core/source/anuga/fit_interpolate/test_interpolate.py

    r7340 r7342  
    17151715        sww = SWW_file(domain)
    17161716        sww.store_connectivity()
    1717         sww.store_timestep(['stage', 'xmomentum', 'ymomentum'])
     1717        sww.store_timestep()
    17181718        domain.set_quantity('stage', 10.0) # This is automatically limited
    17191719        # So it will not be less than the elevation
    17201720        domain.time = 2.
    1721         sww.store_timestep(['stage', 'xmomentum', 'ymomentum'])
     1721        sww.store_timestep()
    17221722
    17231723        # Test the function

  • anuga_core/source/anuga/shallow_water/data_manager.py

    r7340 r7342  
    347347    # @param recursion ??
    348348    # @note Prepare the underlying data file if mode starts with 'w'.
    349     def __init__(self, domain, mode=netcdf_mode_w, max_size=2000000000, recursion=False):
     349    def __init__(self, domain,
     350                 mode=netcdf_mode_w, max_size=2000000000, recursion=False):
    350351        from Scientific.IO.NetCDF import NetCDFFile
    351352
    … …  
    362363            self.minimum_storable_height = default_minimum_storable_height
    363364
    364         # call owning constructor
     365        # Call parent constructor
    365366        Data_format.__init__(self, domain, 'sww', mode)
    366367
     368        # Get static and dynamic quantities from domain
     369        static_quantities = []
     370        dynamic_quantities = []
     371       
     372        for q in domain.quantities_to_be_stored:
     373            flag = domain.quantities_to_be_stored[q]
     374       
     375            msg = 'Quantity %s is requested to be stored ' % q
     376            msg += 'but it does not exist in domain.quantities'
     377            assert q in domain.quantities, msg
     378       
     379            assert flag in [1,2]
     380            if flag == 1: static_quantities.append(q)
     381            if flag == 2: dynamic_quantities.append(q)               
     382                       
     383       
    367384        # NetCDF file definition
    368385        fid = NetCDFFile(self.filename, mode)
    … …  
    370387            description = 'Output from anuga.abstract_2d_finite_volumes ' \
    371388                          'suitable for plotting'
    372             self.writer = Write_sww(['elevation'], domain.conserved_quantities)
     389                         
     390            self.writer = Write_sww(static_quantities, dynamic_quantities)
    373391            self.writer.store_header(fid,
    374392                                     domain.starttime,
    … …  
    421439    # @brief Store connectivity data into the underlying data file.
    422440    def store_connectivity(self):
    423         """Specialisation of store_connectivity for net CDF format
    424 
    425         Writes x,y,z coordinates of triangles constituting
    426         the bed elevation.
     441        """Store information about nodes, triangles and static quantities
     442
     443        Writes x,y coordinates of triangles and their connectivity.
     444       
     445        Store also any quantity that has been identified as static.
    427446        """
    428447
     448        # FIXME: Change name to reflect the fact thta this function
     449        # stores both connectivity (triangulation) and static quantities
     450       
    429451        from Scientific.IO.NetCDF import NetCDFFile
    430452
    … …  
    434456        fid = NetCDFFile(self.filename, netcdf_mode_a)
    435457
    436         # Get X, Y and bed elevation Z
    437         Q = domain.quantities['elevation']
    438         X,Y,Z,V = Q.get_vertex_values(xy=True, precision=self.precision)
     458        # Get X, Y from one (any) of the quantities
     459        Q = domain.quantities.values()[0]
     460        X,Y,_,V = Q.get_vertex_values(xy=True, precision=self.precision)
    439461
    440462        # store the connectivity data
    441         points = num.concatenate( (X[:,num.newaxis],Y[:,num.newaxis]), axis=1 )
     463        points = num.concatenate((X[:,num.newaxis],Y[:,num.newaxis]), axis=1)
    442464        self.writer.store_triangulation(fid,
    443465                                        points,
    444466                                        V.astype(num.float32),
    445                                         Z,
    446467                                        points_georeference=\
    447468                                            domain.geo_reference)
    448469
     470
     471        # Get names of static quantities
     472        static_quantities = {}
     473        for name in self.writer.static_quantities:
     474            Q = domain.quantities[name]
     475            A, _ = Q.get_vertex_values(xy=False,
     476                                       precision=self.precision)
     477            static_quantities[name] = A
     478       
     479        # Store static quantities       
     480        self.writer.store_static_quantities(fid, **static_quantities)
     481                                           
    449482        fid.close()
    450483
    … …  
    452485    # @brief Store time and time dependent quantities
    453486    # to the underlying data file.
    454     def store_timestep(self, names=None):
     487    def store_timestep(self):
    455488        """Store time and time dependent quantities
    456489        """
    … …  
    460493        from time import sleep
    461494        from os import stat
    462 
    463         # Get names of quantities to be stored every time step
    464         names = self.writer.dynamic_quantities   
    465495
    466496        # Get NetCDF
    … …  
    469499        while not file_open and retries < 10:
    470500            try:
    471                 fid = NetCDFFile(self.filename, netcdf_mode_a) # Open existing file
     501                # Open existing file
     502                fid = NetCDFFile(self.filename, netcdf_mode_a)
    472503            except IOError:
    473504                # This could happen if someone was reading the file.
    … …  
    497528            # other modules (I think).
    498529
    499             # Write a filename addon that won't break swollens reader
     530            # Write a filename addon that won't break the anuga viewers
    500531            # (10.sww is bad)
    501532            filename_ext = '_time_%s' % self.domain.time
    … …  
    508539                self.domain.set_name(old_domain_filename + filename_ext)
    509540
    510             # Change the domain starttime to the current time
     541            # Temporarily change the domain starttime to the current time
    511542            old_domain_starttime = self.domain.starttime
    512             self.domain.starttime = self.domain.time
     543            self.domain.starttime = self.domain.get_time()
    513544
    514545            # Build a new data_structure.
    515546            next_data_structure = SWW_file(self.domain, mode=self.mode,
    516                                                   max_size=self.max_size,
    517                                                   recursion=self.recursion+1)
     547                                           max_size=self.max_size,
     548                                           recursion=self.recursion+1)
    518549            if not self.recursion:
    519550                log.critical('    file_size = %s' % file_size)
    … …  
    521552                             % next_data_structure.filename)
    522553
    523             #set up the new data_structure
     554            # Set up the new data_structure
    524555            self.domain.writer = next_data_structure
    525556
    526             #FIXME - could be cleaner to use domain.store_timestep etc.
     557            # Store connectivity and first timestep
    527558            next_data_structure.store_connectivity()
    528             next_data_structure.store_timestep(names)
     559            next_data_structure.store_timestep()
    529560            fid.sync()
    530561            fid.close()
    531562
    532             #restore the old starttime and filename
     563            # Restore the old starttime and filename
    533564            self.domain.starttime = old_domain_starttime
    534565            self.domain.set_name(old_domain_filename)
    … …  
    539570            # Get the variables
    540571            time = fid.variables['time']
    541             stage = fid.variables['stage']
    542             xmomentum = fid.variables['xmomentum']
    543             ymomentum = fid.variables['ymomentum']
    544572            i = len(time)
    545             if type(names) not in [types.ListType, types.TupleType]:
    546                 names = [names]
    547 
    548             if 'stage' in names \
    549                and 'xmomentum' in names \
    550                and 'ymomentum' in names:
    551                 # Get stage, elevation, depth and select only those
    552                 # values where minimum_storable_height is exceeded
     573             
     574            if 'stage' in self.writer.dynamic_quantities:           
     575                # Select only those values for stage,
     576                # xmomentum and ymomentum (if stored) where
     577                # depth exceeds minimum_storable_height
     578                #
     579                # In this branch it is assumed that elevation
     580                # is also available as a quantity           
     581           
    553582                Q = domain.quantities['stage']
    554                 A, _ = Q.get_vertex_values(xy=False, precision=self.precision)
    555                 z = fid.variables['elevation']
    556 
    557                 storable_indices = (A-z[:] >= self.minimum_storable_height)
    558                 stage = num.choose(storable_indices, (z[:], A))
    559 
    560                 # Define a zero vector of same size and type as A
    561                 # for use with momenta
    562                 null = num.zeros(num.size(A), A.dtype.char)
    563 
    564                 # Get xmomentum where depth exceeds minimum_storable_height
    565                 Q = domain.quantities['xmomentum']
    566                 xmom, _ = Q.get_vertex_values(xy=False,
    567                                               precision=self.precision)
    568                 xmomentum = num.choose(storable_indices, (null, xmom))
    569 
    570 
    571                 # Get ymomentum where depth exceeds minimum_storable_height
    572                 Q = domain.quantities['ymomentum']
    573                 ymom, _ = Q.get_vertex_values(xy=False,
    574                                               precision=self.precision)
    575                 ymomentum = num.choose(storable_indices, (null, ymom))
    576 
    577                 # Write quantities to underlying data  file
    578                 self.writer.store_quantities(fid,
    579                                              time=self.domain.time,
    580                                              sww_precision=self.precision,
    581                                              stage=stage,
    582                                              xmomentum=xmomentum,
    583                                              ymomentum=ymomentum)
     583                w, _ = Q.get_vertex_values(xy=False)
     584               
     585                Q = domain.quantities['elevation']
     586                z, _ = Q.get_vertex_values(xy=False)               
     587               
     588                storable_indices = (w-z >= self.minimum_storable_height)
    584589            else:
    585                 msg = 'Quantities stored must be: stage, xmomentum, ymomentum, '
    586                 msg += 'but I got: ' + str(names)
    587                 raise Exception, msg
     590                # Very unlikely branch
     591                storable_indices = None # This means take all
     592           
     593           
     594            # Now store dynamic quantities
     595            dynamic_quantities = {}
     596            for name in self.writer.dynamic_quantities:
     597                netcdf_array = fid.variables[name]
     598               
     599                Q = domain.quantities[name]
     600                A, _ = Q.get_vertex_values(xy=False,
     601                                           precision=self.precision)
     602
     603                if storable_indices is not None:
     604                    if name == 'stage':
     605                        A = num.choose(storable_indices, (z, A))
     606
     607                    if name in ['xmomentum', 'ymomentum']:
     608                        # Get xmomentum where depth exceeds
     609                        # minimum_storable_height
     610                       
     611                        # Define a zero vector of same size and type as A
     612                        # for use with momenta
     613                        null = num.zeros(num.size(A), A.dtype.char)
     614                        A = num.choose(storable_indices, (null, A))
     615               
     616                dynamic_quantities[name] = A
     617               
     618                                       
     619            # Store dynamic quantities
     620            self.writer.store_quantities(fid,
     621                                         time=self.domain.time,
     622                                         sww_precision=self.precision,
     623                                         **dynamic_quantities)
     624
    588625
    589626            # Update extrema if requested
    … …  
    608645
    609646
     647
    610648##
    611649# @brief Class for handling checkpoints data
    … …  
    630668        fid = NetCDFFile(self.filename, mode)
    631669        if mode[0] == 'w':
    632             #Create new file
     670            # Create new file
    633671            fid.institution = 'Geoscience Australia'
    634672            fid.description = 'Checkpoint data'
    … …  
    636674            fid.order = domain.default_order
    637675
    638             # dimension definitions
     676            # Dimension definitions
    639677            fid.createDimension('number_of_volumes', self.number_of_volumes)
    640678            fid.createDimension('number_of_vertices', 3)
    641679
    642             #Store info at all vertices (no smoothing)
     680            # Store info at all vertices (no smoothing)
    643681            fid.createDimension('number_of_points', 3*self.number_of_volumes)
    644682            fid.createDimension('number_of_timesteps', None) #extensible
    645683
    646             # variable definitions
    647 
    648             #Mesh
     684            # Variable definitions
     685
     686            # Mesh
    649687            fid.createVariable('x', self.precision, ('number_of_points',))
    650688            fid.createVariable('y', self.precision, ('number_of_points',))
    … …  
    53225360
    53235361    sww.store_triangulation(outfile, points_utm, volumes,
    5324                             elevation, zone,  new_origin=origin,
     5362                            zone, 
     5363                            new_origin=origin,
    53255364                            verbose=verbose)
     5365    sww.store_static_quantities(outfile, elevation=elevation)
    53265366
    53275367    if verbose: log.critical('Converting quantities')
    … …  
    58715911                     number_of_volumes,
    58725912                     number_of_points,
    5873                      description='Converted from XXX',
     5913                     description='Generated by ANUGA',
    58745914                     smoothing=True,
    58755915                     order=1,
    … …  
    59185958            number_of_times = 0
    59195959            starttime = times
    5920             #times = ensure_numeric([])
     5960
    59215961
    59225962        outfile.starttime = starttime
    … …  
    59926032    # @param points_utm List or array of points in UTM.
    59936033    # @param volumes
    5994     # @param elevation
    59956034    # @param zone
    59966035    # @param new_origin georeference that the points can be set to.
    … …  
    60016040                            points_utm,
    60026041                            volumes,
    6003                             elevation, zone=None, new_origin=None,
    6004                             points_georeference=None, verbose=False):
     6042                            zone=None,
     6043                            new_origin=None,
     6044                            points_georeference=None,
     6045                            verbose=False):
    60056046        """
    60066047        new_origin - qa georeference that the points can be set to. (Maybe
    … …  
    60366077        points_utm = num.array(points_utm)
    60376078
    6038         # given the two geo_refs and the points, do the stuff
     6079        # Given the two geo_refs and the points, do the stuff
    60396080        # described in the method header
    60406081        # if this is needed else where, pull out as a function
    … …  
    60616102        x =  points[:,0]
    60626103        y =  points[:,1]
    6063         z = outfile.variables['z'][:]
    60646104
    60656105        if verbose:
    … …  
    60716111            log.critical('    y in [%f, %f], len(lon) == %d'
    60726112                         % (min(y), max(y), len(y)))
    6073             log.critical('    z in [%f, %f], len(z) == %d'
    6074                          % (min(elevation), max(elevation), len(elevation)))
     6113            #log.critical('    z in [%f, %f], len(z) == %d'
     6114            #             % (min(elevation), max(elevation), len(elevation)))
    60756115            log.critical('geo_ref: %s' % str(geo_ref))
    60766116            log.critical('------------------------------------------------')
    … …  
    60786118        outfile.variables['x'][:] = points[:,0] #- geo_ref.get_xllcorner()
    60796119        outfile.variables['y'][:] = points[:,1] #- geo_ref.get_yllcorner()
    6080         outfile.variables['z'][:] = elevation
    6081         outfile.variables['elevation'][:] = elevation  #FIXME HACK
    60826120        outfile.variables['volumes'][:] = volumes.astype(num.int32) #On Opteron 64
    60836121
    6084         q = 'elevation'
    6085         # This updates the _range values
    6086         outfile.variables[q + Write_sww.RANGE][0] = num.min(elevation)
    6087         outfile.variables[q + Write_sww.RANGE][1] = num.max(elevation)
    6088 
    6089 
     6122
     6123
     6124    # @brief Write the static quantity data to the underlying file.
     6125    # @param outfile Handle to open underlying file.
     6126    # @param sww_precision Format of quantity data to write (default Float32).
     6127    # @param verbose True if this function is to be verbose.
     6128    # @param **quant
     6129    def store_static_quantities(self,
     6130                                outfile,
     6131                                sww_precision=num.float32,
     6132                                verbose=False,
     6133                                **quant):
     6134        """
     6135        Write the static quantity info.
     6136
     6137        **quant is extra keyword arguments passed in. These must be
     6138          the numpy arrays to be stored in the sww file at each timestep.
     6139
     6140        The argument sww_precision allows for storing as either
     6141        * single precision (default): num.float32
     6142        * double precision: num.float64 or num.float
     6143
     6144        Precondition:
     6145            store_triangulation and
     6146            store_header have been called.
     6147        """
     6148
     6149        # The dictionary quant must contain numpy arrays for each name.
     6150        # These will typically be quantities from Domain such as friction
     6151        #
     6152        # Arrays not listed in static_quantitiues will be ignored, silently.
     6153        #
     6154        # This method will also write the ranges for each quantity,
     6155        # e.g. stage_range, xmomentum_range and ymomentum_range
     6156        for q in self.static_quantities:
     6157            if not quant.has_key(q):
     6158                msg = 'Values for quantity %s was not specified in ' % q
     6159                msg += 'store_quantities so they cannot be stored.'
     6160                raise NewQuantity, msg
     6161            else:
     6162                q_values = ensure_numeric(quant[q])
     6163               
     6164                x = q_values.astype(sww_precision)
     6165                outfile.variables[q][:] = x
     6166       
     6167                # This populates the _range values
     6168                outfile.variables[q + Write_sww.RANGE][0] = num.min(x)
     6169                outfile.variables[q + Write_sww.RANGE][1] = num.max(x)               
     6170                   
     6171        outfile.variables['z'][:] = outfile.variables['elevation'][:] #FIXME HACK
     6172
     6173                   
     6174                   
     6175       
     6176       
    60906177    ##
    60916178    # @brief Write the quantity data to the underlying file.
    … …  
    61486235                raise NewQuantity, msg
    61496236            else:
    6150                 q_values = quant[q]
     6237                q_values = ensure_numeric(quant[q])
    61516238               
    61526239                x = q_values.astype(sww_precision)

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

    r7340 r7342  
    120120import anuga.utilities.log as log
    121121
     122import types
    122123from types import IntType, FloatType
    123124from warnings import warn
    … …  
    132133class Domain(Generic_Domain):
    133134
    134     conserved_quantities = ['stage', 'xmomentum', 'ymomentum']
    135     other_quantities = ['elevation', 'friction']
     135    #conserved_quantities = ['stage', 'xmomentum', 'ymomentum']
     136    #other_quantities = ['elevation', 'friction']
    136137
    137138    ##
    … …  
    169170                 number_of_full_triangles=None):
    170171
     172        # Define quantities for the shallow_water domain         
     173        conserved_quantities = ['stage', 'xmomentum', 'ymomentum']         
    171174        other_quantities = ['elevation', 'friction']
     175       
    172176        Generic_Domain.__init__(self,
    173177                                coordinates,
    174178                                vertices,
    175179                                boundary,
    176                                 Domain.conserved_quantities,
    177                                 Domain.other_quantities,
     180                                conserved_quantities,
     181                                other_quantities,
    178182                                tagged_elements,
    179183                                geo_reference,
    … …  
    209213        # Stored output
    210214        self.store = True
    211         #self.format = 'sww'
    212215        self.set_store_vertices_uniquely(False)
    213216        self.minimum_storable_height = minimum_storable_height
    214         self.quantities_to_be_stored = ['stage', 'xmomentum', 'ymomentum']
     217        self.quantities_to_be_stored = {'elevation': 1,
     218                                        'stage': 2,
     219                                        'xmomentum': 2,
     220                                        'ymomentum': 2}
    215221
    216222        # Limiters
    … …  
    317323        self.points_file_block_line_size = points_file_block_line_size
    318324
     325       
     326    # FIXME: Probably obsolete in its curren form   
    319327    ##
    320328    # @brief Set the quantities that will be written to an SWW file.
    … …  
    323331    # @note If 'q' is None, no quantities will be stored in the SWW file.
    324332    def set_quantities_to_be_stored(self, q):
    325         """Specify which quantities will be stored in the sww file.
    326 
     333        """Specify which quantities will be stored in the sww file
     334       
    327335        q must be either:
    328           - the name of a quantity
    329           - a list of quantity names
     336          - a dictionary with quantity names
     337          - a list of quantity names (for backwards compatibility)
    330338          - None
    331339
    332         In the two first cases, the named quantities will be stored at
    333         each yieldstep (This is in addition to the quantities elevation
    334         and friction)
    335 
     340        The format of the dictionary is as follows
     341       
     342        quantity_name: flag where flag must be either 1 or 2.
     343        If flag is 1, the quantity is considered static and will be
     344        stored once at the beginning of the simulation in a 1D array.
     345       
     346        If flag is 2, the quantity is considered time dependent and
     347        it will be stored at each yieldstep by appending it to the
     348        appropriate 2D array in the sww file.   
     349       
    336350        If q is None, storage will be switched off altogether.
     351       
     352        Once the simulation has started and thw sww file opened,
     353        this function will have no effect.
     354       
     355        The format, where q is a list of names is for backwards compatibility only.
     356        It will take the specified quantities to be time dependent and assume
     357        'elevation' to be static regardless.
    337358        """
    338359
    339360        if q is None:
    340             self.quantities_to_be_stored = []
     361            self.quantities_to_be_stored = {}
    341362            self.store = False
    342363            return
    343 
    344         if isinstance(q, basestring):
    345             q = [q] # Turn argument into a list
    346364
    347365        # Check correcness
    … …  
    349367            msg = ('Quantity %s is not a valid conserved quantity'
    350368                   % quantity_name)
    351             assert quantity_name in self.conserved_quantities, msg
    352 
     369            assert quantity_name in self.quantities, msg
     370
     371        if type(q) == types.ListType:
     372
     373            msg = 'List arguments to set_quantities_to_be_stored '
     374            msg += 'has been deprecated and will be removed in future '
     375            msg += 'versions of ANUGA.'
     376            msg += 'Please use dictionary argument instead'
     377            warn(msg, DeprecationWarning)
     378
     379       
     380       
     381            # FIXME: Raise deprecation warning
     382            tmp = {}
     383            for x in q:
     384                tmp[x] = 2
     385            tmp['elevation'] = 1   
     386            q = tmp     
     387           
     388        assert type(q) == types.DictType   
    353389        self.quantities_to_be_stored = q
    354390
    … …  
    658694
    659695        from anuga.shallow_water.data_manager import SWW_file
    660 
     696       
    661697        # Initialise writer
    662698        self.writer = SWW_file(self)

  • anuga_core/source/anuga/shallow_water/test_data_manager.py

    r7340 r7342  
    14511451
    14521452        self.domain.set_datadir('.')
    1453         self.domain.format = 'sww'
    14541453        self.domain.smooth = True
    14551454        self.domain.set_quantity('elevation', lambda x,y: -x-y)
    … …  
    15321531
    15331532        self.domain.set_datadir('.')
    1534         self.domain.format = 'sww'
    15351533        self.domain.smooth = True
    15361534        self.domain.set_quantity('elevation', lambda x,y: -x-y)
    … …  
    36103608        domain.set_boundary({'left': Bd, 'right': Bd, 'top': Bd, 'bottom': Bd})
    36113609
    3612         domain.quantities_to_be_stored.extend(['xmomentum','ymomentum'])
     3610        domain.quantities_to_be_stored['xmomentum'] = 2
     3611        domain.quantities_to_be_stored['ymomentum'] = 2
    36133612        #Initial condition
    36143613        h = 0.05
    … …  
    36413640
    36423641        bits = ['vertex_coordinates']
    3643         for quantity in ['elevation']+domain.quantities_to_be_stored:
    3644             bits.append('get_quantity("%s").get_integral()' %quantity)
    3645             bits.append('get_quantity("%s").get_values()' %quantity)
     3642        for quantity in domain.quantities_to_be_stored:
     3643            bits.append('get_quantity("%s").get_integral()' % quantity)
     3644            bits.append('get_quantity("%s").get_values()' % quantity)
    36463645
    36473646        for bit in bits:
    … …  
    37303729        domain.set_datadir('.')
    37313730        domain.default_order=2
    3732         #domain.quantities_to_be_stored=['stage']
    37333731
    37343732        domain.set_quantity('elevation', lambda x,y: -x/3)
    … …  
    37803778                'vertex_coordinates']
    37813779
    3782         for quantity in ['elevation']+domain.quantities_to_be_stored:
     3780        for quantity in domain.quantities_to_be_stored:
    37833781            bits.append('get_quantity("%s").get_integral()' %quantity)
    37843782            bits.append('get_quantity("%s").get_values()' %quantity)
    … …  
    38513849        domain.set_boundary({'left': Bd, 'right': Bd, 'top': Bd, 'bottom': Bd})
    38523850
    3853         domain.quantities_to_be_stored.extend(['xmomentum','ymomentum'])
     3851        domain.quantities_to_be_stored['xmomentum'] = 2
     3852        domain.quantities_to_be_stored['ymomentum'] = 2       
    38543853        #Initial condition
    38553854        h = 0.05
    … …  
    38843883
    38853884
    3886         bits = []#'vertex_coordinates']
    3887         for quantity in ['elevation']+domain.quantities_to_be_stored:
     3885        bits = []
     3886        for quantity in domain.quantities_to_be_stored:
    38883887            bits.append('quantities["%s"].get_integral()'%quantity)
    38893888
    … …  
    1045610455                         verbose=self.verbose,sww_precision=netcdf_float)
    1045710456        sww.store_triangulation(outfile, points_utm, volumes,
    10458                                 elevation,  new_origin=new_origin,
     10457                                new_origin=new_origin,
    1045910458                                verbose=self.verbose)
     10459        sww.store_static_quantities(outfile, elevation=elevation)                               
     10460                               
    1046010461        outfile.close()
    1046110462        fid = NetCDFFile(filename)

  • anuga_core/source/anuga/shallow_water/test_shallow_water_domain.py

    r7340 r7342  
    37323732        domain.smooth = False
    37333733        domain.default_order = 2
    3734         domain.set_quantities_to_be_stored(['stage', 'xmomentum', 'ymomentum'])
     3734
    37353735
    37363736        # IC
    … …  
    50955095        # FIXME: This is extremely important!
    50965096        # How can we test if they weren't stored?
    5097         domain1.quantities_to_be_stored = ['stage', 'xmomentum', 'ymomentum']
    50985097
    50995098
    … …  
    52005199        # FIXME: This is extremely important!
    52015200        # How can we test if they weren't stored?
    5202         domain1.quantities_to_be_stored = ['stage', 'xmomentum', 'ymomentum']
     5201        domain1.quantities_to_be_stored = {'elevation': 1,
     5202                                           'stage': 2,
     5203                                           'xmomentum': 2,
     5204                                           'ymomentum': 2}
    52035205
    52045206        # Bed-slope and friction at vertices (and interpolated elsewhere)
    … …  
    53735375        # FIXME: This is extremely important!
    53745376        # How can we test if they weren't stored?
    5375         domain1.quantities_to_be_stored = ['stage', 'xmomentum', 'ymomentum']
    53765377
    53775378        # Bed-slope and friction at vertices (and interpolated elsewhere)
    … …  
    55455546        # FIXME: This is extremely important!
    55465547        # How can we test if they weren't stored?
    5547         domain1.quantities_to_be_stored = ['stage', 'xmomentum', 'ymomentum']
     5548        domain1.set_quantities_to_be_stored({'elevation': 1,
     5549                                             'stage': 2,
     5550                                             'xmomentum': 2,
     5551                                             'ymomentum': 2})
    55485552
    55495553
Note: See TracChangeset for help on using the changeset viewer.