Changeset 610

Timestamp:

Nov 22, 2004, 4:43:09 PM (20 years ago)

Author:

ole

Message:

 

Location:

inundation/ga/storm_surge/pyvolution

Files:

• shallow_water.py (modified) (4 diffs)
• test_shallow_water.py (modified) (3 diffs)
• wind_example_scalar.py (moved) (moved from inundation/ga/storm_surge/pyvolution/wind_example.py)
• wind_example_variable.py (moved) (moved from inundation/ga/storm_surge/pyvolution/wind_example_rot.py)

inundation/ga/storm_surge/pyvolution/shallow_water.py

@@ -924 +924 @@
             return phi
 
-        
+        where x and y are vectors.
 
         and then pass the functions in
@@ -950 +950 @@
 
         from math import pi, cos, sin, sqrt
+        from Numeric import Float, ones 
         
         xmom_update = domain.quantities['xmomentum'].explicit_update
@@ -955 +956 @@
         
         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        
-            
-                
+        t = domain.time        
+        
+        if callable(self.s):
+            xc = domain.get_centroid_coordinates()            
+            x = xc[:,0]
+            y = xc[:,1]
+            
+            s_vec = self.s(x,y,t)
         else:
-            #Constants
-            
+            #Assume s is a scalar
+
+            try:
+                s_vec = self.s * ones(N, Float)
+            except:
+                msg = 'Speed must be either callable or a scalar: %s' %self.s
+                raise msg
+
+
+        if callable(self.phi):
+            xc = domain.get_centroid_coordinates()            
+            x = xc[:,0]
+            y = xc[:,1]
+            
+            phi_vec = self.phi(x,y,t)
+        else:
+            #Assume phi is a scalar
+
+            try:
+                phi_vec = self.phi * ones(N, Float)            
+            except:
+                msg = 'Angle must be either callable or a scalar: %s' %self.phi
+                raise msg
+
+        #This is the bit that should be written in C    
+        for k in range(N):
+            s = s_vec[k]
+            phi = phi_vec[k]
+
             #Convert to radians
-            phi = self.phi*pi/180
-        
+            phi = phi*pi/180
+
             #Compute velocity vector (u, v)
-            u = self.s*cos(phi)
-            v = self.s*sin(phi)
+            u = s*cos(phi)
+            v = s*sin(phi)
 
             #Compute wind stress
             S = self.const * sqrt(u**2 + v**2)
-            xmom_update += S*u
-            ymom_update += S*v        
-        
-        
+            xmom_update[k] += S*u
+            ymom_update[k] += S*v        
+            
 
 class Wind_stress_from_file:
@@ -1015 +1014 @@
     
     FIXME: This class may be incorporated in the generic Wind_stress class
-    """
-
+    SUperseded
+    """
+
+    
     def __init__(self, filename):
         """Initialise windfield from a file with the following format

inundation/ga/storm_surge/pyvolution/test_shallow_water.py

@@ -17 +17 @@
     
     from math import sqrt, exp, cos, pi
+
+    N = len(x)
+    s = 0*x  #New array
     
-    r = sqrt(x**2 + y**2)
-
-    factor = exp( -(r-0.15)**2 )
-
-    return 4000 * factor * (cos(t*2*pi/150) + 2)
+    for k in range(N):
+        
+        r = sqrt(x[k]**2 + y[k]**2)
+
+        factor = exp( -(r-0.15)**2 )
+
+        s[k] = 4000 * factor * (cos(t*2*pi/150) + 2)
+
+    return s
 
 
@@ -30 +37 @@
     from math import sqrt, atan, pi
 
-    r = sqrt(x**2 + y**2)    
+
+    N = len(x)
+    a = 0*x  #New array
+
+    for k in range(N):    
+        r = sqrt(x[k]**2 + y[k]**2)    
     
-    angle = atan(y/x)
-
-    if x < 0:
-        angle+=pi #atan in ]-pi/2; pi/2[
-
-    #Take normal direction
-    angle -= pi/2
+        angle = atan(y[k]/x[k])
+
+        if x[k] < 0:
+            angle+=pi #atan in ]-pi/2; pi/2[
+
+        #Take normal direction
+        angle -= pi/2
     
-    #Ensure positive radians    
-    if angle < 0:
-        angle += 2*pi
-
-    return angle/pi*180
-
+        #Ensure positive radians    
+        if angle < 0:
+            angle += 2*pi
+
+        a[k] = angle/pi*180
+        
+    return a
 
         
@@ -791 +804 @@
         t = domain.time
 
+        x = xc[:,0]
+        y = xc[:,1]
+        s_vec = speed(x,y,t)
+        phi_vec = angle(x,y,t)
+        
+
         for k in range(N):
-            x, y = xc[k,:]    
-
-            s = speed(x,y,t)
-            phi = angle(x,y,t)
-
-
             #Convert to radians
-            phi = phi*pi/180
+            phi = phi_vec[k]*pi/180
+            s = s_vec[k]
         
             #Compute velocity vector (u, v)