Changeset 501

Timestamp:

Nov 8, 2004, 6:32:02 PM (20 years ago)

Author:

ole

Message:

First play with wind stress

File:

• inundation/ga/storm_surge/pyvolution/shallow_water.py (modified) (4 diffs)

inundation/ga/storm_surge/pyvolution/shallow_water.py

@@ -804 +804 @@
 
     w = domain.quantities['level'].centroid_values
+    z = domain.quantities['elevation'].centroid_values
+    h = w-z
+    
     uh = domain.quantities['xmomentum'].centroid_values
     vh = domain.quantities['ymomentum'].centroid_values
@@ -816 +819 @@
     
     for k in range(N):
-        if w[k] >= eps:
+        if h[k] >= eps:
             S = -g * eta[k]**2 * sqrt((uh[k]**2 + vh[k]**2))
-            S /= w[k]**(7.0/3)
+            S /= h[k]**(7.0/3)
 
             #Update momentum
@@ -834 +837 @@
         
     w = domain.quantities['level'].centroid_values
+    z = domain.quantities['elevation'].centroid_values
+    h = w-z
     uh = xmom.centroid_values
     vh = ymom.centroid_values
@@ -846 +851 @@
 
     from shallow_water_ext import manning_friction
-    manning_friction(g, eps, w, uh, vh, eta, xmom_update, ymom_update)
-        
+    manning_friction(g, eps, h, uh, vh, eta, xmom_update, ymom_update)
+        
+
+
+        
+def wind_stress(domain):
+    """Apply wind stress to water momentum
+    """
+
+    #FIXME: Under construction
+
+    #Hardwired velocity field - should be initialised as a
+    #callable object
+
+    u = 1.0
+    v = 0.0
+
+    Cw = 3.0e-3 #Wind stress coeffficient
+    rho_a = 1.2e-3 #Atmospheric density
+    rho = 1023 #Density of water
+    
+    from math import sqrt
+
+    xmom_update = domain.quantities['xmomentum'].semi_implicit_update
+    ymom_update = domain.quantities['ymomentum'].semi_implicit_update    
+    
+    N = domain.number_of_elements
+
+    c = Cw*rho_a/rho
+    S = c * sqrt(u**2 + v**2)    
+    for k in range(N):
+        #Update momentum
+        xmom_update[k] += S*u
+        ymom_update[k] += S*v
+            
+        
+
 
 ###########################