Changeset 8281

Timestamp:

Dec 12, 2011, 5:11:47 PM (13 years ago)

Author:

steve

Message:

added in some parallel data to allow sww_merge to create a merged sww file

Location:

trunk/anuga_core/source

Files:

• anuga/abstract_2d_finite_volumes/generic_domain.py (modified) (1 diff)
• anuga/file/sww.py (modified) (2 diffs)
• anuga/utilities/sww_merge.py (modified) (2 diffs)
• anuga/utilities/test_file_utils.py (modified) (2 diffs)
• anuga_parallel/distribute_mesh.py (modified) (13 diffs)
• anuga_parallel/parallel_api.py (modified) (6 diffs)
• anuga_parallel/parallel_shallow_water.py (modified) (4 diffs)
• anuga_parallel/run_parallel_sw_merimbula.py (modified) (4 diffs)
• anuga_parallel/test_parallel_distribute_mesh.py (modified) (4 diffs)

trunk/anuga_core/source/anuga/abstract_2d_finite_volumes/generic_domain.py

@@ -165 +165 @@
         self.fractional_step_operators = []
 
+
+
+        # by default domain is not parallel
+        self.parallel = False
+
+        self.number_of_global_triangles = self.number_of_triangles
+        self.number_of_global_nodes = self.number_of_nodes
+
         # Setup the ghost cell communication
         if full_send_dict is None:

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

@@ -183 +183 @@
                                         points_georeference=\
                                             domain.geo_reference)
+
+
+        if domain.parallel:
+            self.writer.store_parallel_data(fid,
+                                        domain.number_of_global_triangles,
+                                        domain.number_of_global_nodes,
+                                        domain.tri_full_flag,
+                                        domain.tri_l2g,
+                                        domain.node_l2g)
 
 
@@ -676 +685 @@
         outfile.variables['y'][:] = y #- geo_ref.get_yllcorner()
         outfile.variables['volumes'][:] = volumes.astype(num.int32) #On Opteron 64
+
+
+    def store_parallel_data(self,
+                            outfile,
+                            number_of_global_triangles,
+                            number_of_global_nodes,
+                            tri_full_flag = None,
+                            tri_l2g = None,
+                            node_l2g = None,
+                            sww_precision=netcdf_float32,
+                            verbose=False):
+
+
+         # dimension definitions
+        #outfile.createDimension('number_of_volumes', number_of_volumes)
+        #outfile.createDimension('number_of_vertices', 3)
+        #outfile.createDimension('numbers_in_range', 2)
+
+        print 'store parallel data'
+        outfile.number_of_global_triangles = number_of_global_triangles
+        outfile.number_of_global_nodes = number_of_global_nodes
+
+        # variable definitions
+        outfile.createVariable('tri_l2g',  netcdf_int, ('number_of_volumes',))
+        outfile.createVariable('node_l2g', netcdf_int, ('number_of_points',))
+        outfile.createVariable('tri_full_flag', netcdf_int, ('number_of_volumes',))
+
+        print tri_l2g.shape
+        print tri_l2g
+        print outfile.variables['tri_l2g'].shape
+
+        outfile.variables['tri_l2g'][:] = tri_l2g.astype(num.int32)
+
+        print node_l2g.shape
+        print node_l2g
+        print outfile.variables['node_l2g'].shape
+
+        outfile.variables['node_l2g'][:] = node_l2g.astype(num.int32)
+
+        print tri_full_flag.shape
+        print tri_full_flag
+        print outfile.variables['tri_full_flag'].shape
+
+        outfile.variables['tri_full_flag'][:] = tri_full_flag.astype(num.int32)
 
 

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

@@ -14 +14 @@
 
     output = domain_global_name+".sww"
-    swwfiles = [ domain_global_name+"_P"+str(v)+"_"+str(np)+".sww" for v in range(np)]
+    swwfiles = [ domain_global_name+"_P"+str(np)+"_"+str(v)+".sww" for v in range(np)]
 
     _sww_merge(swwfiles, output, verbose)
@@ -132 +132 @@
 
 
+
+    from anuga.coordinate_transforms.geo_reference import Geo_reference
+    geo_reference = Geo_reference()
     
