Changeset 8284

Timestamp:

Dec 15, 2011, 5:58:24 PM (13 years ago)

Author:

steve

Message:

latest ssw_merge (not working though)

Location:

trunk/anuga_core/source

Files:

• anuga/utilities/sww_merge.py (modified) (3 diffs)
• anuga_parallel/run_parallel_sw_merimbula.py (modified) (4 diffs)

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

@@ -17 +17 @@
 
     _sww_merge(swwfiles, output, verbose)
+
+
+def sww_merge_parallel(domain_global_name, np, verbose=False):
+
+    output = domain_global_name+".sww"
+    swwfiles = [ domain_global_name+"_P"+str(np)+"_"+str(v)+".sww" for v in range(np)]
+
+    _sww_merge_parallel(swwfiles, output, verbose)
 
 
@@ -122 +130 @@
 
     if verbose:
-            print 'Writing file ', output, ':'
+        print 'Writing file ', output, ':'
     fido = NetCDFFile(output, netcdf_mode_w)
     sww = Write_sww(static_quantities, dynamic_quantities)
@@ -166 +174 @@
                                         
     fido.close()
-    
+
+
+def _sww_merge_parallel(swwfiles, output, verbose):
+    """
+        Merge a list of sww files into a single file.
+        
+        Use to merge files created by parallel runs.
+
+        The sww files to be merged must have exactly the same timesteps.
+
+        Note that some advanced information and custom quantities may not be
+        exported.
+        
+        swwfiles is a list of .sww files to merge.
+        output is the output filename, including .sww extension.
+        verbose True to log output information
+    """
+
+    if verbose:
+        print "MERGING SWW Files"
+        
+    static_quantities = ['elevation']
+    dynamic_quantities = ['stage', 'xmomentum', 'ymomentum']
+    
+    first_file = True
+    tri_offset = 0
+    for filename in swwfiles:
+        if verbose:
+            print 'Reading file ', filename, ':'    
+    
+        fid = NetCDFFile(filename, netcdf_mode_r)
+         
+        if first_file:
+
+            times    = fid.variables['time'][:]
+            number_of_timesteps = len(times)
+            
+            out_s_quantities = {}
+            out_d_quantities = {}
+
+
+            xllcorner = fid.xllcorner
+            yllcorner = fid.yllcorner
+
+            number_of_global_triangles = int(fid.number_of_global_triangles)
+            number_of_global_nodes     = int(fid.number_of_global_nodes)
+
+            order      = fid.order
+            xllcorner  = fid.xllcorner;
+            yllcorner  = fid.yllcorner ;
+            zone       = fid.zone;
+            false_easting  = fid.false_easting;
+            false_northing = fid.false_northing;
+            datum      = fid.datum;
+            projection = fid.projection;
+
+            g_volumes = num.zeros((number_of_global_triangles,3),num.int)
+            g_x = num.zeros((number_of_global_nodes,),num.float32)
+            g_y = num.zeros((number_of_global_nodes,),num.float32)
+
+            g_points = num.zeros((number_of_global_nodes,2),num.float32)
+
+            for quantity in static_quantities:
+                out_s_quantities[quantity] = num.zeros((number_of_global_nodes,),num.float32)
+
+            # Quantities are stored as a 2D array of timesteps x data.
+            for quantity in dynamic_quantities:
+                out_d_quantities[quantity] = \
+                      num.zeros((number_of_timesteps,number_of_global_nodes),num.float32)
+                 
+            description = 'merged:' + getattr(fid, 'description')          
+            first_file = False
+
+
+        # Read in from files and add to global arrays
+
+        tri_l2g  = fid.variables['tri_l2g'][:]
+        node_l2g = fid.variables['node_l2g'][:]
+        tri_full_flag = fid.variables['tri_full_flag'][:]
+        volumes = fid.variables['volumes'][:]
+
+        # Just pick out the full triangles
+        ftri_l2g = num.compress(tri_full_flag, tri_l2g)
+        g_volumes[ftri_l2g] = num.compress(tri_full_flag,volumes,axis=0)
+
+        #g_x[node_l2g] = fid.variables['x']
+        #g_y[node_l2g] = fid.variables['y']
+
+        g_points[node_l2g,0] = fid.variables['x']
+        g_points[node_l2g,1] = fid.variables['y']
+        
+
+        print number_of_timesteps
+        
+        # Read in static quantities
+        for quantity in static_quantities:
+            out_s_quantities[quantity][node_l2g] = \
+                         num.array(fid.variables[quantity],dtype=num.float32)
+
+        
+        #Collate all dynamic quantities according to their timestep
+        for quantity in dynamic_quantities:
+            q = fid.variables[quantity]
+            print q.shape
+            for i in range(number_of_timesteps):
+                out_d_quantities[quantity][i][node_l2g] = \
+                           num.array(q[i],dtype=num.float32)
+
+
+
+
+        fid.close()
+
+
+    #---------------------------
+    # Write out the SWW file
+    #---------------------------
+    print g_points.shape
+
+    print number_of_global_triangles
+    print number_of_global_nodes
+
+
+    if verbose:
+            print 'Writing file ', output, ':'
+    fido = NetCDFFile(output, netcdf_mode_w)
+    sww = Write_sww(static_quantities, dynamic_quantities)
+    sww.store_header(fido, times,
+                             number_of_global_triangles,
+                             number_of_global_nodes,
+                             description=description,
+                             sww_precision=netcdf_float32)
+
+
+
+    from anuga.coordinate_transforms.geo_reference import Geo_reference
+    geo_reference = Geo_reference()
+    
+    sww.store_triangulation(fido, g_points, g_volumes, points_georeference=geo_reference)
+
+    fido.order      = order
+    fido.xllcorner  = xllcorner;
+    fido.yllcorner  = yllcorner ;
+    fido.zone       = zone;
+    fido.false_easting  = false_easting;
+    fido.false_northing = false_northing;
+    fido.datum      = datum;
+    fido.projection = projection;
+       
+    sww.store_static_quantities(fido, verbose=verbose, **out_s_quantities)
+
+    # Write out all the dynamic quantities for each timestep
+    for q in dynamic_quantities:
+        q_values = out_d_quantities[q]
+        for i in range(number_of_timesteps):
+            fido.variables[q][i] = q_values[i]
+        
+        # This updates the _range values
+        q_range = fido.variables[q + Write_sww.RANGE][:]
+        q_values_min = num.min(q_values)
+        if q_values_min < q_range[0]:
+            fido.variables[q + Write_sww.RANGE][0] = q_values_min
+        q_values_max = num.max(q_values)
+        if q_values_max > q_range[1]:
+            fido.variables[q + Write_sww.RANGE][1] = q_values_max        
+
+                                        
+    print out_s_quantities
+    print out_d_quantities
+    
+    print g_x
+    print g_y
+
+    #print g_volumes
+    
+    fido.close()
+
 
 if __name__ == "__main__":

