Changeset 2625

Timestamp:

Mar 29, 2006, 9:41:58 AM (19 years ago)

Author:

linda

Message:

Added ghost boundary layers

Location:

inundation/parallel

Files:

• build_commun.py (modified) (8 diffs)
• build_local.py (modified) (4 diffs)
• build_submesh.py (modified) (10 diffs)
• parallel_meshes.py (modified) (7 diffs)
• run_parallel_advection.py (modified) (3 diffs)
• run_parallel_merimbula.py (modified) (5 diffs)
• run_parallel_sw_rectangle.py (modified) (2 diffs)

inundation/parallel/build_commun.py

@@ -59 +59 @@
              tagmap[bkey] = counter
              counter = counter+1
+    for b in submesh["ghost_boundary"][p]:
+         bkey = submesh["ghost_boundary"][p][b]
+         if not tagmap.has_key(bkey):
+             tagmap[bkey] = counter
+             counter = counter+1
     pypar.send(tagmap, p)
 
@@ -79 +84 @@
     # send the array sizes so memory can be allocated
 
-    setup_array = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
+    setup_array = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
     setup_array[0] = len(submesh["full_nodes"][p])
     setup_array[1] = len(submesh["ghost_nodes"][p])
@@ -85 +90 @@
     setup_array[3] = len(submesh["ghost_triangles"][p])
     setup_array[4] = len(submesh["full_boundary"][p])
-    setup_array[5] = len(submesh["ghost_commun"][p])
-    setup_array[6] = len(flat_full_commun)
+    setup_array[5] = len(submesh["ghost_boundary"][p])
+    setup_array[6] = len(submesh["ghost_commun"][p])
+    setup_array[7] = len(flat_full_commun)
 
     pypar.send(setup_array, p)
@@ -105 +111 @@
     for b in submesh["full_boundary"][p]:
         bc.append([b[0], b[1], tagmap[submesh["full_boundary"][p][b]]])
+    pypar.send(bc, p, use_buffer=True)
+    bc = []
+    for b in submesh["ghost_boundary"][p]:
+        bc.append([b[0], b[1], tagmap[submesh["ghost_boundary"][p][b]]])
     pypar.send(bc, p, use_buffer=True)
 
@@ -150 +160 @@
     
     tagmap = pypar.receive(p)
+    itagmap = {}
+    for t in tagmap:
+        itagmap[tagmap[t]]=t
 
     # receive the quantities key information
@@ -165 +178 @@
     # recieve information about the array sizes
 
-    setup_array = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
+    setup_array = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
     setup_array = pypar.receive(p, setup_array)
 
@@ -201 +214 @@
     bnd_c = pypar.receive(p, bc)
 
-    itagmap = {}
-    for t in tagmap:
-        itagmap[tagmap[t]]=t
-
     submesh_cell["full_boundary"] = {}
     for b in bnd_c:
         submesh_cell["full_boundary"][b[0],b[1]]=itagmap[b[2]]
 
+    # receive the ghost boundary
+
+    no_ghost_boundary = setup_array[5]
+    bc = []
+    for i in range(no_ghost_boundary):
+        bc.append([0.0, 0.0, 0.0])
+    bnd_c = pypar.receive(p, bc)
+
+    submesh_cell["ghost_boundary"] = {}
+    for b in bnd_c:
+        submesh_cell["ghost_boundary"][b[0],b[1]]=itagmap[b[2]]
+
     # receive the ghost communication pattern
 
-    no_ghost_commun = setup_array[5]
+    no_ghost_commun = setup_array[6]
     ghost_commun = zeros((no_ghost_commun, 2), Int)
     submesh_cell["ghost_commun"] = pypar.receive(p, ghost_commun)
@@ -217 +238 @@
     # receive the full communication pattern
 
-    no_full_commun = setup_array[6]
+    no_full_commun = setup_array[7]
     full_commun = []
     for i in range(no_full_commun):

inundation/parallel/build_local.py

@@ -40 +40 @@
 #########################################################
 
-def build_local_GA(nodes, triangles):
+def build_local_GA(nodes, triangles, boundaries, tri_index):
 
     Nnodes =len(nodes)
