Changeset 955

Timestamp:

Feb 24, 2005, 5:54:54 PM (20 years ago)

Author:

ole

Message:

Unit tested File_function (NetCDF) for interpolation and time updates

Location:

inundation/ga/storm_surge/pyvolution

Files:

• data_manager.py (modified) (3 diffs)
• domain.py (modified) (3 diffs)
• shallow_water.py (modified) (1 diff)
• test_data_manager.py (modified) (1 diff)
• test_util.py (modified) (5 diffs)
• util.py (modified) (9 diffs)

inundation/ga/storm_surge/pyvolution/data_manager.py

@@ -171 +171 @@
         #self.filename = create_filename(domain.get_datadir(), 
         #       domain.get_name(), extension, len(domain))
+
         
         self.filename = create_filename(domain.get_datadir(), 
                         domain.get_name(), extension)   
-                
+
+        #print 'F', self.filename
         self.timestep = 0
         self.number_of_volumes = len(domain)
@@ -214 +216 @@
 
             fid.order = domain.default_order    
-            
+
+            #Reference point
             #Start time in seconds since the epoch (midnight 1/1/1970)
-            fid.starttime = domain.starttime 
+            fid.starttime = domain.starttime
+            fid.xllcorner = domain.xllcorner
+            fid.yllcorner = domain.yllcorner
+            fid.zone = str(domain.zone) #FIXME: ?
 
             # dimension definitions
@@ -1376 +1382 @@
     y = fid.variables['y'][:]
     stage = fid.variables['stage'][:]
-    return [min(x),max(x),min(y),max(y),min(min(stage)),max(max(stage))]
+    #print "stage",stage
+    #print "stage.shap",stage.shape
+    #print "min(stage.flat), mpythonax(stage.flat)",min(stage.flat), max(stage.flat)
+    #print "min(stage)",min(stage)
+    return [min(x),max(x),min(y),max(y),min(stage.flat),max(stage.flat)]
 
 

inundation/ga/storm_surge/pyvolution/domain.py

@@ -16 +16 @@
     def __init__(self, coordinates, vertices, boundary = None,
                  conserved_quantities = None, other_quantities = None,
-                 tagged_elements = None):
+                 tagged_elements = None, zone=None, xllcorner=0, yllcorner=0):
 
         Mesh.__init__(self, coordinates, vertices, boundary, tagged_elements)
@@ -48 +48 @@
         #Create an empty list for explicit forcing terms
         self.forcing_terms = []
-
 
 
@@ -70 +69 @@
         self.finaltime = None
         self.min_timestep = self.max_timestep = 0.0
-        self.starttime = 0    #Physical starttime if any (0 is 1 Jan 1970 00:00:00)
+        self.starttime = 0 #Physical starttime if any (0 is 1 Jan 1970 00:00:00)
+        #Origin in UTM coordinates
+        self.zone = zone
+        self.xllcorner = xllcorner
+        self.yllcorner = yllcorner        
+        
         
         #Checkpointing and storage              

inundation/ga/storm_surge/pyvolution/shallow_water.py

@@ -202 +202 @@
             #Store model data, e.g. for subsequent visualisation    
             if self.store is True:
-                #self.store_timestep(['stage', 'xmomentum', 'ymomentum'])
                 self.store_timestep(self.quantities_to_be_stored)
                 

inundation/ga/storm_surge/pyvolution/test_data_manager.py

@@ -706 +706 @@
         self.domain.reduction = mean 
         self.domain.set_datadir('.')
+
       
         sww = get_dataobject(self.domain)
         sww.store_connectivity()
         sww.store_timestep('stage')
-        self.domain.time = 2. 
+        self.domain.time = 2.
+
+        #Modify stage at second timestep
+        stage = self.domain.quantities['stage'].vertex_values
+        self.domain.set_quantity('stage', stage/2)
+
         sww.store_timestep('stage')
+
         file_and_extension_name = self.domain.filename + ".sww"
         #print "file_and_extension_name",file_and_extension_name

inundation/ga/storm_surge/pyvolution/test_util.py