trunk/anuga_core/source/anuga_parallel/run_parallel_sw_merimbula.py

@@ -45 +45 @@
 #mesh_filename = "merimbula_43200.tsh"   ; x0 = 756000.0 ; x1 = 756500.0
 #mesh_filename = "test-100.tsh" ; x0 = 0.25 ; x1 = 0.5
-mesh_filename = "test-20.tsh" ; x0 = 0.25 ; x1 = 0.5
+mesh_filename = "test-20.tsh" ; x0 = 250.0 ; x1 = 350.0
 yieldstep = 20
 finaltime = 40
@@ -65 +65 @@
         return self.h*((x>self.x0)&(x<self.x1))
 
+
+class Set_Elevation:
+    """Set an elevation
+    """
+
+    def __init__(self, h=1.0):
+        self.x0 = x0
+        self.x1 = x1
+        self.h  = h
+
+    def __call__(self, x, y):
+        return x/self.h
+    
+
 #--------------------------------------------------------------------------
 # Setup Domain only on processor 0
@@ -71 +85 @@
     domain = create_domain_from_file(mesh_filename)
     domain.set_quantity('stage', Set_Stage(x0, x1, 2.0))
+    domain.set_quantity('elevation', Set_Elevation(500.0))
 
     print domain.statistics()
@@ -100 +115 @@
         print domain.get_extent(absolute=True)
         print domain.geo_reference
-        #print domain.tri_map
-        #print domain.inv_tri_map
+        print domain.s2p_map
+        print domain.p2s_map
         print domain.tri_l2g
         print domain.node_l2g