@@ -65 +65 @@
     GAtriangles[:,1] = take(index, triangles[:,1])
     GAtriangles[:,2] = take(index, triangles[:,2])
-    
+
+    # Change the triangle numbering in the boundaries
+
+    GAboundaries = {}
+    for b in boundaries:
+        GAboundaries[tri_index[b[0]], b[1]] = boundaries[b]
+        
     del (index)
-
-    return GAnodes, GAtriangles
+    
+    return GAnodes, GAtriangles, GAboundaries
 
 
@@ -171 +177 @@
     triangles = concatenate((submesh["full_triangles"], gtri))
 
-    # Renumber the boundary edges to correspond to the new
-    # triangle numbering
-
-    GAboundary = {}
-    for b in submesh["full_boundary"]:
-        GAboundary[b[0]-lower_t,b[1]] = submesh["full_boundary"][b]
+    # Combine the full boundaries and ghost boundaries
+
+    boundaries = submesh["full_boundary"]
+    for b in submesh["ghost_boundary"]:
+        boundaries[b]=submesh["ghost_boundary"][b]
 
     # Make note of the new triangle numbers, including the ghost
@@ -190 +195 @@
     for i in range(len(submesh["ghost_triangles"])):
         index[submesh["ghost_triangles"][i][0]] = i+upper_t-lower_t
-
+    
     # Change the node numbering (and update the numbering in the
     # triangles)
 
-    [GAnodes, GAtriangles] = build_local_GA(nodes, triangles)
+    [GAnodes, GAtriangles, GAboundary] = build_local_GA(nodes, triangles, boundaries, index)
 
     # Extract the local quantities

inundation/parallel/build_submesh.py

@@ -113 +113 @@
 # be added to the node list for the current processor
 #
-#  *) The boundary edges for the ghost triangles are
-# ignored.
 #
 # -------------------------------------------------------
@@ -206 +204 @@
     return subnodes, subtriangles
 
+#########################################################
+#
+# Find the edges of the ghost trianlges that do not
+# have a neighbour in the current cell. These are
+# treated as a special type of boundary edge.
+#
+#  *) Given the ghost triangles in a particular
+# triangle, use the mesh to find its neigbours. If
+# the neighbour is not in the processor set it to
+# be a boundary edge
+#
+#  *) The vertices in the ghost triangles must also
+# be added to the node list for the current processor
+#
+#  *) The boundary edges for the ghost triangles are
+# ignored.
+#
+# -------------------------------------------------------
+#
+#  *) The type assigned to the ghost boundary edges is 'ghost'
+#
+#  *)  The boundary information is returned as a directorier
+# with the key = (triangle id, edge no) and the values
+# assigned to the key is 'ghost'
+#
+#
+#########################################################
+def is_in_processor(ghost_list, tlower, tupper, n):
+    return (n in ghost_list) or (tlower <= n and tupper >= n)
+
+def ghost_bnd_layer(ghosttri, tlower, tupper, mesh, p):
+
+    ghost_list = []
+    subboundary = {}
+        
+    for t in ghosttri:
+        ghost_list.append(t[0])
+
+    for t in ghosttri:
+        n = mesh.neighbours[t[0], 0]
+        if not is_in_processor(ghost_list, tlower, tupper, n):
+            subboundary[t[0], 0] = 'ghost'
+
+        n = mesh.neighbours[t[0], 1]
+        if not is_in_processor(ghost_list, tlower, tupper, n):
+            subboundary[t[0], 1] = 'ghost'
+
+        n = mesh.neighbours[t[0], 2]
+        if not is_in_processor(ghost_list, tlower, tupper, n):
+            subboundary[t[0], 2] = 'ghost'
+    
+    return subboundary
 
 #########################################################
@@ -318 +368 @@
 # updates must also be built.
 #
+#  *) Assumes that full triangles, nodes etc have already
+# been found and stored in submesh
+#
 #  *) See the documentation for ghost_layer,
 # ghost_commun_pattern and full_commun_pattern
@@ -327 +380 @@
 # ghost_commun_pattern and full_commun_pattern
 #
-#  *) The ghost_triangles, ghost_nodes, ghost_commun and
-# full_commun is added to submesh
+#  *) The ghost_triangles, ghost_nodes, ghost_boundary,
+# ghost_commun and full_commun is added to submesh
 #########################################################
 
