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

Timestamp:

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

Author:

ole

Message:

Work towards wind_field - it now accepts vect values fuctions or scalars

File:

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

inundation/ga/storm_surge/pyvolution/shallow_water.py

@@ -896 +896 @@
         
         
+def check_forcefield(f):
+    """Check that f is either
+    1: a callable object of x,y,t, where x and y are vectors
+       and that it returns an array or a list of same length
+       as x and y
+    2: a scalar   
+    """
+
+    from Numeric import ones, Float, array
+
+    
+    if callable(f):
+        N = 3
+        x = ones(3, Float)
+        y = ones(3, Float)            
+        try:
+            q = f(x, y, 1.0)
+        except Exception, e:
+            msg = 'Function %s could not be executed:\n%s' %(f, e)
+            raise msg
+                    
+        try:
+            q = array(q).astype(Float)
+        except:
+            msg = 'Return value from vector function %s could ' %f
+            msg += 'not be converted into a Numeric array of floats.\n'
+            msg += 'Specified function should return either list or array.' 
+            raise msg
+
+        msg = 'Return vector from function %s' %f
+        msg += 'must have same lenght as input vectors'
+        assert len(q) == N, msg
+        
+    else:        
+        try:
+            f = float(f)
+        except:
+            msg = 'Force field %s must be either a scalar' %f
+            msg += ' or a vector function'
+            raise msg
+    return f    
+
 
 class Wind_stress:
@@ -938 +980 @@
 
         from config import rho_a, rho_w, eta_w
-        
-        self.s = s
-        self.phi = phi
+        from Numeric import array, Float                
+
+        self.speed = check_forcefield(s)
+        self.phi = check_forcefield(phi)
 
         self.const = eta_w*rho_a/rho_w
@@ -950 +993 @@
 
         from math import pi, cos, sin, sqrt
-        from Numeric import Float, ones 
+        from Numeric import Float, ones, ArrayType 
         
         xmom_update = domain.quantities['xmomentum'].explicit_update
@@ -958 +1001 @@
         t = domain.time        
         
-        if callable(self.s):
+        if callable(self.speed):
             xc = domain.get_centroid_coordinates()            
-            x = xc[:,0]
-            y = xc[:,1]
-            
-            s_vec = self.s(x,y,t)
+            s_vec = self.speed(xc[:,0], xc[:,1], t)
         else:
             #Assume s is a scalar
 
             try:
-                s_vec = self.s * ones(N, Float)
+                s_vec = self.speed * ones(N, Float)
             except:
                 msg = 'Speed must be either callable or a scalar: %s' %self.s
@@ -976 +1016 @@
         if callable(self.phi):
             xc = domain.get_centroid_coordinates()            
-            x = xc[:,0]
-            y = xc[:,1]
-            
-            phi_vec = self.phi(x,y,t)
+            phi_vec = self.phi(xc[:,0], xc[:,1], t)
         else:
             #Assume phi is a scalar
@@ -988 +1025 @@
                 msg = 'Angle must be either callable or a scalar: %s' %self.phi
                 raise msg
+
 
         #This is the bit that should be written in C