Changeset 1424

Timestamp:

May 18, 2005, 5:36:32 PM (20 years ago)

Author:

steve

Message:

 

Location:

inundation/ga/storm_surge

Files:

• parallel/advection.py (modified) (1 diff)
• parallel/parallel_advection.py (modified) (13 diffs)
• parallel/run_parallel_advection.py (modified) (1 diff)
• zeus/pyvolution.zbi (modified) (previous)

inundation/ga/storm_surge/parallel/advection.py

@@ -33 +33 @@
 
 
-        Generic_domain.__init__(self, coordinates, vertices, boundary,
-                                ['stage'])
+        Generic_domain.__init__(self, coordinates, vertices, boundary,['stage'])
 
         if velocity is None:

inundation/ga/storm_surge/parallel/parallel_advection.py

@@ -16 +16 @@
 from advection import *
 Advection_Domain = Domain
-from Numeric import zeros, Float, Int, ones, allclose,  array
+from Numeric import zeros, Float, Int, ones, allclose, array
 import pypar
 
@@ -22 +22 @@
 class Parallel_Domain(Advection_Domain):
 
-    def __init__(self, coordinates, vertices, boundary = None, ghosts = None,
+    def __init__(self, coordinates, vertices, boundary = None,
+                 full_send_dict = None, ghost_recv_dict = None,
                  velocity = None):
-
-        Advection_Domain.__init__(self, coordinates, vertices, boundary, velocity)
-
-        N = self.number_of_elements
 
         self.processor = pypar.rank()
         self.numproc   = pypar.size()
-
-        if  ghosts == None:
-            self.ghosts = None
-            self.full = ones( N, Float)
-        else:
-            self.ghosts = ghosts
-            self.full = ones( N, Int)
-            keys = self.ghosts.keys()
-            for key in keys:
-                self.full[key] = -1
-
-
-
+        print 'Processor %d'%self.processor
+        velocity = [(self.processor+1),0.0]
+
+        print 'velocity',velocity
+
+        Advection_Domain.__init__(self, coordinates, vertices, boundary, velocity)
+
+        N = self.number_of_elements
+
+        self.processor = pypar.rank()
+        self.numproc   = pypar.size()
+
+        self.full_send_dict  = full_send_dict
+        self.ghost_recv_dict = ghost_recv_dict
+
+        #print self.full_send_dict
+        #print self.ghost_recv_dict
 
     def check_integrity(self):
@@ -67 +68 @@
         pypar.raw_reduce(ltimestep, gtimestep, pypar.MIN, 0)
         pypar.broadcast(gtimestep,0)
-        pypar.Barrier()
+        #pypar.Barrier()
 
         self.timestep = gtimestep[0]
@@ -82 +83 @@
         stage_cv = self.quantities['stage'].centroid_values
 
-#        for iproc in full_send_dict:
-#            Ids = full_send_dict[iproc][0]
-#            X   = full_send_dict[iproc][1]
-#            for i _ in enumerate(X):
-#                X[i] = stage_cv[Ids[i]]
-#            pypar.send(X,iproc)
+        # update of non-local ghost cells
+#        for iproc in range(self.numproc):
 #
-#        for iproc in ghost_recv_dict:
-#            Ids = ghost_recv_dict[iproc][0]
-#            X   = ghost_recv_dict[iproc][1]
-#            X = pypar.receive(iproc,X)
-#            for i _ in enumerate(X):
-#                stage_cv[Ids[i]] = X[i]
-
-
-        if self.ghosts is not None:
-            stage_cv = self.quantities['stage'].centroid_values
-            for triangle in self.ghosts:
-                stage_cv[triangle] = stage_cv[self.ghosts[triangle]]
+#            if iproc == self.processor:
+#                #Send data from iproc processor to other processors
+#                for send_proc in full_send_dict:
+#                    if send_proc != iproc:
+#                        Idf  = full_send_dict[iproc][0]
+#                        Xout = full_send_dict[iproc][1]
+#                        for i _ in enumerate(X):
+#                            Xout[i] = stage_cv[Idf[i]]
+#                        pypar.send(Xout,send_proc)
+#            else:
+#                #Receive data from the iproc processor
+#                if  ghost_recv_dict.has_key(iproc):
+#                    Idg = ghost_recv_dict[iproc][0]
+#                    X   = ghost_recv_dict[iproc][1]
+#                    X = pypar.receive(iproc,X)
+#                    for i _ in enumerate(X):
+#                        stage_cv[Idg[i]] = X[i]
+
+        #local update of ghost cells
+        iproc = self.processor
+        if self.full_send_dict.has_key(iproc):
+            Idf  = self.full_send_dict[iproc][0]
+            #print Idf
+            Idg  = self.ghost_recv_dict[iproc][0]
+            #print Idg
+            for i, _ in enumerate(Idf):
+                #print i,Idg[i],Idf[i]
+                stage_cv[Idg[i]] = stage_cv[Idf[i]]
+
+
+#        if self.ghosts is not None:
+#            stage_cv = self.quantities['stage'].centroid_values
+#            for triangle in self.ghosts:
+#                stage_cv[triangle] = stage_cv[self.ghosts[triangle]]
 
 
@@ -137 +156 @@
 
 
-def Parallel_rectangular(m, n, len1=1.0, len2=1.0, origin = (0.0, 0.0)):
+def parallel_rectangular(m, n, len1=1.0, len2=1.0, origin = (0.0, 0.0)):
 
 
@@ -149 +168 @@
     len2: y direction (bottom to top)
 
+    """
+
+    from config import epsilon
+    from Numeric import zeros, Float, Int
+
+    processor = pypar.rank()
+    numproc   = pypar.size()
+
+
+
+    delta1 = float(len1)/m
+    delta2 = float(len2)/n
+
+    #Calculate number of points
+    Np = (m+1)*(n+1)
+
+    class VIndex:
+
+        def __init__(self, n,m):
+            self.n = n
+            self.m = m
+
+        def __call__(self, i,j):
+            return j+i*(self.n+1)
+
+    class EIndex:
+
+        def __init__(self, n,m):
+            self.n = n
+            self.m = m
+
+        def __call__(self, i,j):
+            return 2*(j+i*self.n)
+
+
+    I = VIndex(n,m)
+    E = EIndex(n,m)
+
+    points = zeros( (Np,2), Float)
+
+    for i in range(m+1):
+        for j in range(n+1):
+
+            points[I(i,j),:] = [i*delta1 + origin[0], j*delta2 + origin[1]]
+
+    #Construct 2 triangles per rectangular element and assign tags to boundary
+    #Calculate number of triangles
+    Nt = 2*m*n
+
+
+    elements = zeros( (Nt,3), Int)
+    boundary = {}
+    Idg = []
+    Xg  = []
+    Idf = []
+    Xf  = []
+
+    full_send_dict = {}
+    ghost_recv_dict = {}
+    nt = -1
+    for i in range(m):
+        for j in range(n):
+
+            i1 = I(i,j+1)
+            i2 = I(i,j)
+            i3 = I(i+1,j+1)
+            i4 = I(i+1,j)
+
+            #Lower Element
+            nt = E(i,j)
+            if i == m-1:
+                #print 'nt =',nt
+                Idg.append(nt)
+                Idf.append(E(1,j))
+            if i == 0:
+                Idg.append(nt)
+                Idf.append(E(m-2,j))
+
+            if i == m-1:
+                boundary[nt, 2] = 'right'
+            if j == 0:
+                boundary[nt, 1] = 'bottom'
+            elements[nt,:] = [i4,i3,i2]
+
+            #Upper Element
+            nt = E(i,j)+1
+            if i == m-1:
+                Idg.append(nt)
+                Idf.append(E(1,j)+1)
+            if i == 0:
+                Idg.append(nt)
+                Idf.append(E(m-2,j)+1)
+
+            if i == 0:
+                boundary[nt, 2] = 'left'
+            if j == n-1:
+                boundary[nt, 1] = 'top'
+            elements[nt,:] = [i1,i2,i3]
+
+    Idf = array(Idf,Int)
+    Idg = array(Idg,Int)
+    Xf  = zeros(Idf.shape,Float)
+    Xg  = zeros(Idg.shape,Float)
+
+    #print Idf
+    #print Idg
+    full_send_dict[processor]  = [Idf, Xf]
+    ghost_recv_dict[processor] = [Idg, Xg]
+
+    return  points, elements, boundary, full_send_dict, ghost_recv_dict
+
+
+
+def rectangular_periodic(m, n, len1=1.0, len2=1.0, origin = (0.0, 0.0)):
+
+
+    """Setup a rectangular grid of triangles
+    with m+1 by n+1 grid points
+    and side lengths len1, len2. If side lengths are omitted
+    the mesh defaults to the unit square.
+
+    len1: x direction (left to right)
+    len2: y direction (bottom to top)
+
+    Return to lists: points and elements suitable for creating a Mesh or
+    FVMesh object, e.g. Mesh(points, elements)
     """
 
@@ -213 +358 @@
                 ghosts[nt] = E(m-2,j)
 
+            if j == n-1:
+                ghosts[nt] = E(i,1)
+
+            if j == 0:
+                ghosts[nt] = E(i,n-2)
+
             if i == m-1:
                 boundary[nt, 2] = 'right'
@@ -226 +377 @@
                 ghosts[nt] = E(m-2,j)+1
 
+            if j == n-1:
+                ghosts[nt] = E(i,1)+1
+
+            if j == 0:
+                ghosts[nt] = E(i,n-2)+1
+
             if i == 0:
                 boundary[nt, 2] = 'left'
@@ -232 +389 @@
             elements[nt,:] = [i1,i2,i3]
 
+    #bottom left
+    nt = E(0,0)
+    nf = E(m-2,n-2)
+    ghosts[nt]   = nf
+    ghosts[nt+1] = nf+1
+
+    #bottom right
+    nt = E(m-1,0)
+    nf = E(1,n-2)
+    ghosts[nt]   = nf
+    ghosts[nt+1] = nf+1
+
+    #top left
+    nt = E(0,n-1)
+    nf = E(m-2,1)
+    ghosts[nt]   = nf
+    ghosts[nt+1] = nf+1
+
+    #top right
+    nt = E(m-1,n-1)
+    nf = E(1,1)
+    ghosts[nt]   = nf
+    ghosts[nt+1] = nf+1
 
     return points, elements, boundary, ghosts
 
-def rectangular_periodic(m, n, len1=1.0, len2=1.0, origin = (0.0, 0.0)):
+def rectangular_periodic_lr(m, n, len1=1.0, len2=1.0, origin = (0.0, 0.0)):
 
 
@@ -247 +427 @@
 
     Return to lists: points and elements suitable for creating a Mesh or
-    FVMesh object, e.g. Mesh(points, elements)
+    Domain object, e.g. Mesh(points, elements)
     """
 
@@ -312 +492 @@
                 ghosts[nt] = E(m-2,j)
 
-            if j == n-1:
-                ghosts[nt] = E(i,1)
-
-            if j == 0:
-                ghosts[nt] = E(i,n-2)
-
             if i == m-1:
                 boundary[nt, 2] = 'right'
@@ -331 +505 @@
                 ghosts[nt] = E(m-2,j)+1
 
-            if j == n-1:
-                ghosts[nt] = E(i,1)+1
-
-            if j == 0:
-                ghosts[nt] = E(i,n-2)+1
-
             if i == 0:
                 boundary[nt, 2] = 'left'
@@ -343 +511 @@
             elements[nt,:] = [i1,i2,i3]
 
-    #bottom left
-    nt = E(0,0)
-    nf = E(m-2,n-2)
-    ghosts[nt]   = nf
-    ghosts[nt+1] = nf+1
-
-    #bottom right
-    nt = E(m-1,0)
-    nf = E(1,n-2)
-    ghosts[nt]   = nf
-    ghosts[nt+1] = nf+1
-
-    #top left
-    nt = E(0,n-1)
-    nf = E(m-2,1)
-    ghosts[nt]   = nf
-    ghosts[nt+1] = nf+1
-
-    #top right
-    nt = E(m-1,n-1)
-    nf = E(1,1)
-    ghosts[nt]   = nf
-    ghosts[nt+1] = nf+1
 
     return points, elements, boundary, ghosts
-
-def rectangular_periodic_lr(m, n, len1=1.0, len2=1.0, origin = (0.0, 0.0)):
-
-
-    """Setup a rectangular grid of triangles
-    with m+1 by n+1 grid points
-    and side lengths len1, len2. If side lengths are omitted
-    the mesh defaults to the unit square.
-
-    len1: x direction (left to right)
-    len2: y direction (bottom to top)
-
-    Return to lists: points and elements suitable for creating a Mesh or
-    Domain object, e.g. Mesh(points, elements)
-    """
-
-    from config import epsilon
-    from Numeric import zeros, Float, Int
-
-    delta1 = float(len1)/m
-    delta2 = float(len2)/n
-
-    #Calculate number of points
-    Np = (m+1)*(n+1)
-
-    class VIndex:
-
-        def __init__(self, n,m):
-            self.n = n
-            self.m = m
-
-        def __call__(self, i,j):
-            return j+i*(self.n+1)
-
-    class EIndex:
-
-        def __init__(self, n,m):
-            self.n = n
-            self.m = m
-
-        def __call__(self, i,j):
-            return 2*(j+i*self.n)
-
-
-    I = VIndex(n,m)
-    E = EIndex(n,m)
-
-    points = zeros( (Np,2), Float)
-
-    for i in range(m+1):
-        for j in range(n+1):
-
-            points[I(i,j),:] = [i*delta1 + origin[0], j*delta2 + origin[1]]
-
-    #Construct 2 triangles per rectangular element and assign tags to boundary
-    #Calculate number of triangles
-    Nt = 2*m*n
-
-
-    elements = zeros( (Nt,3), Int)
-    boundary = {}
-    ghosts = {}
-    nt = -1
-    for i in range(m):
-        for j in range(n):
-
-            i1 = I(i,j+1)
-            i2 = I(i,j)
-            i3 = I(i+1,j+1)
-            i4 = I(i+1,j)
-
-            #Lower Element
-            nt = E(i,j)
-            if i == m-1:
-                ghosts[nt] = E(1,j)
-            if i == 0:
-                ghosts[nt] = E(m-2,j)
-
-            if i == m-1:
-                boundary[nt, 2] = 'right'
-            if j == 0:
-                boundary[nt, 1] = 'bottom'
-            elements[nt,:] = [i4,i3,i2]
-
-            #Upper Element
-            nt = E(i,j)+1
-            if i == m-1:
-                ghosts[nt] = E(1,j)+1
-            if i == 0:
-                ghosts[nt] = E(m-2,j)+1
-
-            if i == 0:
-                boundary[nt, 2] = 'left'
-            if j == n-1:
-                boundary[nt, 1] = 'top'
-            elements[nt,:] = [i1,i2,i3]
-
-
-    return points, elements, boundary, ghosts

inundation/ga/storm_surge/parallel/run_parallel_advection.py

@@ -20 +20 @@
 M = 30
 
-points, vertices, boundary, ghosts = rectangular_periodic(N, M)
+points, vertices, boundary, full_send_dict, ghost_recv_dict = parallel_rectangular(N, M)
 
 #Create advection domain with direction (1,-1)
-domain = Parallel_Domain(points, vertices, boundary, ghosts, velocity=[1.0, 0.0])
+domain = Parallel_Domain(points, vertices, boundary,
+                         full_send_dict, ghost_recv_dict, velocity=[1.0, 0.0])
 
 # Initial condition is zero by default