Changeset 614 for inundation/ga/storm_surge/pyvolution/wind_example_variable.py

Timestamp:

Nov 23, 2004, 5:55:46 PM (20 years ago)

Author:

ole

Message:

Played with wind stress + c-extension

File:

• inundation/ga/storm_surge/pyvolution/wind_example_variable.py (modified) (5 diffs)

inundation/ga/storm_surge/pyvolution/wind_example_variable.py

@@ -9 +9 @@
 #
 from mesh_factory import rectangular
-from shallow_water import Domain, Reflective_boundary, Wind_stress
+from shallow_water import Domain, Dirichlet_boundary, Wind_stress
 
 #Create basic mesh
-N = 80
-length = 100
+N = 120
+length = 200
 points, vertices, boundary = rectangular(N, N, length, length)
 
@@ -25 +25 @@
 
 #Set initial conditions
+h = 1.0
 domain.set_quantity('elevation', 0.0)
-domain.set_quantity('level', 2.0)
-domain.set_quantity('friction', 0.07)
+domain.set_quantity('level', h)
+domain.set_quantity('friction', 0.01)
 
 
@@ -36 +37 @@
     """
 
-    from math import sqrt, exp, cos, pi
-
-    N = len(x)
-    s = 0*x  #New array
+    from math import pi
+    from Numeric import sqrt, exp, cos
 
     c = (length/2, length/2)
-    
-    for k in range(N):
-        r = sqrt((x[k] - c[0])**2 + (y[k] - c[1])**2)
+    r = sqrt((x - c[0])**2 + (y - c[1])**2)/length
+    factor = exp( -(r-0.15)**2 )
 
-        r /= max(length,length)
-
-        factor = exp( -(r-0.15)**2 )
-        
-        s[k] = 4000 * factor * (cos(t*2*pi/150) + 2)
-
-    return s
+    #return 9000 * factor
+    return 4000 * factor * (cos(t*2*pi/150) + 2)
 
 
@@ -58 +51 @@
     """Rotating field
     """
-    from math import sqrt, atan, pi
-
-    N = len(x)
-    a = 0*x  #New array
+    
+    from math import pi
+    from Numeric import sqrt, exp, cos, arctan2, choose, less
 
     c = (length/2, length/2)
-    for k in range(N):
-        r = sqrt((x[k] - c[0])**2 + (y[k] - c[1])**2)
+    xx = (x - c[0])/length
+    yy = (y - c[1])/length
+    angle = arctan2(yy,xx)
 
-        xx = (x[k] - c[0])/length
-        yy = (y[k] - c[1])/length
+    #Take normal direction (but reverse direction every 50 seconds)
+    #if sin(t*2*pi/100) < 0:
+    #    sgn = -1
+    #else:
+    #    sgn = 1
+    #angle += sgn*pi/2
+    angle -= pi/2    
 
-        angle = atan(yy/xx)
-
-        if xx < 0:
-            angle+=pi #atan in ]-pi/2; pi/2[
-
-        #Take normal direction
-        angle -= pi/2
-    
-        #Ensure positive radians    
-        if angle < 0:
-            angle += 2*pi
-
-        a[k] = angle/pi*180
-        
-    return a
+    #Convert to degrees and return
+    return angle/pi*180
 
     
@@ -109 +94 @@
 ######################
 # Boundary conditions
-Br = Reflective_boundary(domain)
-domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
+#Br = Reflective_boundary(domain)
+Bd = Dirichlet_boundary([h, 0, 0])
+domain.set_boundary({'left': Bd, 'right': Bd, 'top': Bd, 'bottom': Bd})
 
 ######################