-    
-    sww.store_triangulation(fido, points, out_tris)
+    sww.store_triangulation(fido, points, out_tris, points_georeference=geo_reference)
 
     fido.order      = order

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

@@ -7 +7 @@
 from anuga.utilities.file_utils import copy_code_files, get_all_swwfiles
 from anuga.utilities.file_utils import del_dir
-from sww_merge import sww_merge
+from sww_merge import sww_merge, _sww_merge
 
 
@@ -139 +139 @@
                 
         outfile = 'test_out.sww'
-        sww_merge(['test1.sww', 'test2.sww'], outfile)
+        _sww_merge(['test1.sww', 'test2.sww'], outfile)
         self.failUnless(os.access(outfile, os.F_OK))  
         

trunk/anuga_core/source/anuga_parallel/distribute_mesh.py

@@ -147 +147 @@
         del npart
 
-        # Sometimes (usu. on x86_64), partMeshNodal returnes an array of zero
+        # Sometimes (usu. on x86_64), partMeshNodal returns an array of zero
         # dimensional arrays. Correct this.
         if type(epart[0]) == num.ndarray:
@@ -808 +808 @@
 #########################################################
 
-def build_local_GA(nodes, triangles, boundaries, tri_index):
+def build_local_GA(nodes, triangles, boundaries, tri_map):
 
     Nnodes =len(nodes)
@@ -823 +823 @@
         if nodes[i][0] > NGlobal:
             NGlobal = nodes[i][0]
-    index = num.zeros(int(NGlobal)+1, num.int)
-    num.put(index, num.take(nodes, (0,), 1).astype(num.int), \
+
+    node_map = -1*num.ones(int(NGlobal)+1, num.int)
+
+    num.put(node_map, num.take(nodes, (0,), 1).astype(num.int), \
         num.arange(Nnodes))
         
@@ -830 +832 @@
 
     GAtriangles = num.zeros((Ntriangles, 3), num.int)
-    GAtriangles[:,0] = num.take(index, triangles[:,0])
-    GAtriangles[:,1] = num.take(index, triangles[:,1])
-    GAtriangles[:,2] = num.take(index, triangles[:,2])
+    GAtriangles[:,0] = num.take(node_map, triangles[:,0])
+    GAtriangles[:,1] = num.take(node_map, triangles[:,1])
+    GAtriangles[:,2] = num.take(node_map, triangles[:,2])
 
     # Change the triangle numbering in the boundaries
@@ -838 +840 @@
     GAboundaries = {}
     for b in boundaries:
-        GAboundaries[tri_index[b[0]], b[1]] = boundaries[b]
+        GAboundaries[tri_map[b[0]], b[1]] = boundaries[b]
         
-    del (index)
-    
-    return GAnodes, GAtriangles, GAboundaries
+    
+    return GAnodes, GAtriangles, GAboundaries, node_map
 
 
@@ -868 +869 @@
 #########################################################
 
-def build_local_commun(index, ghostc, fullc, nproc):
+def build_local_commun(tri_map, ghostc, fullc, nproc):
 
     # Initialise
@@ -886 +887 @@
             ghost_recv[c] = [0, 0]
             ghost_recv[c][1] = d
-            ghost_recv[c][0] = num.take(index, d)
+            ghost_recv[c][0] = num.take(tri_map, d)
             
     # Build a temporary copy of the full_send dictionary
@@ -899 +900 @@
                 tmp_send[neigh] = []
             tmp_send[neigh].append([global_id, \
-                                    index[global_id]])
+                                    tri_map[global_id]])
 
     # Extract the full send information and put it in the form
@@ -958 +959 @@
         if id > NGlobal:
             NGlobal = id
-    index = num.zeros(int(NGlobal)+1, num.int)
-    index[lower_t:upper_t]=num.arange(upper_t-lower_t)
+    #index = num.zeros(int(NGlobal)+1, num.int)
+    tri_map = -1*num.ones(int(NGlobal)+1, num.int)
+    tri_map[lower_t:upper_t]=num.arange(upper_t-lower_t)
     for i in range(len(submesh["ghost_triangles"])):
