Changeset 518 for inundation/ga/storm_surge/pyvolution/shallow_water.py

Timestamp:

Nov 10, 2004, 5:22:27 PM (20 years ago)

Author:

ole

Message:

Looked at friction

File:

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

inundation/ga/storm_surge/pyvolution/shallow_water.py

@@ -745 +745 @@
 #
 def gravity(domain):
-    """Implement forcing function for bed slope working with
-    consecutive data structures of class Volume
+    """Apply gravitational pull in the presence of bed slope
     """
 
@@ -855 +854 @@
         
 
-
+class Wind_stress:
+    """Apply wind stress to water momentum
+    """
+
+    def __init__(self, s, phi):
+        """Initialise windfield from wind speed s [m/s]
+        and wind direction phi [degrees]
+        
+        Inputs v and phi can be either scalars or Python functions, e.g.
+
+        W = Wind_stress(10, 178)
+
+        #FIXME - 'normal' degrees are assumed for now, i.e. the
+        vector (1,0) has zero degrees. 
+        We may need to convert from 'compass' degrees later on and also
+        map from True north to grid north.
+
+        Arguments can also be Python functions of x,y,t as in
+
+        def speed(x,y,t):
+            ...
+            return s
+        
+        def angle(x,y,t):
+            ...
+            return phi
+
+        
+
+        and then pass the functions in
+
+        W = Wind_stress(speed, angle)
+
+        The instantiated object W can be appended to the list of
+        forcing_terms as in
+
+        domain.forcing_terms.append(W)
+        
+        """
+
+        self.s = s
+        self.phi = phi
+
+        #FIXME: Maybe move to config or read from domain
+        Cw = 3.0e-3 #Wind stress coefficient
+        rho_a = 1.2e-3 #Atmospheric density
+        rho = 1023 #Density of water
+        
+        self.const = Cw*rho_a/rho
+        
+
+    def __call__(self, domain):
+        """Evaluate windfield based on values found in domain
+        """
+
+        from math import pi, cos, sin, sqrt
+        
+        xmom_update = domain.quantities['xmomentum'].explicit_update
+        ymom_update = domain.quantities['ymomentum'].explicit_update    
+        
+        N = domain.number_of_elements
+        
+        if callable(self.s) or callable(self.phi):
+            xc = domain.get_centroid_coordinates()
+            t = domain.time
+
+            maxphi = -100
+            minphi = 100
+            for k in range(N):
+                x, y = xc[k,:]    
+
+                if callable(self.s):
+                    s = self.s(x,y,t)
+                else:
+                    s = self.s
+
+                if callable(self.phi):
+                    phi = self.phi(x,y,t)
+                else:
+                    phi = self.phi                  
+
+                #Convert to radians
+                phi = phi*pi/180
+        
+                #Compute velocity vector (u, v)
+                u = s*cos(phi)
+                v = s*sin(phi)
+
+                #Compute wind stress
+                S = self.const * sqrt(u**2 + v**2)
+                xmom_update[k] += S*u
+                ymom_update[k] += S*v        
+            
+                
+        else:
+            #Constants
+            
+            #Convert to radians
+            phi = self.phi*pi/180
+        
+            #Compute velocity vector (u, v)
+            u = self.s*cos(phi)
+            v = self.s*sin(phi)
+
+            #Compute wind stress
+            S = self.const * sqrt(u**2 + v**2)
+            xmom_update += S*u
+            ymom_update += S*v        
+        
         
 def wind_stress(domain):