@@ -338 +391 @@
     ghost_nodes = []
     ghost_commun = []
+    ghost_bnd = []
 
     # Loop over the processors
@@ -354 +408 @@
         ghost_nodes.append(subnodes)
 
+        # Find the boundary layer formed by the ghost triangles
+        
+        subbnd = ghost_bnd_layer(subtri, tupper, tlower, mesh, p)
+        ghost_bnd.append(subbnd)
+        
         # Build the communication pattern for the ghost nodes
 
@@ -369 +428 @@
     submesh["ghost_triangles"] = ghost_triangles
     submesh["ghost_commun"] = ghost_commun
-
+    submesh["ghost_boundary"] = ghost_bnd
+    
     # Build the communication pattern for the full triangles
 
@@ -449 +509 @@
 #  *) A dictionary containing the full_triangles,
 # full_nodes, full_boundary, ghost_triangles, ghost_nodes,
-# ghost_commun and full_commun is returned.
+# ghost_boundary, ghost_commun and full_commun is returned.
 #
 #########################################################
@@ -489 +549 @@
 #  *) A dictionary containing the full_triangles,
 # full_nodes, full_boundary, ghost_triangles, ghost_nodes,
-# ghost_commun and full_commun belonging to processor zero
-# are returned.
+# ghost_boundary, ghost_commun and full_commun belonging
+# to processor zero are returned.
 #
 #########################################################
@@ -501 +561 @@
     submesh_cell["ghost_triangles"] = submesh["ghost_triangles"][0]
     submesh_cell["full_boundary"] = submesh["full_boundary"][0]
+    submesh_cell["ghost_boundary"] = submesh["ghost_boundary"][0]
     submesh_cell["ghost_commun"] = submesh["ghost_commun"][0]
     submesh_cell["full_commun"] = submesh["full_commun"][0]

inundation/parallel/parallel_meshes.py

@@ -12 +12 @@
 Ole Nielsen, Stephen Roberts, Duncan Gray, Christopher Zoppou
 Geoscience Australia, 2005
