Changeset 1393

Timestamp:

May 15, 2005, 10:25:56 PM (19 years ago)

Author:

steve

Message:

Added ghosts to rectangle

Location:

inundation/ga/storm_surge

Files:

• Merimbula/run_merimbula_lake.py (modified) (1 diff)
• parallel/parallel_advection.py (modified) (3 diffs)
• parallel/run_parallel_advection.py (modified) (2 diffs)
• pyvolution/data_manager.py (modified) (2 diffs)
• pyvolution/netherlands.py (modified) (3 diffs)
• pyvolution/normalDate.py (added)
• pyvolution/quantity.py (modified) (1 diff)
• pyvolution/test_quantity.py (modified) (3 diffs)
• zeus/anuga-workspace.zwi (modified) (1 diff)
• zeus/pyvolution.zbi (modified) (previous)
• zeus/pyvolution.zpi (modified) (1 diff)

inundation/ga/storm_surge/Merimbula/run_merimbula_lake.py

@@ -79 +79 @@
 domain.store = True     #Store for visualisation purposes
 domain.format = 'sww'   #Native netcdf visualisation format
-domain.filename = 'Merimbula_2003_4days_dry'
+from normalDate import ND
+domain.filename = 'Merimbula_2003_4days_dry_%s'%ND()
+
+print domain.filename
 
 #----------------------

inundation/ga/storm_surge/parallel/parallel_advection.py

