Changeset 750

Timestamp:

Jan 20, 2005, 11:18:13 AM (19 years ago)

Author:

ole

Message:

Test that initial values are indeed output at t == 0 (just to make sure).
Also a check that bounsry is created before starting evolve.

Location:

inundation/ga/storm_surge/pyvolution

Files:

• domain.py (modified) (2 diffs)
• shallow_water.py (modified) (5 diffs)
• test_data_manager.py (modified) (4 diffs)
• test_shallow_water.py (modified) (5 diffs)
• util.py (modified) (3 diffs)

inundation/ga/storm_surge/pyvolution/domain.py

@@ -262 +262 @@
             msg = 'Conserved quantities must be a subset of all quantities'
             assert quantity in self.quantities, msg
+
+        ##assert hasattr(self, 'boundary_objects')    
     
     def write_time(self):
@@ -313 +315 @@
         from config import min_timestep, max_timestep, epsilon
 
-
+        #FIXME: Maybe lump into a larger check prior to evolving
+        msg = 'Boundary tags must be bound to boundary objects before evolving system, '
+        msg += 'e.g. using the method set_boundary.\n'        
+        msg += 'This system has the boundary tags %s ' %self.get_boundary_tags()
+        assert hasattr(self, 'boundary_objects'), msg
+        
         ##self.set_defaults()
 

inundation/ga/storm_surge/pyvolution/shallow_water.py

@@ -196 +196 @@
 
         #Store vertices and connectivity
-        self.writer.store_connectivity()                    
+        self.writer.store_connectivity()
+
 
     def store_timestep(self, name):
@@ -206 +207 @@
         self.writer.store_timestep(name)
         
+
 #Rotation of momentum vector
 def rotate(q, normal, direction = 1):
@@ -949 +951 @@
     """
 
-    def __init__(self, s, phi):
+    def __init__(self, *args, **kwargs):
         """Initialise windfield from wind speed s [m/s]
         and wind direction phi [degrees]
@@ -983 +985 @@
         Alternatively, one vector valued function for (speed, angle)
         can be applied, providing both quantities simultaneously.
-        FIXME: Not yet implemented
+        As in
+        W = Wind_stress(F), where returns (speed, angle) for each t.         
 
         domain.forcing_terms.append(W)
@@ -991 +994 @@
         from Numeric import array, Float                
 
+        if len(args) == 2:
+            s = args[0]
+            phi = args[1]
+        elif len(args) == 1:
+            #Assume vector function returning (s, phi)(t,x,y)
+            vector_function = args[0]
+            s = lambda t,x,y: vector_function(t,x,y)[0]
+            phi = lambda t,x,y: vector_function(t,x,y)[1]
+        else:
+           #Assume info is in 2 keyword arguments
+        
+           if len(kwargs) == 2:
+               s = kwargs['s']
+               phi = kwargs['phi']
+           else:
+               raise 'Assumes two keyword arguments: s=..., phi=....'   
+            
         self.speed = check_forcefield(s)
         self.phi = check_forcefield(phi)
-
+            
         self.const = eta_w*rho_a/rho_w
         

inundation/ga/storm_surge/pyvolution/test_data_manager.py

@@ -3 +3 @@
 
 import unittest
+import copy
 from Numeric import zeros, array, allclose, Float
 from util import mean
@@ -50 +51 @@
                 
         domain.set_quantity('level', level)
+        self.initial_level = copy.copy(domain.quantities['level'].vertex_values)
 
         domain.distribute_to_vertices_and_edges()    
@@ -366 +368 @@
         #Check contents
         #Get NetCDF
