Changeset 1414

Timestamp:

May 17, 2005, 6:28:13 PM (20 years ago)

Author:

steve

Message:

 

Location:

inundation/ga/storm_surge

Files:

• Hobart/run_hobart_buildings.py (modified) (1 diff)
• parallel/parallel_advection.py (modified) (12 diffs)
• parallel/run_parallel_advection.py (modified) (2 diffs)
• pyvolution/realtime_visualisation_new.py (modified) (1 diff)
• zeus/pyvolution.zbi (modified) (previous)

inundation/ga/storm_surge/Hobart/run_hobart_buildings.py

@@ -42 +42 @@
 domain.default_order = 1
 domain.smooth = True
-domain.visualise = False
+domain.visualise = True
 
 #------------------------------

inundation/ga/storm_surge/parallel/parallel_advection.py

@@ -16 +16 @@
 from advection import *
 Advection_Domain = Domain
-from Numeric import zeros, Float, Int, ones
+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, velocity = None, processor=0):
+    def __init__(self, coordinates, vertices, boundary = None, ghosts = None,
+                 velocity = None):
 
         Advection_Domain.__init__(self, coordinates, vertices, boundary, velocity)
@@ -28 +29 @@
         N = self.number_of_elements
 
-        self.processor = processor
+        self.processor = pypar.rank()
+        self.numproc   = pypar.size()
 
         if  ghosts == None:
@@ -52 +54 @@
 
     def update_timestep(self, yieldstep, finaltime):
+
+        # Calculate local timestep
         Advection_Domain.update_timestep(self, yieldstep, finaltime)
 
-        timestep = 0.0
-        pypar.raw_reduce(self.timestep,timestep,pypar.MIN,0)
-        print 'Processor %d, Timestep %.4f'%(self.processor,timestep)
-
+        # For some reason it looks like pypar only reduces numeric arrays
+        # hence we need to create some dummy arrays for communication
+        ltimestep = ones( 1, Float )
+        ltimestep[0] = self.timestep
+
+        gtimestep = zeros( 1, Float) # Buffer for results
+
+        pypar.raw_reduce(ltimestep, gtimestep, pypar.MIN, 0)
+        pypar.broadcast(gtimestep,0)
+        pypar.Barrier()
+
+        self.timestep = gtimestep[0]
 
 
 
     def update_ghosts(self):
+
+        # We must send the information from the full cells and
+        # receive the information for the ghost cells
+        # We have a dictionary of lists with ghosts expecting updates from
+        # the separate processors
+
+        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)
+#
+#        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
@@ -101 +136 @@
 
 
-def rectangular_periodic(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)):
 
 
@@ -107 +143 @@
     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.
+    the mesh defaults to the unit square, divided between all the
+    processors
 
     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)
     """
 
@@ -178 +213 @@
                 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'
@@ -197 +226 @@
                 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'
@@ -209 +232 @@
             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)):
+def rectangular_periodic(m, n, len1=1.0, len2=1.0, origin = (0.0, 0.0)):
 
 
@@ -312 +312 @@
                 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'
@@ -325 +331 @@
                 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'
@@ -331 +343 @@
             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

@@ -15 +15 @@
 processor_name = pypar.Get_processor_name()
 
+
+
 N = 30
 M = 30
@@ -21 +23 @@
 
 #Create advection domain with direction (1,-1)
-domain = Parallel_Domain(points, vertices, boundary, ghosts, velocity=[1.0, 0.0],processor = myid)
+domain = Parallel_Domain(points, vertices, boundary, ghosts, velocity=[1.0, 0.0])
 
 # Initial condition is zero by default

inundation/ga/storm_surge/pyvolution/realtime_visualisation_new.py

@@ -21 +21 @@
         self.z_models = {}
         keys = self.domain.quantities.keys()
-        print keys
+        #print keys
         for key in keys:
             self.z_models[key] = faces(frame= self.frame)
 
-        print self.z_models
+        #print self.z_models
 
         self.max_x = max(max(self.vertices[:,0]),max(self.vertices[:,2]),max(self.vertices[:,4]))