@@ -195 +195 @@
         
 
+
+
+    def test_spatio_temporal_file_function_time(self):
+        """Test that File function interpolates correctly
+        between given times.
+        NetCDF version (x,y,t dependency)
+        """
+
+        #Create NetCDF (sww) file to be read
+        # x: 0, 5, 10, 15
+        # y: -20, -10, 0, 10
+        # t: 0, 60, 120, ...., 1200
+        #
+        # test quantities (arbitrary but non-trivial expressions):
+        #
+        #   stage     = 3*x - y**2 + 2*t
+        #   xmomentum = exp( -((x-7)**2 + (y+5)**2)/20 ) * t**2
+        #   ymomentum = x**2 + y**2 * sin(t*pi/600)
+
+        #Nice test that may render some of the others redundant.
+        
+        import os, time
+        from config import time_format
+        from Numeric import sin, pi, exp
+        from mesh_factory import rectangular
+        from shallow_water import Domain
+        import data_manager
+
+        finaltime = 1200
+        filename = 'test_file_function'
+
+        #Create a domain to hold test grid
+
+        points, vertices, boundary =\
+                rectangular(4, 4, 15, 30, origin = (0, -20))
+
+        
+        #print 'Number of elements', len(vertices)
+        domain = Domain(points, vertices, boundary)
+        domain.smooth = False
+        domain.default_order = 2
+        domain.set_datadir('.')
+        domain.set_name(filename)
+        domain.store = True
+        domain.format = 'sww'   #Native netcdf visualisation format
+
+        #print 'E', domain.get_extent()
+        #print points
+        start = time.mktime(time.strptime('2000', '%Y'))
+        domain.starttime = start
+
+        
+        #Store structure
+        domain.initialise_storage()
+
+        #Compute artificial time steps and store
+        dt = 60  #One minute intervals
+        t = 0.0
+        while t <= finaltime:
+            #Compute quantities
+            f1 = lambda x,y: 3*x - y**2 + 2*t + 4
+            domain.set_quantity('stage', f1)
+
+            f2 = lambda x,y: x+y+t**2
+            domain.set_quantity('xmomentum', f2)
+
+            f3 = lambda x,y: x**2 + y**2 * sin(t*pi/600)
+            domain.set_quantity('ymomentum', f3)                        
+
+            #Store and advance time
+            domain.time = t
+            domain.store_timestep(domain.conserved_quantities)
+            t += dt
+
+
+        interpolation_points = [[0,-20], [1,0], [0,1], [1.1, 3.14], [10,-12.5]]
+
+
+
+        #Set domain.starttime to too early
+        domain.starttime = start - 1                
+
+        #Create file function
+        F = file_function(filename + '.sww', domain,
+                          quantities = domain.conserved_quantities,
+                          interpolation_points = interpolation_points)
+
+        #Check that FF updates fixes domain starttime
+        assert allclose(domain.starttime, start)
+
+        #Check that domain.starttime isn't updated if later
+        domain.starttime = start + 1                
+        F = file_function(filename + '.sww', domain,
+                          quantities = domain.conserved_quantities,
+                          interpolation_points = interpolation_points)
+        assert allclose(domain.starttime, start+1)                
+
+        domain.starttime = start
+
+
+        
+        #Check linear interpolation in time
+        #for id in range(len(interpolation_points)):
+        for id in [1]:
+            x = interpolation_points[id][0]
+            y = interpolation_points[id][1]            
+            
+            for i in range(20):
+                t = i*10
+                k = i%6
+            
+                if k == 0:
+                    q0 = F(t, point_id=id)                                
+                    q1 = F(t+60, point_id=id)                                
+            
+            
+                q = F(t, point_id=id)
+                assert allclose(q, (k*q1 + (6-k)*q0)/6)
+            
+ 
+        #Another check of linear interpolation in time
+        for id in range(len(interpolation_points)):
+            q60 = F(60, point_id=id)
+            q120 = F(120, point_id=id)        
+        
+            t = 90 #Halfway between 60 and 120
+            q = F(t,point_id=id)
+            assert allclose( (q120+q60)/2, q )
+
+            t = 100 #Two thirds of the way between between 60 and 120
+            q = F(t, point_id=id)
+            assert allclose(q60/3 + 2*q120/3, q)
+
+
+
+        #Check that domain.starttime isn't updated if later than file starttime but earlier
+        #than file end time
+        delta = 23
+        domain.starttime = start + delta                
+        F = file_function(filename + '.sww', domain,
+                          quantities = domain.conserved_quantities,
+                          interpolation_points = interpolation_points)
+        assert allclose(domain.starttime, start+delta)                
+
+        
+        
+       
+        #Now try interpolation with delta offset
+        for id in [1]:
+            x = interpolation_points[id][0]
+            y = interpolation_points[id][1]            
+            
+            for i in range(20):
+                t = i*10
+                k = i%6
+            
+                if k == 0:
+                    q0 = F(t-delta, point_id=id)      
+                    q1 = F(t+60-delta, point_id=id)                            
+            
+                q = F(t-delta, point_id=id)
+                assert allclose(q, (k*q1 + (6-k)*q0)/6)
+            
+        
+        os.remove(filename + '.sww')
+        
+
+
     def test_file_function_time_with_domain(self):
         """Test that File function interpolates correctly
@@ -238 +406 @@
 
         #Check that domain.starttime isn't updated if later
-        #domain.starttime = start + 1                
-        #F = file_function(filename, domain)
-        #assert allclose(domain.starttime, start+1)                
-
+        domain.starttime = start + 1                
+        F = file_function(filename, domain)
+        assert allclose(domain.starttime, start+1)                
+
+        domain.starttime = start
         
        
@@ -275 +444 @@
         between given times. No x,y dependency here.
         Use domain with a starttime later than that of file
+
+        ASCII version
         """
 