-        fid = NetCDFFile(sww.filename, 'r')  #Open existing file for append
+        fid = NetCDFFile(sww.filename, 'r')
       
 
@@ -400 +402 @@
 
 
+    def test_sync(self):
+        """Test info stored at each timestep is as expected (incl initial condition)
+        """
+
+        import time, os, config
+        from Numeric import array, zeros, allclose, Float, concatenate
+        from Scientific.IO.NetCDF import NetCDFFile        
+
+        self.domain.filename = 'synctest'
+        self.domain.format = 'sww'
+        self.domain.smooth = False
+        self.domain.store = True
+
+        #Evolution
+        print
+        for t in self.domain.evolve(yieldstep = 1.0, finaltime = 4.0):
+            level = self.domain.quantities['level'].vertex_values
+
+            #Get NetCDF
+            fid = NetCDFFile(self.domain.writer.filename, 'r')
+            stage = fid.variables['stage']        
+
+            if t == 0.0:
+                assert allclose(level, self.initial_level)
+                assert allclose(stage[:], level.flat) 
+            else:
+                assert not allclose(level, self.initial_level)
+                assert not allclose(stage[:], level.flat)                 
+
+
+
+
     def test_sww_DSG(self):
         """Not a test, rather a look at the sww format

inundation/ga/storm_surge/pyvolution/test_shallow_water.py

@@ -1 +1 @@
 #!/usr/bin/env python
 
-import unittest
+import unittest, os
 from math import sqrt, pi
 
@@ -18 +18 @@
     from math import sqrt, exp, cos, pi
 
+    x = array(x)
+    y = array(y)    
+
     N = len(x)
     s = 0*x  #New array
@@ -45 +48 @@
     from math import sqrt, atan, pi
 
+    x = array(x)
+    y = array(y)    
 
     N = len(x)
@@ -649 +654 @@
             pass
 
+
+
+
+    #####################################################    
+    def test_initial_condition(self):
+        """Test that initial condition is output at time == 0
+        """
+
+        from config import g
+        import copy
+        
+        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
+
+        domain.set_quantity('elevation', slope)
+        domain.set_quantity('level', level)
+
+        initial_level = copy.copy(domain.quantities['level'].vertex_values)
+
+        domain.set_boundary({'exterior': Reflective_boundary(domain)})
+
+        #Evolution
+        for t in domain.evolve(yieldstep = 1.0, finaltime = 2.0):
+            level = domain.quantities['level'].vertex_values
+
+            if t == 0.0:
+                assert allclose(level, initial_level)
+            else:
+                assert not allclose(level, initial_level)
+                
+      
+        
+
     #####################################################    
     def test_gravity(self):
@@ -879 +935 @@
 
 
+
+    def test_windfield_from_file(self):
+        from config import rho_a, rho_w, eta_w
+        from math import pi, cos, sin, sqrt
+        from config import time_format
+        from util import File_function
+        import time
+        
+
+        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 = 7 #Take a time that is represented in file (not zero)
+
+        #Write wind stress file (ensure that domain.time is covered)
+        #Take x=1 and y=0
+        filename = 'test_windstress_from_file.txt'
+        start = time.mktime(time.strptime('2000', '%Y'))        
+        fid = open(filename, 'w')
+        dt = 1  #One second interval
+        t = 0.0
+        while t <= 10.0:
+            t_string = time.strftime(time_format, time.gmtime(t+start))
+
+            fid.write('%s, %f %f\n' %(t_string,
+                                      speed(t,[1],[0])[0],
+                                      angle(t,[1],[0])[0]))
+            t += dt
+    
+        fid.close()
+
+        #Setup wind stress
+        F = File_function(filename)
+        W = Wind_stress(F) 
+        domain.forcing_terms = []
+        domain.forcing_terms.append(W)
+        
+        domain.compute_forcing_terms()
+
+        #Compute reference solution
+        const = eta_w*rho_a/rho_w
+        
+        N = domain.number_of_elements
+
+        t = domain.time
+
+        s = speed(t,[1],[0])[0]
+        phi = angle(t,[1],[0])[0]
+        
+        #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)
+
+        for k in range(N):
+            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)
+
+        os.remove(filename)
 
     def test_wind_stress_error_condition(self):

inundation/ga/storm_surge/pyvolution/util.py

@@ -301 +301 @@
         time [DD/MM/YY hh:mm:ss], x, y, value0 value1 value2 ...
 
-    In either case, the callable obejct will return a tuple of
+    In either case, the callable object will return a tuple of
     interpolated values, one each value stated in the file.
 
@@ -442 +442 @@
         """Evaluate f(t,x,y)
 
-        FIXME: x, y dependency not yet implemented
+        FIXME: x, y dependency not yet implemented except that
+        result is a vector of same length as x and y
         """
 
@@ -473 +474 @@
                 (self.T[self.index+1] - self.T[self.index])
 
-        #Compute interpolated values for values
+        #Compute interpolated values
         q = self.Q[self.index,:] +\
             ratio*(self.Q[self.index+1,:] - self.Q[self.index,:]) 
 
         #Return vector of interpolated values
-        return q
+        if x == None and y == None:
+            return q
+        else:
+            from Numeric import ones, Float
+            #Create one constant column for each value
+            N = len(x)
+            assert len(y) == N, 'x and y must have same length'
+
+            res = []
+            for col in q:
+                res.append(col*ones(N, Float))
+            
+            return res