Changeset 1587

Timestamp:

Jul 7, 2005, 3:47:06 PM (19 years ago)

Author:

steve

Message:

 

Location:

inundation/ga/storm_surge

Files:

• Hobart/Landslide_1D.for (modified) (6 diffs)
• Hobart/Landslide_Tadmor.for (modified) (1 diff)
• pyvolution/general_mesh.py (modified) (1 diff)
• pyvolution/mesh_factory.py (modified) (1 diff)

inundation/ga/storm_surge/Hobart/Landslide_1D.for

@@ -39 +39 @@
 c      
       nx = 101
-      dtmax  = 1.0d0
+      dtmax  = 0.5d0
       Tout   = 1.0d0
-      tf = 200.0d0
+      tf = 1500.0d0
       length = 1000.0d0
 
@@ -49 +49 @@
       
       limith = 1
-      dx = 10.d0
+      
+      dx = 10000.d0/dble(nx-1)
       K = 0.0d0            
 
       if (nx.gt.nxmax) write(*,*)'nx too large'
       toll = 1.0d-10
-      cfl = 0.01d0
+      cfl = 0.2d0
       theta = 1.2d0
       version = 10
@@ -310 +311 @@
          w1 = h1 + z(i)
          w2 = h2 + z(i+1)
-         
+ 
+         d1h = h1 - h0
+         d2h = h2 - h1
          d1w = w1 - w0
          d2w = w2 - w1
@@ -316 +319 @@
                 
 c ------- W = H+Z           
-         h_x = xmin( d1w, d2w) - dz
+c         h_x = xmin( d1w, d2w) - dz
+         h_x = xmic(theta, d1h, d2h)
          
 c ------- Return the left and right values of h in an interval
@@ -323 +327 @@
          
 c ------- Test for h_l or h_r negative
-         hmin = dmin1(h0,h1,h2)
-         if ( h_l(i) .le. hmin ) then
-            h_l(i) = hmin
-            h_r(i) = h(i) + (h(i) - h_l(i))
-         elseif ( h_r(i) .le. hmin ) then
-            h_r(i) = hmin
-            h_l(i) = h(i) + (h(i) - h_r(i))
-         endif
+c$$$         hmin = dmin1(h0,h1,h2)
+c$$$         if ( h_l(i) .le. hmin ) then
+c$$$            h_l(i) = hmin
+c$$$            h_r(i) = h(i) + (h(i) - h_l(i))
+c$$$         elseif ( h_r(i) .le. hmin ) then
+c$$$            h_r(i) = hmin
+c$$$            h_l(i) = h(i) + (h(i) - h_r(i))
+c$$$         endif
       enddo
 c
 c ---- Enforce a minimum height
 c
-      do i=md,nx+md+2
-         hmax = dmax1(h_l(i),h_r(i))
-         if (hmax .lt. 1.0d-3 ) then
-           h_r(i)  = 0.0d0
-           uh_r(i) = 0.0d0
-           h_l(i)  = 0.0d0
-           uh_l(i) = 0.0d0
-           h(i)    = 0.0d0
-           uh(i)   = 0.0d0
-         endif
-      enddo
+c$$$      do i=md,nx+md+2
+c$$$         hmax = dmax1(h_l(i),h_r(i))
+c$$$         if (hmax .lt. 1.0d-3 ) then
+c$$$c           h_r(i)  = 0.0d0
+c$$$           uh_r(i) = 0.0d0
+c$$$c           h_l(i)  = 0.0d0
+c$$$           uh_l(i) = 0.0d0
+c$$$c           h(i)    = 0.0d0
+c$$$           uh(i)   = 0.0d0
+c$$$         endif
+c$$$      enddo
 c
 c ---- Try for a flat subinterval 
@@ -410 +414 @@
           d2 = u(i+1) - u(i)
           u_x = xmic( theta, d1, d2 ) 
+c          u_x = xmin( d1, d2 ) 
           u_l(i) = u(i) - 0.5d0*u_x
           u_r(i) = u(i) + 0.5d0*u_x          

