Ignore:
Timestamp:
Aug 9, 2009, 1:20:42 AM (14 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.

File:
1 edited

Legend:

Unmodified
Added
Removed
  • 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)
Note: See TracChangeset for help on using the changeset viewer.