Changeset 518

Timestamp:

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

Author:

ole

Message:

Looked at friction

Location:

inundation/ga/storm_surge/pyvolution

Files:

• mesh.py (modified) (1 diff)
• shallow_water.py (modified) (2 diffs)
• shallow_water_ext.c (modified) (1 diff)
• treenode.py (modified) (1 diff)

inundation/ga/storm_surge/pyvolution/mesh.py

@@ -433 +433 @@
 
 
+    def get_centroid_coordinates(self):
+        """Return all centroid coordinates.
+        Return all centroid coordinates for all triangles as an Nx2 array
+        (ordered as x0, y0 for each triangle)        
+        """
+        return self.centroid_coordinates

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):

inundation/ga/storm_surge/pyvolution/shallow_water_ext.c

@@ -168 +168 @@
       if (h >= eps) {
         S = -g * eta[k]*eta[k] * sqrt((uh[k]*uh[k] + vh[k]*vh[k]));
-        S /= pow(h, 7.0/3);      //Expensive
-        //S /= h*h*(1 + h/3.0 - h*h/9.0);  //FIXME: Use a Taylor expansion
-        //FIXME: This will save a lot of time
+        S /= pow(h, 7.0/3);      //Expensive (on Ole's home computer)
+        //S /= h*h*(1 + h/3.0 - h*h/9.0); //FIXME: Could use a Taylor expansion
+
 
         //Update momentum

inundation/ga/storm_surge/pyvolution/treenode.py

@@ -39 +39 @@
 
 
+    #Make treenode elements appear as sequences such thta on
+    #can iterate over them             
     def __iter__(self):
         self.index = -1   # incremented on first call to next()