Changeset 609

Timestamp:

Nov 22, 2004, 4:05:25 PM (20 years ago)

Author:

ole

Message:

First unit tests of wind field

Location:

inundation/ga/storm_surge/pyvolution

Files:

• shallow_water.py (modified) (1 diff)
• test_shallow_water.py (modified) (2 diffs)

inundation/ga/storm_surge/pyvolution/shallow_water.py

@@ -975 +975 @@
                     phi = self.phi                  
 
+        
+
                 #Convert to radians
                 phi = phi*pi/180
-        
+
+                
                 #Compute velocity vector (u, v)
                 u = s*cos(phi)

inundation/ga/storm_surge/pyvolution/test_shallow_water.py

@@ -7 +7 @@
 from Numeric import allclose, array, zeros, ones, Float
 from shallow_water import *
+
+
+
+#Variable windfield implemented using functions
+def speed(x,y,t):
+    """Large speeds halfway between center and edges
+    Low speeds at center and edges
+    """
+    
+    from math import sqrt, exp, cos, pi
+    
+    r = sqrt(x**2 + y**2)
+
+    factor = exp( -(r-0.15)**2 )
+
+    return 4000 * factor * (cos(t*2*pi/150) + 2)
+
+
+def angle(x,y,t):
+    """Rotating field
+    """
+    from math import sqrt, atan, pi
+
+    r = sqrt(x**2 + y**2)    
+    
+    angle = atan(y/x)
+
+    if x < 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
+
 
         
@@ -586 +625 @@
         assert allclose(domain.quantities['ymomentum'].explicit_update, 0)
         
+        
+    def test_manning_friction(self):
+        from config import g
+        
+        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)
+
+        #Set up for a gradient of (3,0) at mid triangle         
+        def slope(x, y):
+            return 3*x
+
+        h = 0.1
+        def level(x,y):
+            return slope(x,y)+h
+
+        eta = 0.07
+        domain.set_quantity('elevation', slope)
+        domain.set_quantity('level', level)
+        domain.set_quantity('friction', eta)        
+
+        for name in domain.conserved_quantities:
+            assert allclose(domain.quantities[name].explicit_update, 0)
+            assert allclose(domain.quantities[name].semi_implicit_update, 0)
+            
+        domain.compute_forcing_terms()                                    
+
+        assert allclose(domain.quantities['level'].explicit_update, 0)
+        assert allclose(domain.quantities['xmomentum'].explicit_update, -g*h*3)
+        assert allclose(domain.quantities['ymomentum'].explicit_update, 0)
+
+        assert allclose(domain.quantities['level'].semi_implicit_update, 0)
+        assert allclose(domain.quantities['xmomentum'].semi_implicit_update, 0)
+        assert allclose(domain.quantities['ymomentum'].semi_implicit_update, 0)        
+
+        #Create some momentum for friction to work with
+        domain.set_quantity('xmomentum', 1)
+        S = -g * eta**2 / h**(7.0/3)
+
+        domain.compute_forcing_terms()                                   
+        assert allclose(domain.quantities['level'].semi_implicit_update, 0)
+        assert allclose(domain.quantities['xmomentum'].semi_implicit_update, S)
+        assert allclose(domain.quantities['ymomentum'].semi_implicit_update, 0)        
+
+        #A more complex example 
+        domain.quantities['level'].semi_implicit_update[:] = 0.0
+        domain.quantities['xmomentum'].semi_implicit_update[:] = 0.0
+        domain.quantities['ymomentum'].semi_implicit_update[:] = 0.0        
+       
+        domain.set_quantity('xmomentum', 3)
+        domain.set_quantity('ymomentum', 4)        
+
+        S = -g * eta**2 * 5 / h**(7.0/3)
+
+
+        domain.compute_forcing_terms()
+
+        assert allclose(domain.quantities['level'].semi_implicit_update, 0)
+        assert allclose(domain.quantities['xmomentum'].semi_implicit_update, 3*S)
+        assert allclose(domain.quantities['ymomentum'].semi_implicit_update, 4*S)        
+
+    def test_constant_wind_stress(self):
+        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})
+
+        #Setup only one forcing term, constant wind stress
+        s = 100
+        phi = 135
+        domain.forcing_terms = []
+        domain.forcing_terms.append( Wind_stress(s, phi) )
+        
+        domain.compute_forcing_terms()
+
+
+        const = eta_w*rho_a/rho_w
+        
+        #Convert to radians
+        phi = phi*pi/180
+        
+        #Compute velocity vector (u, v)
+        u = s*cos(phi)
+        v = s*sin(phi)
+
+        #Compute wind stress
+        S = const * sqrt(u**2 + v**2)
+
+        assert allclose(domain.quantities['level'].explicit_update, 0)
+        assert allclose(domain.quantities['xmomentum'].explicit_update, S*u)
+        assert allclose(domain.quantities['ymomentum'].explicit_update, S*v)
+
+
+    def test_variable_wind_stress(self):
+        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, constant wind stress
+        s = 100
+        phi = 135
+        domain.forcing_terms = []
+        domain.forcing_terms.append( Wind_stress(s = speed, phi = angle) )
+        
+        domain.compute_forcing_terms()
+
+        #Compute reference solution
+        const = eta_w*rho_a/rho_w
+        
+        N = domain.number_of_elements
+
+        xc = domain.get_centroid_coordinates()
+        t = domain.time
+
+        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
+        
+            #Compute velocity vector (u, v)
+            u = s*cos(phi)
+            v = s*sin(phi)
+
+            #Compute wind stress
+            S = const * sqrt(u**2 + v**2)
+
+            assert allclose(domain.quantities['level'].explicit_update[k], 0)
+            assert allclose(domain.quantities['xmomentum'].explicit_update[k], S*u)
+            assert allclose(domain.quantities['ymomentum'].explicit_update[k], S*v)
+
+
+
+
+        #FIXME: Do this when I get a minute of peace - goddammit
+        
+