-        index[submesh["ghost_triangles"][i][0]] = i+upper_t-lower_t
+        tri_map[submesh["ghost_triangles"][i][0]] = i+upper_t-lower_t
     
     # Change the node numbering (and update the numbering in the
     # triangles)
 
-    [GAnodes, GAtriangles, GAboundary] = build_local_GA(nodes, triangles, boundaries, index)
+    [GAnodes, GAtriangles, GAboundary, node_map] = \
+    build_local_GA(nodes, triangles, boundaries, tri_map)
 
     # Extract the local quantities
@@ -984 +987 @@
     fcommun = submesh["full_commun"]
     [ghost_rec, full_send] = \
-                build_local_commun(index, gcommun, fcommun, nproc)
-
-    # Clean up before exiting
-
-    del(index)
+                build_local_commun(tri_map, gcommun, fcommun, nproc)
+
 
     return GAnodes, GAtriangles, GAboundary, quantities, ghost_rec, \
-           full_send
+           full_send, tri_map, node_map
 
 
@@ -1284 +1284 @@
     # datastructure
 
-    [GAnodes, GAtriangles, boundary, quantities, ghost_rec, full_send] = \
+    [GAnodes, GAtriangles, boundary, quantities, \
+     ghost_rec, full_send, tri_map, node_map] = \
               build_local_mesh(submesh_cell, lower_t, upper_t, \
                                numproc)
@@ -1290 +1291 @@
     return GAnodes, GAtriangles, boundary, quantities,\
            ghost_rec, full_send,\
-           number_of_full_nodes, number_of_full_triangles
+           number_of_full_nodes, number_of_full_triangles, tri_map, node_map
+          
 
 
@@ -1327 +1329 @@
 
     numprocs = len(triangles_per_proc)
-    points, vertices, boundary, quantities, ghost_recv_dict, full_send_dict = \
+    points, vertices, boundary, quantities, ghost_recv_dict, \
+            full_send_dict, tri_map, node_map = \
             build_local_mesh(submesh_cell, 0, triangles_per_proc[0], numprocs)
 
 
     return  points, vertices, boundary, quantities, ghost_recv_dict, \
-           full_send_dict
-
-
-
-
+           full_send_dict, tri_map, node_map
+           
+
+
+
+

trunk/anuga_core/source/anuga_parallel/parallel_api.py

@@ -99 +99 @@
                 ghost_recv_dict, full_send_dict,\
                 number_of_full_nodes, number_of_full_triangles,\
-                s2p_map, p2s_map =\
+                s2p_map, p2s_map, tri_map, node_map =\
                 distribute_mesh(domain, verbose=verbose)
+
+        number_of_global_triangles = len(tri_map)
+        number_of_global_nodes = len(node_map)
+
+        # Extract l2g maps
+        tri_l2g  = extract_l2g_map(tri_map)
+        node_l2g = extract_l2g_map(node_map)
 
         # Send serial to parallel (s2p) and parallel to serial (p2s) triangle mapping to proc 1 .. numprocs
@@ -115 +122 @@
         points, vertices, boundary, quantities,\
                 ghost_recv_dict, full_send_dict,\
-                number_of_full_nodes, number_of_full_triangles =\
+                number_of_full_nodes, number_of_full_triangles, \
+                tri_map, node_map =\
                 rec_submesh(0, verbose)
 
+        number_of_global_triangles = len(tri_map)
+        number_of_global_nodes = len(node_map)
+
+        # Extract l2g maps
+        tri_l2g  = extract_l2g_map(tri_map)
+        node_l2g = extract_l2g_map(node_map)
+        
         # Recieve serial to parallel (s2p) and parallel to serial (p2s) triangle mapping
         s2p_map = receive(0)
@@ -136 +151 @@
                              number_of_full_triangles=number_of_full_triangles,
                              geo_reference=georef,
-                             tri_map = s2p_map,
-                             inv_tri_map = p2s_map) ## jj added this
+                             number_of_global_triangles = number_of_global_triangles,
+                             number_of_global_nodes = number_of_global_nodes,
+                             s2p_map = s2p_map,
+                             p2s_map = p2s_map, ## jj added this
+                             tri_l2g = tri_l2g, ## SR added this
+                             node_l2g = node_l2g)
 
     #------------------------------------------------------------------------
