Changeset 1402

Timestamp:

May 16, 2005, 5:25:56 PM (20 years ago)

Author:

steve

Message:

Periodic boundary working via ghosts

Location:

inundation/ga/storm_surge

Files:

• parallel/parallel_advection.py (modified) (6 diffs)
• parallel/run_parallel_advection.py (modified) (2 diffs)
• pyvolution/domain.py (modified) (3 diffs)
• zeus/pyvolution.zbi (modified) (previous)

inundation/ga/storm_surge/parallel/parallel_advection.py

@@ -16 +16 @@
 from advection import *
 Advection_Domain = Domain
-from Numeric import zeros, Float, Int
+from Numeric import zeros, Float, Int, ones
 
 class Parallel_Domain(Advection_Domain):
 
-    def __init__(self, coordinates, vertices, boundary = None, velocity = None,
-                 processor = 0, global_ids = None):
+    def __init__(self, coordinates, vertices, boundary = None, ghosts = None, velocity = None):
 
         Advection_Domain.__init__(self, coordinates, vertices, boundary, velocity)
 
-        self.processor = processor
-
         N = self.number_of_elements
 
-        if  global_ids == None:
-            self.global_ids = global_ids
+        if  ghosts == None:
+            self.ghosts = None
+            self.full = ones( N, Float)
         else:
-            self.global_ids = zeros(N, Int)
-
+            self.ghosts = ghosts
+            self.full = ones( N, Int)
+            keys = self.ghosts.keys()
+            for key in keys:
+                self.full[key] = -1
 
     def check_integrity(self):
@@ -42 +43 @@
 
 
-
+    def update_ghosts(self):
+        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]]
 
     def evolve(self, yieldstep = None, finaltime = None):
@@ -64 +69 @@
 
 
-def rectangular_with_ghosts(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)):
 
 
@@ -141 +146 @@
                 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'
@@ -154 +165 @@
                 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'
@@ -160 +177 @@
             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
+    FVMesh 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

@@ -9 +9 @@
 from Numeric import array
 
-from mesh_factory import rectangular
+N = 40
+M = 40
 
-N = 4
-M = 4
-
-points, vertices, boundary, ghosts = rectangular_with_ghosts(N, M)
-
-
-print ghosts
+points, vertices, boundary, ghosts = rectangular_periodic_lr(N, M)
 
 #Create advection domain with direction (1,-1)
-domain = Parallel_Domain(points, vertices, boundary, velocity=[1.0, 0.0])
+domain = Parallel_Domain(points, vertices, boundary, ghosts, velocity=[1.0, 0.0])
 
 # Initial condition is zero by default
@@ -53 +48 @@
 
 #Check that the boundary value gets propagated to all elements
-for t in domain.evolve(yieldstep = 0.01, finaltime = 0.01):
+for t in domain.evolve(yieldstep = 0.01, finaltime = 1.0):
     domain.write_time()

inundation/ga/storm_surge/pyvolution/domain.py

@@ -431 +431 @@
         self.distribute_to_vertices_and_edges()
 
+        #Initial update boundary values
+        self.update_boundary()
 
         #Or maybe restore from latest checkpoint
@@ -436 +438 @@
             self.goto_latest_checkpoint()
 
-
         yield(self.time)  #Yield initial values
 
         while True:
+
+            #Compute fluxes across each element edge
+            self.compute_fluxes()
+
+            #Update timestep to fit yieldstep and finaltime
+            self.update_timestep(yieldstep, finaltime)
+
+            #Update conserved quantities
+            self.update_conserved_quantities()
+
+            #update ghosts
+            self.update_ghosts()
+
+            #Update vertex and edge values
+            self.distribute_to_vertices_and_edges()
+
             #Update boundary values
             self.update_boundary()
-
-            #Compute fluxes across each element edge
-            self.compute_fluxes()
-
-            #Update timestep to fit yieldstep and finaltime
-            self.update_timestep(yieldstep, finaltime)
-
-            #Update conserved quantities
-            self.update_conserved_quantities()
-
-            #Update vertex and edge values
-            self.distribute_to_vertices_and_edges()
 
             #Update time
@@ -604 +609 @@
             Q.semi_implicit_update[:] = 0.0
 
+    def update_ghosts(self):
+        pass