@@ -304 +475 @@
         domain = Domain(points, vertices) 
 
-        #Check that domain.starttime isn't updated if later
+        #Check that domain.starttime isn't updated if later than file starttime but earlier
+        #than file end time
         delta = 23
         domain.starttime = start + delta                
@@ -525 +697 @@
         fid.close()
 
+
     def test_xya_ascii(self):
         import time, os

inundation/ga/storm_surge/pyvolution/util.py

@@ -133 +133 @@
     """
     
-    def  __init__(self, filename, domain=None, quantities=None, interpolation_points=None, verbose = False):
+    def  __init__(self, filename, domain=None, quantities=None,
+                  interpolation_points=None, verbose = False):
         """Initialise callable object from a file.
 
@@ -171 +172 @@
         #Get first timestep
         try:
-            starttime = fid.starttime[0]
+            self.starttime = fid.starttime[0]
         except ValueError:    
             msg = 'Could not read starttime from file %s' %filename
             raise msg
 
+        #Get origin
+        self.xllcorner = fid.xllcorner[0]
+        self.yllcorner = fid.yllcorner[0]            
     
         self.number_of_values = len(quantities)
         self.fid = fid
         self.filename = filename
-        self.starttime = starttime
         self.quantities = quantities
         self.domain = domain
@@ -237 +240 @@
             
             
-            ####self.Q = zeros( (len(time),
-            self.T = time[:] 
-            self.index = 0 #Initial index       
+            self.T = time[:]  #Time assumed to be relative to starttime
+            self.index = 0    #Initial time index       
 
             self.values = {}                
@@ -268 +270 @@
                     interpol.interpolate(fid.variables[name][i,:])   
 
+            #Report
+            if verbose:
+                print '------------------------------------------------'
+                print 'File_function statistics:'
+                print '  Name: %s' %self.filename
+                print '  Reference:'
+                print '    Lower left corner: [%f, %f]'\
+                      %(self.xllcorner, self.yllcorner) 
+                print '    Start time: %f' %self.starttime
+                print '  Extent:'
+                print '    x in [%f, %f], len(x) == %d'\
+                      %(min(x.flat), max(x.flat), len(x.flat))
+                print '    y in [%f, %f], len(y) == %d'\
+                      %(min(y.flat), max(y.flat), len(y.flat))
+                print '    t in [%f, %f], len(t) == %d'\
+                      %(min(self.T), max(self.T), len(self.T))
+                print '  Quantities:'
+                for name in self.quantities:
+                    q = fid.variables[name][:].flat
+                    print '    %s in [%f, %f]' %(name, min(q), max(q))
+
+                    
+                print '  Interpolation points (xi, eta):'
+                print '    xi in [%f, %f], len(xi) == %d'\
+                      %(min(P[:,0]), max(P[:,0]), P.shape[0])
+                print '    eta in [%f, %f], len(eta) == %d'\
+                      %(min(P[:,1]), max(P[:,1]), P.shape[1])
+                print '  Interpolated quantities:'
+                for name in self.quantities:
+                    q = self.values[name][:].flat                    
+                    print '    %s at interpolation points in [%f, %f]'\
+                          %(name, min(q), max(q))                
+                
+                
+                
+
+                print '------------------------------------------------'
         else:
-            pass 
+            msg = 'File_function_NetCDF must be invoked with '
+            msg += 'a list of interpolation points'
+            raise msg
 
     def __repr__(self):
@@ -285 +326 @@
           FIXME: point_id could also be a slice     
           FIXME: One could allow arbitrary x, y coordinates 
-          (requires computation of a new interpolator) 
+          (requires computation of a new interpolator)
+          Maybe not,.,.
           FIXME: What if x and y are vectors? 
         """
@@ -291 +333 @@
         from math import pi, cos, sin, sqrt
         from Numeric import zeros, Float
+
+
+        if point_id is None:
+            msg = 'NetCDF File function needs a point_id when invoked'
+            raise msg
+
 
         #Find time tau such that
@@ -300 +348 @@
         else:
             tau = t
+
+        #print 'D start', self.domain.starttime
+        #print 'S start', self.starttime
+        #print 't', t
+        #print 'tau', tau             
                 
-        msg = 'Time interval derived from file %s (%s:%s) does not match model time: %s'\
-              %(self.filename, self.T[0], self.T[1], tau)
+        msg = 'Time interval derived from file %s [%s:%s]'\
+              %(self.filename, self.T[0], self.T[1])
+        msg += ' does not match model time: %s' %tau
+
         if tau < self.T[0]: raise msg
         if tau > self.T[-1]: raise msg        
@@ -323 +378 @@
                     (self.T[self.index+1] - self.T[self.index])
 
+
+
                 
         #Compute interpolated values
@@ -337 +394 @@
                 q[i] = Q[self.index, point_id]
 
-            
         #Return vector of interpolated values
         return q