@@ -177 +196 @@
     # Subdivide the mesh
     if verbose: print 'Subdivide mesh'
-    nodes, triangles, boundary, triangles_per_proc, quantities, s2p_map, p2s_map = \
+    nodes, triangles, boundary, triangles_per_proc, quantities, \
+           s2p_map, p2s_map = \
            pmesh_divide_metis_with_map(domain, numprocs)
 
     #PETE: s2p_map (maps serial domain triangles to parallel domain triangles)
-    #p2_map (maps parallel domain triangles to domain triangles)
+    #      sp2_map (maps parallel domain triangles to domain triangles)
 
 
@@ -204 +224 @@
 
     # Build the local mesh for processor 0
-    points, vertices, boundary, quantities, ghost_recv_dict, full_send_dict =\
+    points, vertices, boundary, quantities, \
+            ghost_recv_dict, full_send_dict, tri_map, node_map =\
               extract_hostmesh(submesh, triangles_per_proc)
 
@@ -240 +261 @@
     return points, vertices, boundary, quantities,\
            ghost_recv_dict, full_send_dict,\
-           number_of_full_nodes, number_of_full_triangles, s2p_map, p2s_map
-    
-
-
-
+           number_of_full_nodes, number_of_full_triangles, \
+           s2p_map, p2s_map, tri_map, node_map
+    
+
+
+def extract_l2g_map(map):
+    # Extract l2g_map
+
+    import numpy as num
+    
+    b = num.arange(len(map))
+
+    l_ids = num.extract(map>-1,map)
+    g_ids = num.extract(map>-1,b)
+
+#    print len(g_ids)
+#    print len(l_ids)
+#    print l_ids
+
+    l2g = num.zeros_like(g_ids)
+    l2g[l_ids] = g_ids
+
+    return l2g
+

trunk/anuga_core/source/anuga_parallel/parallel_shallow_water.py

@@ -33 +33 @@
                  number_of_full_triangles=None,
                  geo_reference=None,
-                 tri_map=None,
-                 inv_tri_map=None): #jj added this
+                 number_of_global_triangles=None, ## SR added this
+                 number_of_global_nodes= None, ## SR added this
+                 s2p_map=None,
+                 p2s_map=None, #jj added this
+                 tri_l2g = None, ## SR added this
+                 node_l2g = None): ## SR added this
 
         Domain.__init__(self,
@@ -48 +52 @@
                         geo_reference=geo_reference) #jj added this
         
+
+
+        self.parallel = True
+
         # PETE: Find the number of full nodes and full triangles, this is a temporary fix
         # until the bug with get_number_of_full_[nodes|triangles]() is fixed.
@@ -64 +72 @@
 
         self.global_name = 'domain'
-        self.tri_map = tri_map
-        self.inv_tri_map = inv_tri_map
+
+        self.number_of_global_triangles=number_of_global_triangles
+        self.number_of_global_nodes = number_of_global_nodes
+
+        self.s2p_map = s2p_map
+        self.p2s_map = p2s_map
+        self.tri_l2g = tri_l2g
+        self.node_l2g = node_l2g
 
 
@@ -78 +92 @@
 
         # Call parents method with processor number attached.
-        Domain.set_name(self, name + '_P%d_%d' %(self.processor, self.numproc))
+        Domain.set_name(self, name + '_P%d_%d' %(self.numproc, self.processor))
 
 

trunk/anuga_core/source/anuga_parallel/run_parallel_sw_merimbula.py

@@ -35 +35 @@
 
 
-from anuga_parallel import distribute, myid, numprocs, finalize
+from anuga_parallel import distribute, myid, numprocs, finalize, barrier
 
 
@@ -42 +42 @@
 #--------------------------------------------------------------------------
 