+
+Modified by Linda Stals, March 2006, to include ghost boundaries
+
 """
 
@@ -124 +127 @@
 
             if i == m-1:
-                boundary[nt, 2] = 'right'
+                if processor == numproc-1:
+                    boundary[nt, 2] = 'right'
+                else:
+                    boundary[nt, 2] = 'ghost'
+        
             if j == 0:
                 boundary[nt, 1] = 'bottom'
@@ -144 +151 @@
 
             if i == 0:
-                boundary[nt, 2] = 'left'
+                if processor == 0:
+                    boundary[nt, 2] = 'left'
+                else:
+                    boundary[nt, 2] = 'ghost'
             if j == n-1:
                 boundary[nt, 1] = 'top'
@@ -263 +273 @@
 
             if i == m-1:
-                boundary[nt, 2] = 'right'
+                if processor == numproc-1:
+                    boundary[nt, 2] = 'right'
+                else:
+                    boundary[nt, 2] = 'ghost'
             if j == 0:
                 boundary[nt, 1] = 'bottom'
@@ -282 +295 @@
 
             if i == 0:
-                boundary[nt, 2] = 'left'
+                if processor == 0:
+                    boundary[nt, 2] = 'left'
+                else:
+                    boundary[nt, 2] = 'ghost'
             if j == n-1:
                 boundary[nt, 1] = 'top'
@@ -391 +407 @@
 
             if i == m-1:
-                boundary[nt, 2] = 'right'
+                if processor == numproc-1:
+                    boundary[nt, 2] = 'right'
+                else:
+                    boundary[nt, 2] = 'ghost'
             if j == 0:
                 boundary[nt, 1] = 'bottom'
@@ -404 +423 @@
 
             if i == 0:
-                boundary[nt, 2] = 'left'
+                if processor == 0:
+                    boundary[nt, 2] = 'left'
+                else:
+                    boundary[nt, 2] = 'ghost'
             if j == n-1:
                 boundary[nt, 1] = 'top'

inundation/parallel/run_parallel_advection.py

@@ -5 +5 @@
 sys.path.append('..'+sep+'pyvolution')
 
+
 #========================================================================
+
+import pypar
+
 from config import g, epsilon
 from Numeric import allclose, array, zeros, ones, Float
@@ -14 +18 @@
 from parallel_advection import Transmissive_boundary, Dirichlet_boundary
 
-import pypar
-
 numprocs = pypar.size()
 myid = pypar.rank()
 processor_name = pypar.Get_processor_name()
-
 
 
@@ -41 +42 @@
 #Boundaries
 T = Transmissive_boundary(domain)
-D = Dirichlet_boundary(array([1.0]))
+D = Dirichlet_boundary(array([1.0])) 
 
 domain.default_order = 2
 
-
-domain.set_boundary( {'left': T, 'right': T, 'bottom': T, 'top': T} )
+domain.set_boundary( {'left': T, 'right': T, 'bottom': T, 'top': T, 'ghost': None} )
 domain.check_integrity()
 

inundation/parallel/run_parallel_merimbula.py

@@ -87 +87 @@
     # read in the test files
 
-    filename = 'test-100.tsh'
-#    filename = 'merimbula_10785.tsh'
+#    filename = 'test-100.tsh'
+    filename = 'merimbula_10785.tsh'
     nx = numprocs
     ny = 1
@@ -95 +95 @@
 
     domain_full = pmesh_to_domain_instance(filename, Advection_Domain)
-    domain_full.set_quantity('stage', Set_Stage(200.0,300.0,1.0))
-#    domain.set_quantity('stage', Set_Stage(756000.0,756500.0,4.0))
+#    domain_full.set_quantity('stage', Set_Stage(200.0,300.0,1.0))
+    domain_full.set_quantity('stage', Set_Stage(756000.0,756500.0,4.0))
 
     nodes, triangles, boundary, triangles_per_proc, quantities  =\
             pmesh_divide(domain_full, nx, ny)
-
+    
     # subdivide the mesh
     rect = array(domain_full.xy_extent, Float)
-
-    print rect
 
     rect = array(rect, Float)
@@ -128 +126 @@
              rec_submesh(0)
 
-#if myid == 0:
-#    print 'ghost'
-#    print ghost_recv_dict
-#
-#if myid == 0:
-#    print 'full'
-#    print full_send_dict
-
-
 pypar.broadcast(rect,0)
-print rect
 
 domain = Parallel_Domain(points, vertices, boundary,
@@ -150 +138 @@
 
 T = Transmissive_boundary(domain)
-R = Reflective_boundary(domain)
-domain.set_boundary( {'outflow': R, 'inflow': R, 'inner':R, 'exterior': R, 'open':R} )
+#R = Reflective_boundary(domain)
+domain.set_boundary( {'outflow': T, 'inflow': T, 'inner':T, 'exterior': T, 'open':T, 'ghost':None} )
 
 
@@ -170 +158 @@
 
 
+from norms import linf_norm
 
 #yieldstep = 1
+#finaltime = 4000
+
 yieldstep = 1000
 finaltime = 50000
-# finaltime = 4000
+
 for t in domain.evolve(yieldstep = yieldstep, finaltime = finaltime):
     if myid == 0:
         domain.write_time()
+##    Stage = domain.quantities['stage']
+##    max_stage = max(max(Stage.edge_values))
+##    print max_stage
+##    if max_stage > 4.0:
+##        print "######################", max_stage
 
 if myid == 0:

inundation/parallel/run_parallel_sw_rectangle.py

@@ -46 +46 @@
 N = M*numprocs
 
-#N = M = 250
-
 if myid == 0:
     print 'N == %d' %N
 
 points, vertices, boundary, full_send_dict, ghost_recv_dict =  \
-    parallel_rectangle(N, M, len1_g=1.0*numprocs, len2_g=1.0)
+    parallel_rectangle(N, M, len1_g=1.0*numprocs, len2_g = 1.0)
 
 
@@ -84 +82 @@
 
 
-domain.set_boundary( {'left': R, 'right': R, 'bottom': R, 'top': R} )
+domain.set_boundary( {'left': R, 'right': R, 'bottom': R, 'top': R, 'ghost': None} )
 domain.check_integrity()