Changeset 8284


Ignore:
Timestamp:
Dec 15, 2011, 5:58:24 PM (8 years ago)
Author:
steve
Message:

latest ssw_merge (not working though)

Location:
trunk/anuga_core/source
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • trunk/anuga_core/source/anuga/utilities/sww_merge.py

    r8281 r8284  
    1717
    1818    _sww_merge(swwfiles, output, verbose)
     19
     20
     21def sww_merge_parallel(domain_global_name, np, verbose=False):
     22
     23    output = domain_global_name+".sww"
     24    swwfiles = [ domain_global_name+"_P"+str(np)+"_"+str(v)+".sww" for v in range(np)]
     25
     26    _sww_merge_parallel(swwfiles, output, verbose)
    1927
    2028
     
    122130
    123131    if verbose:
    124             print 'Writing file ', output, ':'
     132        print 'Writing file ', output, ':'
    125133    fido = NetCDFFile(output, netcdf_mode_w)
    126134    sww = Write_sww(static_quantities, dynamic_quantities)
     
    166174                                       
    167175    fido.close()
    168    
     176
     177
     178def _sww_merge_parallel(swwfiles, output, verbose):
     179    """
     180        Merge a list of sww files into a single file.
     181       
     182        Use to merge files created by parallel runs.
     183
     184        The sww files to be merged must have exactly the same timesteps.
     185
     186        Note that some advanced information and custom quantities may not be
     187        exported.
     188       
     189        swwfiles is a list of .sww files to merge.
     190        output is the output filename, including .sww extension.
     191        verbose True to log output information
     192    """
     193
     194    if verbose:
     195        print "MERGING SWW Files"
     196       
     197    static_quantities = ['elevation']
     198    dynamic_quantities = ['stage', 'xmomentum', 'ymomentum']
     199   
     200    first_file = True
     201    tri_offset = 0
     202    for filename in swwfiles:
     203        if verbose:
     204            print 'Reading file ', filename, ':'   
     205   
     206        fid = NetCDFFile(filename, netcdf_mode_r)
     207         
     208        if first_file:
     209
     210            times    = fid.variables['time'][:]
     211            number_of_timesteps = len(times)
     212           
     213            out_s_quantities = {}
     214            out_d_quantities = {}
     215
     216
     217            xllcorner = fid.xllcorner
     218            yllcorner = fid.yllcorner
     219
     220            number_of_global_triangles = int(fid.number_of_global_triangles)
     221            number_of_global_nodes     = int(fid.number_of_global_nodes)
     222
     223            order      = fid.order
     224            xllcorner  = fid.xllcorner;
     225            yllcorner  = fid.yllcorner ;
     226            zone       = fid.zone;
     227            false_easting  = fid.false_easting;
     228            false_northing = fid.false_northing;
     229            datum      = fid.datum;
     230            projection = fid.projection;
     231
     232            g_volumes = num.zeros((number_of_global_triangles,3),num.int)
     233            g_x = num.zeros((number_of_global_nodes,),num.float32)
     234            g_y = num.zeros((number_of_global_nodes,),num.float32)
     235
     236            g_points = num.zeros((number_of_global_nodes,2),num.float32)
     237
     238            for quantity in static_quantities:
     239                out_s_quantities[quantity] = num.zeros((number_of_global_nodes,),num.float32)
     240
     241            # Quantities are stored as a 2D array of timesteps x data.
     242            for quantity in dynamic_quantities:
     243                out_d_quantities[quantity] = \
     244                      num.zeros((number_of_timesteps,number_of_global_nodes),num.float32)
     245                 
     246            description = 'merged:' + getattr(fid, 'description')         
     247            first_file = False
     248
     249
     250        # Read in from files and add to global arrays
     251
     252        tri_l2g  = fid.variables['tri_l2g'][:]
     253        node_l2g = fid.variables['node_l2g'][:]
     254        tri_full_flag = fid.variables['tri_full_flag'][:]
     255        volumes = fid.variables['volumes'][:]
     256
     257        # Just pick out the full triangles
     258        ftri_l2g = num.compress(tri_full_flag, tri_l2g)
     259        g_volumes[ftri_l2g] = num.compress(tri_full_flag,volumes,axis=0)
     260
     261        #g_x[node_l2g] = fid.variables['x']
     262        #g_y[node_l2g] = fid.variables['y']
     263
     264        g_points[node_l2g,0] = fid.variables['x']
     265        g_points[node_l2g,1] = fid.variables['y']
     266       
     267
     268        print number_of_timesteps
     269       
     270        # Read in static quantities
     271        for quantity in static_quantities:
     272            out_s_quantities[quantity][node_l2g] = \
     273                         num.array(fid.variables[quantity],dtype=num.float32)
     274
     275       
     276        #Collate all dynamic quantities according to their timestep
     277        for quantity in dynamic_quantities:
     278            q = fid.variables[quantity]
     279            print q.shape
     280            for i in range(number_of_timesteps):
     281                out_d_quantities[quantity][i][node_l2g] = \
     282                           num.array(q[i],dtype=num.float32)
     283
     284
     285
     286
     287        fid.close()
     288
     289
     290    #---------------------------
     291    # Write out the SWW file
     292    #---------------------------
     293    print g_points.shape
     294
     295    print number_of_global_triangles
     296    print number_of_global_nodes
     297
     298
     299    if verbose:
     300            print 'Writing file ', output, ':'
     301    fido = NetCDFFile(output, netcdf_mode_w)
     302    sww = Write_sww(static_quantities, dynamic_quantities)
     303    sww.store_header(fido, times,
     304                             number_of_global_triangles,
     305                             number_of_global_nodes,
     306                             description=description,
     307                             sww_precision=netcdf_float32)
     308
     309
     310
     311    from anuga.coordinate_transforms.geo_reference import Geo_reference
     312    geo_reference = Geo_reference()
     313   
     314    sww.store_triangulation(fido, g_points, g_volumes, points_georeference=geo_reference)
     315
     316    fido.order      = order
     317    fido.xllcorner  = xllcorner;
     318    fido.yllcorner  = yllcorner ;
     319    fido.zone       = zone;
     320    fido.false_easting  = false_easting;
     321    fido.false_northing = false_northing;
     322    fido.datum      = datum;
     323    fido.projection = projection;
     324       
     325    sww.store_static_quantities(fido, verbose=verbose, **out_s_quantities)
     326
     327    # Write out all the dynamic quantities for each timestep
     328    for q in dynamic_quantities:
     329        q_values = out_d_quantities[q]
     330        for i in range(number_of_timesteps):
     331            fido.variables[q][i] = q_values[i]
     332       
     333        # This updates the _range values
     334        q_range = fido.variables[q + Write_sww.RANGE][:]
     335        q_values_min = num.min(q_values)
     336        if q_values_min < q_range[0]:
     337            fido.variables[q + Write_sww.RANGE][0] = q_values_min
     338        q_values_max = num.max(q_values)
     339        if q_values_max > q_range[1]:
     340            fido.variables[q + Write_sww.RANGE][1] = q_values_max       
     341
     342                                       
     343    print out_s_quantities
     344    print out_d_quantities
     345   
     346    print g_x
     347    print g_y
     348
     349    #print g_volumes
     350   
     351    fido.close()
     352
    169353
    170354if __name__ == "__main__":
  • trunk/anuga_core/source/anuga_parallel/run_parallel_sw_merimbula.py

    r8283 r8284  
    4545#mesh_filename = "merimbula_43200.tsh"   ; x0 = 756000.0 ; x1 = 756500.0
    4646#mesh_filename = "test-100.tsh" ; x0 = 0.25 ; x1 = 0.5
    47 mesh_filename = "test-20.tsh" ; x0 = 0.25 ; x1 = 0.5
     47mesh_filename = "test-20.tsh" ; x0 = 250.0 ; x1 = 350.0
    4848yieldstep = 20
    4949finaltime = 40
     
    6565        return self.h*((x>self.x0)&(x<self.x1))
    6666
     67
     68class Set_Elevation:
     69    """Set an elevation
     70    """
     71
     72    def __init__(self, h=1.0):
     73        self.x0 = x0
     74        self.x1 = x1
     75        self.h  = h
     76
     77    def __call__(self, x, y):
     78        return x/self.h
     79   
     80
    6781#--------------------------------------------------------------------------
    6882# Setup Domain only on processor 0
     
    7185    domain = create_domain_from_file(mesh_filename)
    7286    domain.set_quantity('stage', Set_Stage(x0, x1, 2.0))
     87    domain.set_quantity('elevation', Set_Elevation(500.0))
    7388
    7489    print domain.statistics()
     
    100115        print domain.get_extent(absolute=True)
    101116        print domain.geo_reference
    102         #print domain.tri_map
    103         #print domain.inv_tri_map
     117        print domain.s2p_map
     118        print domain.p2s_map
    104119        print domain.tri_l2g
    105120        print domain.node_l2g
Note: See TracChangeset for help on using the changeset viewer.