Changeset 9019


Ignore:
Timestamp:
Nov 9, 2013, 6:13:58 PM (11 years ago)
Author:
steve
Message:

Added storing of centroids to the sww file

Location:
trunk/anuga_core/source
Files:
5 edited

Legend:

Unmodified
Added
Removed

  • trunk/anuga_core/source/anuga/file/sww.py

    r9018 r9019  
    9393        static_quantities = []
    9494        dynamic_quantities = []
     95        static_c_quantities = []
     96        dynamic_c_quantities = []
    9597       
    9698        for q in domain.quantities_to_be_stored:
    … …  
    102104       
    103105            assert flag in [1,2]
    104             if flag == 1: static_quantities.append(q)
    105             if flag == 2: dynamic_quantities.append(q)               
     106            if flag == 1:
     107                static_quantities.append(q)
     108                if self.store_centroids: static_c_quantities.append(q+'_c')
     109               
     110            if flag == 2:
     111                dynamic_quantities.append(q)
     112                if self.store_centroids: dynamic_c_quantities.append(q+'_c')
    106113                       
    107114       
    … …  
    112119                          'suitable for plotting'
    113120                         
    114             self.writer = Write_sww(static_quantities, dynamic_quantities, store_centroids=self.store_centroids)
     121            self.writer = Write_sww(static_quantities,
     122                                    dynamic_quantities,
     123                                    static_c_quantities,
     124                                    dynamic_c_quantities)
     125           
    115126            self.writer.store_header(fid,
    116127                                     domain.starttime,
    … …  
    187198                                        V.astype(num.float32),
    188199                                        points_georeference=\
    189                                             domain.geo_reference)
     200                                        domain.geo_reference)
    190201
    191202
    … …  
    202213        static_quantities = {}
    203214        static_quantities_centroid = {}
     215       
    204216        for name in self.writer.static_quantities:
    205217            Q = domain.quantities[name]
    … …  
    207219                                       precision=self.precision)
    208220            static_quantities[name] = A
    209            
    210             if self.store_centroids:
    211                 static_quantities_centroid[name] = Q.centroid_values
     221
     222        #print domain.quantities
     223        #print self.writer.static_c_quantities
     224
     225        for name in self.writer.static_c_quantities:
     226            Q = domain.quantities[name[:-2]]  # rip off _c from name
     227            static_quantities_centroid[name] = Q.centroid_values
    212228       
    213229        # Store static quantities       
    214230        self.writer.store_static_quantities(fid, **static_quantities)
    215        
    216         if  self.store_centroids:
    217             self.writer.store_static_quantities_centroid(fid, **static_quantities_centroid)
     231        self.writer.store_static_quantities_centroid(fid, **static_quantities_centroid)
    218232                                       
    219233        fid.close()
    … …  
    357371                dynamic_quantities[name] = A
    358372               
    359                 if self.store_centroids:
    360                     dynamic_quantities_centroid[name] = Q.centroid_values
     373            for name in self.writer.dynamic_c_quantities:
     374                Q = domain.quantities[name[:-2]]
     375                dynamic_quantities_centroid[name] = Q.centroid_values
    361376               
    362377                                       
    … …  
    493508    """
    494509
    495     def __init__(self, static_quantities, dynamic_quantities, store_centroids=False):
    496         """Initialise Write_sww with two list af quantity names:
     510    def __init__(self,
     511                 static_quantities,
     512                 dynamic_quantities,
     513                 static_c_quantities = [],
     514                 dynamic_c_quantities = []):
     515       
     516        """Initialise Write_sww with two (or 4) list af quantity names:
    497517       
    498518        static_quantities (e.g. elevation or friction):
    … …  
    501521        dynamic_quantities (e.g stage):
    502522            Stored every timestep in a 2D array with
    503             dimensions number_of_points X number_of_timesteps       
     523            dimensions number_of_points X number_of_timesteps
     524
     525        static_c_quantities (e.g. elevation_c or friction_c):
     526            Stored once at the beginning of the simulation in a 1D array
     527            of length number_of_triangles 
     528        dynamic_c_quantities (e.g stage_c):
     529            Stored every timestep in a 2D array with
     530            dimensions number_of_triangles X number_of_timesteps
    504531       
    505532        """
    506533        self.static_quantities = static_quantities   
    507534        self.dynamic_quantities = dynamic_quantities
    508         self.store_centroids = store_centroids
     535        self.static_c_quantities = static_c_quantities
     536        self.dynamic_c_quantities = dynamic_c_quantities
     537
     538        self.store_centroids = False
     539        if static_c_quantities or dynamic_c_quantities:
     540            self.store_centroids = True
    509541
    510542
    … …  
    606638                                   ('numbers_in_range',))
    607639           
    608             if self.store_centroids:
    609                 outfile.createVariable(q + '_c', sww_precision,
    610                                    ('number_of_volumes',))
    611                                    
    612640            # Initialise ranges with small and large sentinels.
    613641            # If this was in pure Python we could have used None sensibly
    … …  
    615643            outfile.variables[q+Write_sww.RANGE][1] = -max_float # Max
    616644
    617        
    618         self.write_dynamic_quantities(outfile, self.dynamic_quantities, times, \
    619                                         precis = sww_precision)
     645
     646        for q in self.static_c_quantities:
     647            outfile.createVariable(q, sww_precision,
     648                                   ('number_of_volumes',))
     649                                   
     650
     651        self.write_dynamic_quantities(outfile, times, precis = sww_precision)
    620652
    621653
    … …  
    734766
    735767
    736     def write_dynamic_quantities(self, outfile, quantities,
    737                     times, precis = netcdf_float32, verbose = False):   
     768    def write_dynamic_quantities(self, outfile,
     769                                 times, precis = netcdf_float32, verbose = False):   
    738770        """
    739771            Write out given quantities to file.
    … …  
    741773       
    742774
    743         for q in quantities:
     775        for q in self.dynamic_quantities:
    744776            outfile.createVariable(q, precis, ('number_of_timesteps',
    745                                                       'number_of_points'))
     777                                               'number_of_points'))
    746778            outfile.createVariable(q + Write_sts.RANGE, precis,
    747779                                   ('numbers_in_range',))
    748780           
    749             if self.store_centroids:
    750                 outfile.createVariable(q + '_c', precis,
    751                                    ('number_of_timesteps','number_of_volumes'))
    752 
    753781            # Initialise ranges with small and large sentinels.
    754782            # If this was in pure Python we could have used None sensibly
    755783            outfile.variables[q+Write_sts.RANGE][0] = max_float  # Min
    756784            outfile.variables[q+Write_sts.RANGE][1] = -max_float # Max
     785
     786        for q in self.dynamic_c_quantities:
     787            outfile.createVariable(q, precis, ('number_of_timesteps',
     788                                                    'number_of_volumes'))
    757789
    758790        # Doing sts_precision instead of Float gives cast errors.
    … …  
    888920        # This method will also write the ranges for each quantity,
    889921        # e.g. stage_range, xmomentum_range and ymomentum_range
    890         for q in self.static_quantities:
     922
     923        #print outfile.variables.keys()
     924        #print self.static_c_quantities
     925       
     926        for q in self.static_c_quantities:
    891927            if not quant.has_key(q):
    892928                msg = 'Values for quantity %s was not specified in ' % q
    … …  
    897933               
    898934                x = q_values.astype(sww_precision)
    899                 outfile.variables[q+'_c'][:] = x
     935                outfile.variables[q][:] = x
    900936       
    901937
    … …  
    10141050        # This method will also write the ranges for each quantity,
    10151051        # e.g. stage_range, xmomentum_range and ymomentum_range
    1016         for q in self.dynamic_quantities:
     1052       
     1053        #print 50*"="
     1054        #print quant
     1055        #print self.dynamic_c_quantities
     1056
     1057        for q in self.dynamic_c_quantities:
    10171058            if not quant.has_key(q):
    10181059                msg = 'Values for quantity %s was not specified in ' % q
    … …  
    10231064               
    10241065                q_retyped = q_values.astype(sww_precision)
    1025                 outfile.variables[q+'_c'][slice_index] = q_retyped
     1066                outfile.variables[q][slice_index] = q_retyped
    10261067
    10271068

  • trunk/anuga_core/source/anuga/file/test_read_sww.py

    r9018 r9019  
    177177
    178178        # Clean up
    179         os.remove(source)
     179        #os.remove(source)
    180180       
    181181    def test_read_sww_with_centroids(self):
    … …  
    277277        assert num.allclose(sww_file.get_bounds(), [0.0, length, 0.0, width])
    278278
     279        #print 50*"="
     280        #print sww_file.quantities.keys()
     281
    279282        assert 'stage'     in sww_file.quantities.keys()
    280283        assert 'friction'  in sww_file.quantities.keys()

  • trunk/anuga_core/source/anuga/utilities/plot_utils.py

    r8950 r9019  
    1 """ Random utilities for reading sww file data and for plotting (in ipython, or
    2     in scripts)
     1""" Random utilities for reading sww file data and for plotting
     2(in ipython, or in scripts)
    33
    44    Functionality of note:
    55
    6     util.get_outputs -- read the data from a single sww file into a single object
    7    
    8     util.combine_outputs -- read the data from a list of sww files into a single object
    9    
    10     util.near_transect -- for finding the indices of points 'near' to a given
    11                        line, and assigning these points a coordinate along that line.
    12                        This is useful for plotting outputs which are 'almost' along a
    13                        transect (e.g. a channel cross-section) -- see example below
     6    util.get_outputs -- read the data from a single sww file
     7    into a single object
     8   
     9    util.combine_outputs -- read the data from a list of sww
     10    files into a single object
     11   
     12    util.near_transect -- for finding the indices of points
     13                          'near' to a given line, and
     14                          assigning these points a
     15                          coordinate along that line.
     16
     17    This is useful for plotting outputs which are 'almost' along a
     18    transect (e.g. a channel cross-section) -- see example below
    1419
    1520    util.sort_sww_filenames -- match sww filenames by a wildcard, and order
    1621                               them according to their 'time'. This means that
    17                                they can be stuck together using 'combine_outputs' correctly
    18 
    19     util.triangle_areas -- compute the areas of every triangle in a get_outputs object [ must be vertex-based]
    20 
    21     util.water_volume -- compute the water volume at every time step in an sww file (needs both vertex and centroid value input).
     22                               they can be stuck together using
     23                               'combine_outputs' correctly
     24
     25    util.triangle_areas -- compute the areas of every triangle
     26                           in a get_outputs object [ must be vertex-based]
     27
     28    util.water_volume -- compute the water volume at every
     29                         time step in an sww file (needs both
     30                         vertex and centroid value input).
    2231 
    2332    Here is an example ipython session which uses some of these functions:
    … …  
    417426###
    418427
    419 def water_volume(p,p2, per_unit_area=False, subset=None):
     428def water_volume(p, p2, per_unit_area=False, subset=None):
    420429    # Compute the water volume from p(vertex values) and p2(centroid values)
    421430

  • trunk/anuga_core/source/anuga/utilities/sww_merge.py

    r8809 r9019  
    240240            out_d_quantities = {}
    241241
     242            out_s_c_quantities = {}
     243            out_d_c_quantities = {}
     244
    242245
    243246            xllcorner = fid.xllcorner
    … …  
    262265            g_points = num.zeros((number_of_global_nodes,2),num.float32)
    263266
     267            #=====================================
     268            # Deal with the vertex based variables
     269            #=====================================
    264270            quantities = set(['elevation', 'friction', 'stage', 'xmomentum',
    265271                              'ymomentum', 'xvelocity', 'yvelocity', 'height'])
    … …  
    272278
    273279            for quantity in quantities:
    274                 # Test if elevation is static
     280                # Test if quantity is static
    275281                if n_steps == fid.variables[quantity].shape[0]:
    276282                    dynamic_quantities.append(quantity)
    … …  
    285291                out_d_quantities[quantity] = \
    286292                      num.zeros((n_steps,number_of_global_nodes),num.float32)
     293
     294            #=======================================
     295            # Deal with the centroid based variables
     296            #=======================================
     297            quantities = set(['elevation_c', 'friction_c', 'stage_c', 'xmomentum_c',
     298                              'ymomentum_c', 'xvelocity_c', 'yvelocity_c', 'height_c'])
     299            variables = set(fid.variables.keys())
     300
     301            quantities = list(quantities & variables)
     302           
     303            static_c_quantities = []
     304            dynamic_c_quantities = []
     305
     306            for quantity in quantities:
     307                # Test if quantity is static
     308                if n_steps == fid.variables[quantity].shape[0]:
     309                    dynamic_c_quantities.append(quantity)
     310                else:
     311                    static_c_quantities.append(quantity)
     312               
     313            for quantity in static_c_quantities:
     314                out_s_c_quantities[quantity] = num.zeros((number_of_global_triangles,),num.float32)
     315
     316            # Quantities are stored as a 2D array of timesteps x data.
     317            for quantity in dynamic_c_quantities:
     318                out_d_c_quantities[quantity] = \
     319                      num.zeros((n_steps,number_of_global_triangles),num.float32)
    287320                 
    288321            description = 'merged:' + getattr(fid, 'description')         
    … …  
    321354
    322355        # Just pick out the full triangles
     356        ftri_ids = num.where(tri_full_flag>0)
    323357        ftri_l2g = num.compress(tri_full_flag, tri_l2g)
    324358
    … …  
    379413
    380414
     415        # Read in static c quantities
     416        for quantity in static_c_quantities:
     417            #out_s_quantities[quantity][node_l2g] = \
     418            #             num.array(fid.variables[quantity],dtype=num.float32)
     419            q = fid.variables[quantity]
     420            out_s_c_quantities[quantity][ftri_l2g] = \
     421                         num.array(q,dtype=num.float32)[ftri_ids]
     422
     423       
     424        #Collate all dynamic c quantities according to their timestep
     425        for quantity in dynamic_c_quantities:
     426            q = fid.variables[quantity]
     427            #print q.shape
     428            for i in range(n_steps):
     429                out_d_c_quantities[quantity][i][ftri_l2g] = \
     430                           num.array(q[i],dtype=num.float32)[ftri_ids]
    381431
    382432
    … …  
    396446            print 'Writing file ', output, ':'
    397447    fido = NetCDFFile(output, netcdf_mode_w)
    398     sww = Write_sww(static_quantities, dynamic_quantities)
     448
     449    sww = Write_sww(static_quantities, dynamic_quantities, static_c_quantities, dynamic_c_quantities)
    399450    sww.store_header(fido, starttime,
    400451                             number_of_global_triangles,
    … …  
    420471       
    421472    sww.store_static_quantities(fido, verbose=verbose, **out_s_quantities)
     473    sww.store_static_quantities_centroid(fido, verbose=verbose, **out_s_c_quantities)
    422474
    423475    # Write out all the dynamic quantities for each timestep
    … …  
    439491        if q_values_max > q_range[1]:
    440492            fido.variables[q + Write_sww.RANGE][1] = q_values_max       
     493
     494    for q in dynamic_c_quantities:
     495        q_values = out_d_c_quantities[q]
     496        for i in range(n_steps):
     497            fido.variables[q][i] = q_values[i]
    441498
    442499                                       

  • trunk/anuga_core/source/anuga_parallel/run_parallel_sw_merimbula.py

    r8950 r9019  
    4343#--------------------------------------------------------------------------
    4444
    45 mesh_filename = "merimbula_10785_1.tsh" ; x0 = 756000.0 ; x1 = 756500.0; yieldstep = 10; finaltime = 10
     45#mesh_filename = "merimbula_10785_1.tsh" ; x0 = 756000.0 ; x1 = 756500.0; yieldstep = 10; finaltime = 10
    4646#mesh_filename = "merimbula_17156.tsh"   ; x0 = 756000.0 ; x1 = 756500.0; yieldstep = 50; finaltime = 500
    4747#mesh_filename = "merimbula_43200_1.tsh"   ; x0 = 756000.0 ; x1 = 756500.0; yieldstep = 50; finaltime = 500
    48 #mesh_filename = "test-100.tsh" ; x0 = 200.0 ; x1 = 300.0; yieldstep = 1; finaltime = 50
     48mesh_filename = "test-100.tsh" ; x0 = 200.0 ; x1 = 300.0; yieldstep = 1; finaltime = 10
    4949#mesh_filename = "test-20.tsh" ; x0 = 250.0 ; x1 = 350.0; yieldstep = 1; finaltime = 50
    5050
    … …  
    8585if myid == 0:
    8686    domain = create_domain_from_file(mesh_filename)
    87     #domain.set_quantity('stage', Set_Stage(x0, x1, 1.0))
    88     domain.set_quantity('stage', 0.0)
     87    domain.set_quantity('stage', Set_Stage(x0, x1, 1.0))
     88    #domain.set_quantity('stage', 0.0)
    8989    #domain.set_datadir('.')
    9090    domain.set_name('merimbula_new')
    … …  
    113113
    114114domain.set_flow_algorithm('2_0')
    115 
     115domain.set_store_centroids()
    116116
    117117#for p in range(numprocs):
Note: See TracChangeset for help on using the changeset viewer.