inundation/ga/storm_surge/Hobart/Landslide_Tadmor.for

@@ -37 +37 @@
       theta = 1.0d0
       g = 9.81d0
-      hl=10.0d0
+      hl= 5.0d0
       hr = 0.d0
       dx = 10000.d0/dble(nx-1)

inundation/ga/storm_surge/pyvolution/general_mesh.py

@@ -85 +85 @@
         self.number_of_elements = N = self.triangles.shape[0]
 
-        self.xy_extent = [ min(self.coordinates[:,0]), min(self.coordinates[:,1]) ,
-                           max(self.coordinates[:,0]), max(self.coordinates[:,1]) ]
+        xy_extent = [ min(self.coordinates[:,0]), min(self.coordinates[:,1]) ,
+                      max(self.coordinates[:,0]), max(self.coordinates[:,1]) ]
+
+        self.xy_extent = array(xy_extent, Float)
 
 

inundation/ga/storm_surge/pyvolution/mesh_factory.py

@@ -517 +517 @@
 
     return points, elements, boundary
+
+
+def contracting_channel_cross(m, n, W_upstream = 1., W_downstream = 0.75,
+                              L_1 = 5.0, L_2 = 2.0, L_3 = 10, origin = (0.0, 0.0)):
+    """Setup a contracting channel grid of triangles
+    with m segments in the x-direction and n segments in the y-direction
+
+    """
+
+    from Numeric import array
+    import math
+
+    from config import epsilon
+
+
+    lenx = L_1 + L_2 + L_3
+    leny = W_upstream
+    deltax = lenx/float(m)
+    deltay = leny/float(n)
+
+    x1 = 0
+    y1 = 0
+    x2 = L_1
+    y2 = 0
+    x3 = L_1 + L_2
+    y3 = (W_upstream - W_downstream)/2
+    x4 = L_1 + L_2 + L_3
+    y4 = y3
+    x5 = x4
+    y5 = y4 + W_downstream
+    x6 = L_1 + L_2
+    y6 = y5
+    x7 = L_1
+    y7 = W_upstream
+    x8 = 0
+    y8 = W_upstream
+    a1 = 0
+    a2 = (W_upstream - W_downstream)/(2*L_2)
+    a3 = 1
+    a4 = (W_downstream - W_upstream)/(L_2*W_upstream)
+
+    # Dictionary of vertex objects
+    vertices = {}
+    points = []
+
+    for i in range(m+1):
+        x = deltax*i
+        for j in range(n+1):
+            y = deltay*j
+            if x > L_1 and x <= (L_1 + L_2):
+                y = a1 + a2*(x - L_1) + a3*y + a4*(x - L_1)*y
+            elif x > L_1 + L_2:
+                y = (W_upstream - W_downstream)/2 + deltay*j*W_downstream/W_upstream
+
+            vertices[i,j] = len(points)
+            points.append([x + origin[0], y + origin[1]])
+
+    # Construct 4 triangles per element
+    elements = []
+    boundary = {}
+    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]
+            x = (points[v1][0]+points[v2][0]+points[v3][0]+points[v4][0])*0.25
+            y = (points[v1][1]+points[v2][1]+points[v3][1]+points[v4][1])*0.25
+            v5 = len(points)
+            points.append([x + origin[0], y + origin[1]])
+
+            #Create left triangle
+            if i == 0:
+                boundary[(len(elements), 1)] = 'left'
+            elements.append([v2,v5,v1])
+
+            #Create bottom triangle
+            if j == 0:
+                boundary[(len(elements), 1)] = 'bottom'
+            elements.append([v4,v5,v2])
+
+            #Create right triangle
+            if i == m-1:
+                boundary[(len(elements), 1)] = 'right'
+            elements.append([v3,v5,v4])
+
+            #Create top triangle
+            if j == n-1:
+                boundary[(len(elements), 1)] = 'top'
+            elements.append([v1,v5,v3])
+
+
+    return points, elements, boundary