Changeset 612

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

Location:

inundation/ga/storm_surge/pyvolution

Files:

• shallow_water.py (modified) (6 diffs)
• test_shallow_water.py (modified) (2 diffs)
• wind_example_variable.py (modified) (4 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    

inundation/ga/storm_surge/pyvolution/test_shallow_water.py

@@ -30 +30 @@
 
     return s
+
+
+def scalar_func(x,y,t):
+    """Function that returns a scalar.
+    Used to test error message when Numeric array is expected
+    """
+
+    return 17.7
 
 
@@ -829 +837 @@
 
 
-        #FIXME: Do this when I get a minute of peace - goddammit
-        
-
+    def test_wind_stress_error_condition(self):
+        """Test that windstress reacts properly when forcing functions
+        are wrong - e.g. returns a scalar
+        """
+        
+        from config import rho_a, rho_w, eta_w
+        from math import pi, cos, sin, sqrt
+
+        a = [0.0, 0.0]
+        b = [0.0, 2.0]
+        c = [2.0, 0.0]
+        d = [0.0, 4.0]
+        e = [2.0, 2.0]
+        f = [4.0, 0.0]
+
+        points = [a, b, c, d, e, f]
+        #bac, bce, ecf, dbe
+        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
+
+        domain = Domain(points, vertices)
+
+        #Flat surface with 1m of water
+        domain.set_quantity('elevation', 0)
+        domain.set_quantity('level', 1.0)
+        domain.set_quantity('friction', 0)        
+
+        Br = Reflective_boundary(domain)
+        domain.set_boundary({'exterior': Br})
+
+
+        domain.time = 5.54 #Take a random time (not zero)
+
+        #Setup only one forcing term, bad func
+        domain.forcing_terms = []
+
+        try:
+            domain.forcing_terms.append(Wind_stress(s = scalar_func,
+                                                    phi = angle))
+        except AssertionError:
+            pass
+        else:
+            msg = 'Should have raised exception'
+            raise msg
+        
+
+        try:
+            domain.forcing_terms.append(Wind_stress(s = speed,
+                                                    phi = scalar_func))
+        except AssertionError:
+            pass
+        else:
+            msg = 'Should have raised exception'
+            raise msg
+        
+        try:
+            domain.forcing_terms.append(Wind_stress(s = speed,
+                                                    phi = 'xx'))
+        except:
+            pass
+        else:
+            msg = 'Should have raised exception'
+            raise msg
         
 

inundation/ga/storm_surge/pyvolution/wind_example_variable.py

@@ -13 +13 @@
 #Create basic mesh
 N = 80
-len = 100
-points, vertices, boundary = rectangular(N, N, len, len)
+length = 100
+points, vertices, boundary = rectangular(N, N, length, length)
 
 #Create shallow water domain
@@ -35 +35 @@
     Low speeds at center and edges
     """
+
+    from math import sqrt, exp, cos, pi
+
+    N = len(x)
+    s = 0*x  #New array
+
+    c = (length/2, length/2)
     
-    from math import sqrt, exp, cos, pi
-    
-    c = (len1/2, len2/2)
-    r = sqrt((x - c[0])**2 + (y - c[1])**2)
+    for k in range(N):
+        r = sqrt((x[k] - c[0])**2 + (y[k] - c[1])**2)
 
-    r /= max(len1,len2)
+        r /= max(length,length)
 
-    factor = exp( -(r-0.15)**2 )
+        factor = exp( -(r-0.15)**2 )
         
-    #print x,y,factor
+        s[k] = 4000 * factor * (cos(t*2*pi/150) + 2)
 
-    return 4000 * factor * (cos(t*2*pi/150) + 2)
+    return s
 
 
@@ -55 +60 @@
     from math import sqrt, atan, pi
 
-    c = (len1/2, len2/2)
-    r = sqrt((x - c[0])**2 + (y - c[1])**2)    
+    N = len(x)
+    a = 0*x  #New array
+
+    c = (length/2, length/2)
+    for k in range(N):
+        r = sqrt((x[k] - c[0])**2 + (y[k] - c[1])**2)
+
+        xx = (x[k] - c[0])/length
+        yy = (y[k] - c[1])/length
+
+        angle = atan(yy/xx)
+
+        if xx < 0:
+            angle+=pi #atan in ]-pi/2; pi/2[
+
+        #Take normal direction
+        angle -= pi/2
     
-    xx = (x - c[0])/len1
-    yy = (y - c[1])/len2
+        #Ensure positive radians    
+        if angle < 0:
+            angle += 2*pi
 
-    angle = atan(yy/xx)
-
-    if xx < 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
+        a[k] = angle/pi*180
+        
+    return a
 
     
@@ -82 +93 @@
 def gust(x,y,t): 
     from math import sin, pi
+    from Numeric import zeros, ones, Float
+
+    N = len(x)
 
     tt = sin(2*pi*t/200)
 
     if tt > 0.9:
-        return 6000*tt
+        return 6000*tt*ones(N, Float)
     else:
-        return 0
+        return zeros(N, Float)
     
 domain.forcing_terms.append(Wind_stress(gust, 25))