Changeset 7819


Ignore:
Timestamp:
Jun 11, 2010, 10:59:25 AM (10 years ago)
Author:
hudson
Message:

Semi-working code for SWW merge util.

File:
1 edited

Legend:

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

    r7817 r7819  
    11import numpy as num
     2from anuga.utilities.numerical_tools import ensure_numeric
    23
    34from Scientific.IO.NetCDF import NetCDFFile
    4 from anuga.config import netcdf_mode_r, netcdf_mode_w, netcdf_mode_a
     5from anuga.config import netcdf_mode_r, netcdf_mode_w, \
     6                            netcdf_mode_a, netcdf_float
    57from anuga.file.sww import SWW_file, Write_sww
     8
    69
    710def sww_merge(swwfiles, output, verbose = False):
     
    912        Merge a list of sww files into a single file.
    1013    """
     14   
     15    static_quantities = ['elevation']
     16    dynamic_quantities = ['stage', 'xmomentum', 'ymomentum']
    1117   
    1218    first_file = True
     
    2430            y = []
    2531            out_tris = list(tris) 
     32            out_s_quantities = {}
     33            out_d_quantities = {}
     34           
     35            for quantity in static_quantities:
     36                out_s_quantities[quantity] = []
     37
     38            for quantity in dynamic_quantities:
     39                out_d_quantities[quantity] = [ [] for _ in range(len(times))]
     40                 
    2641            description = 'merged:' + getattr(fid, 'description')         
    2742            first_file = False
     
    3146                out_tris.append(verts)
    3247
    33         tri_offset += fid.dimensions['number_of_points']
     48        num_pts = fid.dimensions['number_of_points']
     49        tri_offset += num_pts
    3450       
    3551        if verbose:
     
    3854        x.extend(list(fid.variables['x'][:]))
    3955        y.extend(list(fid.variables['y'][:]))
    40      
    41     points = [[xx, yy] for xx, yy in zip(x, y)]   
    42 
    43 
     56       
     57        for quantity in static_quantities:
     58            out_s_quantities[quantity].extend(fid.variables[quantity][:])
     59           
     60        for quantity in dynamic_quantities:
     61            time_chunks = fid.variables[quantity][:]
     62            for i, time_chunk in enumerate(time_chunks):
     63                out_d_quantities[quantity][i].extend(time_chunk)           
     64       
     65    points = [[xx, yy] for xx, yy in zip(x, y)]
     66    fid.close()
     67   
    4468    # NetCDF file definition
    4569    fido = NetCDFFile(output, netcdf_mode_w)
    46     sww = Write_sww(['elevation'], ['stage', 'xmomentum', 'ymomentum'])
     70    sww = Write_sww(static_quantities, dynamic_quantities)
    4771    sww.store_header(fido, times,
    4872                             len(out_tris),
    4973                             len(points),
    50                              description=description)
    51 
    52     sww.store_triangulation(fido,
    53                                     points,
    54                                     num.array(out_tris).astype(num.float32))
     74                             description=description,
     75                             sww_precision=netcdf_float)
    5576
    5677
    57     # Get names of static quantities
    58     static_quantities = {}
    59     for name in ['elevation']:
    60         static_quantities[name] = fid.variables[name][:]
    61    
    62     # Store static quantities       
    63     self.writer.store_static_quantities(fido, **static_quantities)
     78    sww.store_triangulation(fido, points, out_tris)
     79       
     80    sww.store_static_quantities(fido, verbose=verbose, **out_s_quantities)
     81
     82    for q in dynamic_quantities:
     83        q_values = ensure_numeric(out_d_quantities[quantity])
     84        x = q_values.astype(netcdf_float)
     85        fido.variables[q] = x
     86
     87        # This updates the _range values
     88        q_range = fido.variables[q + Write_sww.RANGE][:]
     89        q_values_min = num.min(q_values)
     90        if q_values_min < q_range[0]:
     91            fido.variables[q + Write_sww.RANGE][0] = q_values_min
     92        q_values_max = num.max(q_values)
     93        if q_values_max > q_range[1]:
     94            fido.variables[q + Write_sww.RANGE][1] = q_values_max       
     95
    6496                                       
    65     fid.close()
    66 
    67 
    68 
    69     domain = Domain(points, out_tris)
    70     domain.set_name(output)
    71     sww = SWW_file(domain)
    72     sww.store_connectivity()
    73    
    74     for _ in times:
    75         sww.store_timestep()
     97    fido.close()
    7698   
    7799
     
    80102from anuga.shallow_water.shallow_water_domain import Domain
    81103from anuga.file.sww import SWW_file
     104from anuga.abstract_2d_finite_volumes.generic_boundary_conditions import \
     105    Dirichlet_boundary
     106
     107Bd = Dirichlet_boundary([0.5, 0., 0.])
    82108
    83109# Create shallow water domain
    84110domain = Domain(*rectangular_cross(2, 2))
    85111domain.set_name('test1')
    86 sww = SWW_file(domain)
    87 sww.store_connectivity()
    88 
     112domain.set_quantity('elevation', 2)
     113domain.set_quantity('stage', 5)
     114domain.set_boundary({'left': Bd, 'right': Bd, 'top': Bd, 'bottom': Bd})
     115for t in domain.evolve(yieldstep=0.5, finaltime=1):
     116    pass
     117   
    89118domain = Domain(*rectangular(3, 3))
    90119domain.set_name('test2')
    91 sww = SWW_file(domain)
    92 sww.store_connectivity()
    93 
     120domain.set_quantity('elevation', 3)
     121domain.set_quantity('stage', 50)
     122domain.set_boundary({'left': Bd, 'right': Bd, 'top': Bd, 'bottom': Bd})
     123for t in domain.evolve(yieldstep=0.5, finaltime=1):
     124    pass
    94125       
    95126sww_merge(['test1.sww', 'test2.sww'], 'test_out.sww', verbose = True)
Note: See TracChangeset for help on using the changeset viewer.