-mesh_filename = "merimbula_10785_1.tsh" ; x0 = 756000.0 ; x1 = 756500.0
+#mesh_filename = "merimbula_10785_1.tsh" ; x0 = 756000.0 ; x1 = 756500.0
 #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-100.tsh" ; x0 = 0.25 ; x1 = 0.5
 yieldstep = 20
 finaltime = 500
@@ -70 +70 @@
     domain = create_domain_from_file(mesh_filename)
     domain.set_quantity('stage', Set_Stage(x0, x1, 2.0))
+
+    print domain.statistics()
+    print domain.get_extent()
+    print domain.get_extent(absolute=True)
+    print domain.geo_reference
 else:
     domain = None
@@ -86 +91 @@
 #domain.set_CFL(0.7)
 #domain.set_beta(1.5)
-domain.set_name('meribula')
+domain.set_name('merimbula')
+
+for p in range(numprocs):
+    if myid == p:
+        print 'P%d'%p
+        print domain.get_extent()
+        print domain.get_extent(absolute=True)
+        print domain.geo_reference
+        #print domain.tri_map
+        #print domain.inv_tri_map
+        print domain.tri_l2g
+        print domain.node_l2g
+    else:
+        pass
+    
+    barrier()
+
 
 #------------------------------------------------------------------------------

trunk/anuga_core/source/anuga_parallel/test_parallel_distribute_mesh.py

@@ -195 +195 @@
         # Test extract_hostmesh
         #----------------------------------------------------------------------------------
-        points, vertices, boundary, quantities, ghost_recv_dict, full_send_dict  =\
+        points, vertices, boundary, quantities, \
+                ghost_recv_dict, full_send_dict, tri_map, node_map  =\
         extract_hostmesh(submesh, triangles_per_proc)
 
@@ -222 +223 @@
         # Test rec_submesh
         #----------------------------------------------------------------------------------
-        points, vertices, boundary, quantities, ghost_recv_dict, full_send_dict, no_full_nodes, no_full_trigs = rec_submesh(0, verbose=False)    
+        points, vertices, boundary, quantities, \
+                ghost_recv_dict, full_send_dict, \
+                no_full_nodes, no_full_trigs, tri_map, node_map  = \
+                rec_submesh(0, verbose=False)    
 
         if myid == 1:
@@ -235 +239 @@
             true_full_send =  {0: [num.array([0, 1, 3, 4, 5]), num.array([ 5,  6,  8,  9, 10])], 2: [num.array([0, 1]), num.array([5, 6])]}
 
+            true_tri_map  = num.array([ 6,  7,  8, -1,  9,  0, 1,  2,  3,  4,  5, 10, 11])
+
+            true_node_map = num.array([ 0,  1,  2,  7,  3,  4, -1,  8,  9,  5,  6, 10, 11])
+
+            assert_(num.allclose(tri_map,   true_tri_map))
+            assert_(num.allclose(node_map,   true_node_map))
             assert_(num.allclose(points,   true_points))
             assert_(num.allclose(vertices, true_vertices))
@@ -249 +259 @@
             true_vertices =  [[1, 5, 0], [1, 6, 2], [3, 6, 1], [4, 6, 3], [2, 6, 4], [2, 5, 1], [2, 10, 8], [4, 10, 2], [9, 10, 4], [0, 5, 7], [7, 5, 2]]
 
-
             true_ghost_recv =   {0: [num.array([5, 6, 7, 8]), num.array([0, 1, 2, 3])], 1: [num.array([ 9, 10]), num.array([5, 6])]}
 
             true_full_send =   {0: [num.array([0, 1, 2, 3, 4]), num.array([11, 12, 13, 14, 15])], 1: [num.array([0, 1]), num.array([11, 12])]}
 
-
+            true_tri_map  = num.array([5, 6, 7, 8, -1, 9, 10, -1, -1, -1, -1, 0, 1, 2, 3, 4, -1])
+
+            true_node_map = num.array([ 0,  7, -1,  1,  2,  8 , 3,  4,  9,  5, -1,  6, 10])
+
+            assert_(num.allclose(tri_map,   true_tri_map))
+            assert_(num.allclose(node_map,   true_node_map))
             assert_(num.allclose(points,   true_points))
             assert_(num.allclose(vertices, true_vertices))