@@ -66 +66 @@
 def rectangular_with_ghosts(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
@@ -74 +75 @@
     len2: y direction (bottom to top)
 
-    Also returns a list of
-
     Return to lists: points and elements suitable for creating a Mesh or
     FVMesh object, e.g. Mesh(points, elements)
@@ -81 +80 @@
 
     from config import epsilon
-
-    #E = m*n*2        #Number of triangular elements
-    #P = (m+1)*(n+1)  #Number of initial vertices
+    from Numeric import zeros, Float, Int
 
     delta1 = float(len1)/m
     delta2 = float(len2)/n
 
-    #Dictionary of vertex objects
-    vertices = {}
-    points = []
+    #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):
-            vertices[i,j] = len(points)
-            points.append([i*delta1 + origin[0], j*delta2 + origin[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
 
 
-
-    #Construct 2 triangles per rectangular element and assign tags to boundary
-    elements = []
+    elements = zeros( (Nt,3), Int)
     boundary = {}
+    ghosts = {}
+    nt = -1
     for i in range(m):
         for j in range(n):
-            v1 = vertices[i,j+1]
-            v2 = vertices[i,j]
-            v3 = vertices[i+1,j+1]
-            v4 = vertices[i+1,j]
 
-            #Update boundary dictionary and create elements
+            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:
-                boundary[(len(elements), 2)] = 'right'
+                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[(len(elements), 1)] = 'bottom'
-            elements.append([v4,v3,v2]) #Lower element
+                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[(len(elements), 2)] = 'left'
+                boundary[nt, 2] = 'left'
             if j == n-1:
-                boundary[(len(elements), 1)] = 'top'
-            elements.append([v1,v2,v3]) #Upper element
+                boundary[nt, 1] = 'top'
+            elements[nt,:] = [i1,i2,i3]
 
-    ghosts = {}
-    i=0
-    for j in range(n):
-        v1 = vertices[i,j+1]
-        v2 = vertices[i,j]
-        v3 = vertices[i+1,j+1]
-        v4 = vertices[i+1,j]
-        ghosts.append(elements.index([
+    return points, elements, boundary, ghosts
 
-    return points, elements, boundary

inundation/ga/storm_surge/parallel/run_parallel_advection.py

@@ -11 +11 @@
 from mesh_factory import rectangular
 
-points, vertices, boundary = rectangular(30, 30)
+N = 4
+M = 4
+
+points, vertices, boundary, ghosts = rectangular_with_ghosts(N, M)
+
+
+print ghosts
 
 #Create advection domain with direction (1,-1)
@@ -47 +53 @@
 
 #Check that the boundary value gets propagated to all elements
-for t in domain.evolve(yieldstep = 0.01, finaltime = 1.5):
+for t in domain.evolve(yieldstep = 0.01, finaltime = 0.01):
     domain.write_time()

inundation/ga/storm_surge/pyvolution/data_manager.py

@@ -342 +342 @@
         from Scientific.IO.NetCDF import NetCDFFile
         import types
-        from time import sleep
+        from time import sleep
         from os import stat
 
 
         #Get NetCDF
-        retries = 0
-        file_open = False
-        while not file_open and retries < 10:
-            try:
+        retries = 0
+        file_open = False
+        while not file_open and retries < 10:
+            try:
                 fid = NetCDFFile(self.filename, 'a') #Open existing file
-            except IOError:
-                #This could happen if someone was reading the file.
-                #In that case, wait a while and try again
-                msg = 'Warning (store_timestep): File %s could not be opened'\
-                %self.filename
-                msg += ' - trying again'
-                print msg
-                retries += 1
-                sleep(1)
-            else:
-               file_open = True
-
-        if not file_open:
-            msg = 'File %s could not be opened for append' %self.filename
-            raise msg
+            except IOError:
+                #This could happen if someone was reading the file.
+                #In that case, wait a while and try again
+                msg = 'Warning (store_timestep): File %s could not be opened'\
+                      %self.filename
+                msg += ' - trying again'
+                print msg
+                retries += 1
+                sleep(1)
+            else:
+                file_open = True
+
+        if not file_open:
+            msg = 'File %s could not be opened for append' %self.filename
+            raise msg
 
 
@@ -1995 +1995 @@
         q = (1-ratio)*Q[index]+ ratio*Q[index+1]
     else:
-         q = Q[index]
+        q = Q[index]
     #Return vector of interpolated values
     return q

inundation/ga/storm_surge/pyvolution/netherlands.py

@@ -65 +65 @@
 
             #Wall
-            #if x[i] > 0.995:
-            #    z[i] = -x[i]/20+0.3
+            if x[i] > 0.995:
+                z[i] = -x[i]/20+0.3
 
         return z/2
@@ -103 +103 @@
 
 #Set bed-slope and friction
-inflow_stage = 0.2
-manning = 0.00
+inflow_stage = 0.1
+manning = 0.03
 Z = Weir(inflow_stage)
 
@@ -156 +156 @@
 t0 = time.time()
 
-for t in domain.evolve(yieldstep = 1.0, finaltime = 100.0):
+from visual import *
+
+for t in domain.evolve(yieldstep = 0.1, finaltime = 15.0):
     domain.write_time()
     #domain.visualiser.scale_z = 1.0
-    domain.visualiser.update_quantity_color('xmomentum',scale_z = 4.0)
+    #domain.visualiser.update_quantity_color('xmomentum',scale_z = 4.0)
+    #integral_label.text='Integral=%10.5e'%domain.quantities['stage'].get_integral()
 
 

inundation/ga/storm_surge/pyvolution/quantity.py

@@ -536 +536 @@
 
             # Create a V like [[0 1 2], [3 4 5]....[3*m-2 3*m-1 3*m]]
-            # These vert_id's will relate to the verts created bellow
+            # These vert_id's will relate to the verts created below
             m = len(self.domain)  #Number of volumes
             M = 3*m        #Total number of unique vertices

inundation/ga/storm_surge/pyvolution/test_quantity.py

@@ -201 +201 @@
         quantity = Quantity(self.mesh4)
 
-        #Try constants first
-        const = 5
+        #Try constants first
+        const = 5
         quantity.set_values(const, location = 'vertices')
-        #print 'Q', quantity.get_integral()
-
-        assert allclose(quantity.get_integral(),
-                        self.mesh4.get_area() * const)
-
-        #Try with a linear function
+        #print 'Q', quantity.get_integral()
+
+        assert allclose(quantity.get_integral(), self.mesh4.get_area() * const)
+
+        #Try with a linear function
         def f(x, y):
             return x+y
@@ -216 +215 @@
 
 
-        ref_integral = (4.0/3 + 8.0/3 + 10.0/3 + 10.0/3) * 2
-
-        assert allclose (quantity.get_integral(), ref_integral)
+        ref_integral = (4.0/3 + 8.0/3 + 10.0/3 + 10.0/3) * 2
+
+        assert allclose (quantity.get_integral(), ref_integral)
 
 
@@ -929 +928 @@
     runner = unittest.TextTestRunner()
     runner.run(suite)
-

inundation/ga/storm_surge/zeus/anuga-workspace.zwi

@@ -2 +2 @@
 <workspace name="anuga-workspace">
     <mode>Debug</mode>
-    <active>steve</active>
+    <active>pyvolution</active>
     <project name="Merimbula">Merimbula.zpi</project>
     <project name="parallel">parallel.zpi</project>

inundation/ga/storm_surge/zeus/pyvolution.zpi

@@ -99 +99 @@
     <file>..\pyvolution\netherlands.py</file>
     <file>..\pyvolution\netherlands-chris.py</file>
+    <file>..\pyvolution\normalDate.py</file>
     <file>..\pyvolution\pmesh2domain.py</file>
     <file>..\pyvolution\pressure_force.py</file>