Changeset 1664

Timestamp:

Aug 1, 2005, 12:15:19 PM (19 years ago)

Author:

ole

Message:

Played around with file_function

Location:

inundation/ga/storm_surge/pyvolution

Files:

• test_util.py (modified) (4 diffs)
• util.py (modified) (10 diffs)

inundation/ga/storm_surge/pyvolution/test_util.py

@@ -275 +275 @@
 
         
+    def test_spatio_temporal_file_function(self):
+        """Test that spatio temporal file function performs the correct
+        interpolations in both time and space
+        NetCDF version (x,y,t dependency)        
+        """
+        import time
+
+        #Create sww file of simple propagation from left to right
+        #through rectangular domain
+        from shallow_water import Domain, Dirichlet_boundary
+        from mesh_factory import rectangular
+        from Numeric import take, concatenate, reshape
+
+        #Create basic mesh and shallow water domain
+        points, vertices, boundary = rectangular(3, 3)
+        domain1 = Domain(points, vertices, boundary)
+
+        from util import mean
+        domain1.reduction = mean
+        domain1.smooth = True #NOTE: Mimic sww output where each vertex has
+                              # only one value.
+
+        domain1.default_order = 2
+        domain1.store = True
+        domain1.set_datadir('.')
+        domain1.set_name('spatio_temporal_boundary_source_%d' %(id(self)))
+        domain1.quantities_to_be_stored = ['stage', 'xmomentum', 'ymomentum']
+
+        #Bed-slope, friction and IC at vertices (and interpolated elsewhere)
+        domain1.set_quantity('elevation', 0)
+        domain1.set_quantity('friction', 0)
+        domain1.set_quantity('stage', 0)
+
+        # Boundary conditions
+        B0 = Dirichlet_boundary([0,0,0])
+        B6 = Dirichlet_boundary([0.6,0,0])
+        domain1.set_boundary({'left': B6, 'top': B6, 'right': B0, 'bottom': B0})
+        domain1.check_integrity()
+
+        finaltime = 8
+        #Evolution
+        for t in domain1.evolve(yieldstep = 0.1, finaltime = finaltime):
+            pass
+            #domain1.write_time()
+
+
+        #Now read data from sww and check
+        from Scientific.IO.NetCDF import NetCDFFile
+        filename = domain1.get_name() + '.' + domain1.format
+        fid = NetCDFFile(filename)
+
+        x = fid.variables['x'][:]
+        y = fid.variables['y'][:]
+        stage = fid.variables['stage'][:]
+        xmomentum = fid.variables['xmomentum'][:]
+        ymomentum = fid.variables['ymomentum'][:]
+        time = fid.variables['time'][:]
+
+        #Take stage vertex values at last timestep on diagonal
+        #Diagonal is identified by vertices: 0, 5, 10, 15
+
+        timestep = len(time)-1 #Last timestep
+        d_stage = reshape(take(stage[timestep, :], [0,5,10,15]), (4,1))
+        d_uh = reshape(take(xmomentum[timestep, :], [0,5,10,15]), (4,1))
+        d_vh = reshape(take(ymomentum[timestep, :], [0,5,10,15]), (4,1))
+        D = concatenate( (d_stage, d_uh, d_vh), axis=1)
+
+        #Reference interpolated values at midpoints on diagonal at
+        #this timestep are
+        r0 = (D[0] + D[1])/2
+        r1 = (D[1] + D[2])/2
+        r2 = (D[2] + D[3])/2
+
+        #And the midpoints are found now
+        Dx = take(reshape(x, (16,1)), [0,5,10,15])
+        Dy = take(reshape(y, (16,1)), [0,5,10,15])
+
+        diag = concatenate( (Dx, Dy), axis=1)
+        d_midpoints = (diag[1:] + diag[:-1])/2
+
+        #Let us see if the file function can find the correct
+        #values at the midpoints at the last timestep:
+        f = file_function(filename, domain1,
+                          interpolation_points = d_midpoints)
+
+        q = f(timestep/10., point_id=0); assert allclose(r0, q)
+        q = f(timestep/10., point_id=1); assert allclose(r1, q)
+        q = f(timestep/10., point_id=2); assert allclose(r2, q)
+
+
+        ##################
+        #Now do the same for the first timestep
+
+        timestep = 0 #First timestep
+        d_stage = reshape(take(stage[timestep, :], [0,5,10,15]), (4,1))
+        d_uh = reshape(take(xmomentum[timestep, :], [0,5,10,15]), (4,1))
+        d_vh = reshape(take(ymomentum[timestep, :], [0,5,10,15]), (4,1))
+        D = concatenate( (d_stage, d_uh, d_vh), axis=1)
+
+        #Reference interpolated values at midpoints on diagonal at
+        #this timestep are
+        r0 = (D[0] + D[1])/2
+        r1 = (D[1] + D[2])/2
+        r2 = (D[2] + D[3])/2
+
+        #Let us see if the file function can find the correct
+        #values
+        q = f(0, point_id=0); assert allclose(r0, q)
+        q = f(0, point_id=1); assert allclose(r1, q)
+        q = f(0, point_id=2); assert allclose(r2, q)
+
+
+        ##################
+        #Now do it again for a timestep in the middle
+
+        timestep = 33
+        d_stage = reshape(take(stage[timestep, :], [0,5,10,15]), (4,1))
+        d_uh = reshape(take(xmomentum[timestep, :], [0,5,10,15]), (4,1))
+        d_vh = reshape(take(ymomentum[timestep, :], [0,5,10,15]), (4,1))
+        D = concatenate( (d_stage, d_uh, d_vh), axis=1)
+
+        #Reference interpolated values at midpoints on diagonal at
+        #this timestep are
+        r0 = (D[0] + D[1])/2
+        r1 = (D[1] + D[2])/2
+        r2 = (D[2] + D[3])/2
+
+        q = f(timestep/10., point_id=0); assert allclose(r0, q)
+        q = f(timestep/10., point_id=1); assert allclose(r1, q)
+        q = f(timestep/10., point_id=2); assert allclose(r2, q)
+
+
+        ##################
+        #Now check temporal interpolation
+        #Halfway between timestep 15 and 16
+
+        timestep = 15
+        d_stage = reshape(take(stage[timestep, :], [0,5,10,15]), (4,1))
+        d_uh = reshape(take(xmomentum[timestep, :], [0,5,10,15]), (4,1))
+        d_vh = reshape(take(ymomentum[timestep, :], [0,5,10,15]), (4,1))
+        D = concatenate( (d_stage, d_uh, d_vh), axis=1)
+
+        #Reference interpolated values at midpoints on diagonal at
+        #this timestep are
+        r0_0 = (D[0] + D[1])/2
+        r1_0 = (D[1] + D[2])/2
+        r2_0 = (D[2] + D[3])/2
+
+        #
+        timestep = 16
+        d_stage = reshape(take(stage[timestep, :], [0,5,10,15]), (4,1))
+        d_uh = reshape(take(xmomentum[timestep, :], [0,5,10,15]), (4,1))
+        d_vh = reshape(take(ymomentum[timestep, :], [0,5,10,15]), (4,1))
+        D = concatenate( (d_stage, d_uh, d_vh), axis=1)
+
+        #Reference interpolated values at midpoints on diagonal at
+        #this timestep are
+        r0_1 = (D[0] + D[1])/2
+        r1_1 = (D[1] + D[2])/2
+        r2_1 = (D[2] + D[3])/2
+
+        # The reference values are
+        r0 = (r0_0 + r0_1)/2
+        r1 = (r1_0 + r1_1)/2
+        r2 = (r2_0 + r2_1)/2
+
+        q = f((timestep - 0.5)/10., point_id=0); assert allclose(r0, q)
+        q = f((timestep - 0.5)/10., point_id=1); assert allclose(r1, q)
+        q = f((timestep - 0.5)/10., point_id=2); assert allclose(r2, q)
+
+        ##################
+        #Finally check interpolation 2 thirds of the way
+        #between timestep 15 and 16
+
+        # The reference values are
+        r0 = (r0_0 + 2*r0_1)/3
+        r1 = (r1_0 + 2*r1_1)/3
+        r2 = (r2_0 + 2*r2_1)/3
+
+        #And the file function gives
+        q = f((timestep - 1.0/3)/10., point_id=0); assert allclose(r0, q)
+        q = f((timestep - 1.0/3)/10., point_id=1); assert allclose(r1, q)
+        q = f((timestep - 1.0/3)/10., point_id=2); assert allclose(r2, q)
+
+        fid.close()
+        import os
+        os.remove(filename)
+
         
 
@@ -294 +482 @@
         #   ymomentum = x**2 + y**2 * sin(t*pi/600)
 
-        #Nice test that may render some of the others redundant.
+        #NOTE: Nice test that may render some of the others redundant.
 
         import os, time
@@ -354 +542 @@
 
 
-        #Set domain.starttime to too early
+        #Deliberately set domain.starttime to too early
         domain.starttime = start - 1
 
@@ -591 +779 @@
         os.remove(filename)
 
-
-    def test_spatio_temporal_file_function(self):
-        """Test that spatio temporal file function performs the correct
-        interpolations in both time and space
-        """
-        import time
-
-        #Create sww file of simple propagation from left to right
-        #through rectangular domain
-        from shallow_water import Domain, Dirichlet_boundary
-        from mesh_factory import rectangular
-        from Numeric import take, concatenate, reshape
-
-        #Create basic mesh and shallow water domain
-        points, vertices, boundary = rectangular(3, 3)
-        domain1 = Domain(points, vertices, boundary)
-
-        from util import mean
-        domain1.reduction = mean
-        domain1.smooth = True #NOTE: Mimic sww output where each vertex has
-                              # only one value.
-
-        domain1.default_order = 2
-        domain1.store = True
-        domain1.set_datadir('.')
-        domain1.set_name('spatio_temporal_boundary_source_%d' %(id(self)))
-        domain1.quantities_to_be_stored = ['stage', 'xmomentum', 'ymomentum']
-
-        #Bed-slope, friction and IC at vertices (and interpolated elsewhere)
-        domain1.set_quantity('elevation', 0)
-        domain1.set_quantity('friction', 0)
-        domain1.set_quantity('stage', 0)
-
-        # Boundary conditions
-        B0 = Dirichlet_boundary([0,0,0])
-        B6 = Dirichlet_boundary([0.6,0,0])
-        domain1.set_boundary({'left': B6, 'top': B6, 'right': B0, 'bottom': B0})
-        domain1.check_integrity()
-
-        finaltime = 8
-        #Evolution
-        for t in domain1.evolve(yieldstep = 0.1, finaltime = finaltime):
-            pass
-            #domain1.write_time()
-
-
-        #Now read data from sww and check
-        from Scientific.IO.NetCDF import NetCDFFile
-        filename = domain1.get_name() + '.' + domain1.format
-        fid = NetCDFFile(filename)
-
-        x = fid.variables['x'][:]
-        y = fid.variables['y'][:]
-        stage = fid.variables['stage'][:]
-        xmomentum = fid.variables['xmomentum'][:]
-        ymomentum = fid.variables['ymomentum'][:]
-        time = fid.variables['time'][:]
-
-        #Take stage vertex values at last timestep on diagonal
-        #Diagonal is identified by vertices: 0, 5, 10, 15
-
-        timestep = len(time)-1 #Last timestep
-        d_stage = reshape(take(stage[timestep, :], [0,5,10,15]), (4,1))
-        d_uh = reshape(take(xmomentum[timestep, :], [0,5,10,15]), (4,1))
-        d_vh = reshape(take(ymomentum[timestep, :], [0,5,10,15]), (4,1))
-        D = concatenate( (d_stage, d_uh, d_vh), axis=1)
-
-        #Reference interpolated values at midpoints on diagonal at
-        #this timestep are
-        r0 = (D[0] + D[1])/2
-        r1 = (D[1] + D[2])/2
-        r2 = (D[2] + D[3])/2
-
-        #And the midpoints are found now
-        Dx = take(reshape(x, (16,1)), [0,5,10,15])
-        Dy = take(reshape(y, (16,1)), [0,5,10,15])
-
-        diag = concatenate( (Dx, Dy), axis=1)
-        d_midpoints = (diag[1:] + diag[:-1])/2
-
-        #Let us see if the file function can find the correct
-        #values at the midpoints at the last timestep:
-        f = file_function(filename, domain1,
-                          interpolation_points = d_midpoints)
-
-        q = f(timestep/10., point_id=0); assert allclose(r0, q)
-        q = f(timestep/10., point_id=1); assert allclose(r1, q)
-        q = f(timestep/10., point_id=2); assert allclose(r2, q)
-
-
-        ##################
-        #Now do the same for the first timestep
-
-        timestep = 0 #First timestep
-        d_stage = reshape(take(stage[timestep, :], [0,5,10,15]), (4,1))
-        d_uh = reshape(take(xmomentum[timestep, :], [0,5,10,15]), (4,1))
-        d_vh = reshape(take(ymomentum[timestep, :], [0,5,10,15]), (4,1))
-        D = concatenate( (d_stage, d_uh, d_vh), axis=1)
-
-        #Reference interpolated values at midpoints on diagonal at
-        #this timestep are
-        r0 = (D[0] + D[1])/2
-        r1 = (D[1] + D[2])/2
-        r2 = (D[2] + D[3])/2
-
-        #Let us see if the file function can find the correct
-        #values
-        q = f(0, point_id=0); assert allclose(r0, q)
-        q = f(0, point_id=1); assert allclose(r1, q)
-        q = f(0, point_id=2); assert allclose(r2, q)
-
-
-        ##################
-        #Now do it again for a timestep in the middle
-
-        timestep = 33
-        d_stage = reshape(take(stage[timestep, :], [0,5,10,15]), (4,1))
-        d_uh = reshape(take(xmomentum[timestep, :], [0,5,10,15]), (4,1))
-        d_vh = reshape(take(ymomentum[timestep, :], [0,5,10,15]), (4,1))
-        D = concatenate( (d_stage, d_uh, d_vh), axis=1)
-
-        #Reference interpolated values at midpoints on diagonal at
-        #this timestep are
-        r0 = (D[0] + D[1])/2
-        r1 = (D[1] + D[2])/2
-        r2 = (D[2] + D[3])/2
-
-        q = f(timestep/10., point_id=0); assert allclose(r0, q)
-        q = f(timestep/10., point_id=1); assert allclose(r1, q)
-        q = f(timestep/10., point_id=2); assert allclose(r2, q)
-
-
-        ##################
-        #Now check temporal interpolation
-        #Halfway between timestep 15 and 16
-
-        timestep = 15
-        d_stage = reshape(take(stage[timestep, :], [0,5,10,15]), (4,1))
-        d_uh = reshape(take(xmomentum[timestep, :], [0,5,10,15]), (4,1))
-        d_vh = reshape(take(ymomentum[timestep, :], [0,5,10,15]), (4,1))
-        D = concatenate( (d_stage, d_uh, d_vh), axis=1)
-
-        #Reference interpolated values at midpoints on diagonal at
-        #this timestep are
-        r0_0 = (D[0] + D[1])/2
-        r1_0 = (D[1] + D[2])/2
-        r2_0 = (D[2] + D[3])/2
-
-        #
-        timestep = 16
-        d_stage = reshape(take(stage[timestep, :], [0,5,10,15]), (4,1))
-        d_uh = reshape(take(xmomentum[timestep, :], [0,5,10,15]), (4,1))
-        d_vh = reshape(take(ymomentum[timestep, :], [0,5,10,15]), (4,1))
-        D = concatenate( (d_stage, d_uh, d_vh), axis=1)
-
-        #Reference interpolated values at midpoints on diagonal at
-        #this timestep are
-        r0_1 = (D[0] + D[1])/2
-        r1_1 = (D[1] + D[2])/2
-        r2_1 = (D[2] + D[3])/2
-
-        # The reference values are
-        r0 = (r0_0 + r0_1)/2
-        r1 = (r1_0 + r1_1)/2
-        r2 = (r2_0 + r2_1)/2
-
-        q = f((timestep - 0.5)/10., point_id=0); assert allclose(r0, q)
-        q = f((timestep - 0.5)/10., point_id=1); assert allclose(r1, q)
-        q = f((timestep - 0.5)/10., point_id=2); assert allclose(r2, q)
-
-        ##################
-        #Finally check interpolation 2 thirds of the way
-        #between timestep 15 and 16
-
-        # The reference values are
-        r0 = (r0_0 + 2*r0_1)/3
-        r1 = (r1_0 + 2*r1_1)/3
-        r2 = (r2_0 + 2*r2_1)/3
-
-        #And the file function gives
-        q = f((timestep - 1.0/3)/10., point_id=0); assert allclose(r0, q)
-        q = f((timestep - 1.0/3)/10., point_id=1); assert allclose(r1, q)
-        q = f((timestep - 1.0/3)/10., point_id=2); assert allclose(r2, q)
-
-        fid.close()
-        import os
-        os.remove(filename)
 
 

inundation/ga/storm_surge/pyvolution/util.py

@@ -17 +17 @@
     v2 = v[1]/l
 
+    theta = acos(v1)
     try:
         theta = acos(v1)
@@ -25 +26 @@
         # why is it happening?
         # is it catching something we should avoid?
+        #FIXME (Ole): When did this happen? We need a unit test to
+        #reveal this condition or otherwise remove the hack
+        
         s = 1e-6
         if (v1+s >  1.0) and (v1-s < 1.0) :
@@ -154 +158 @@
     #Choose format
     #FIXME: Maybe these can be merged later on
-    #FIXME: We shoulp probably not waste enrgy on ASCII files -
+    #FIXME: We shoulp probably not waste energy on ASCII files -
     #the different formats are a pain. E.g. the time field.
     #
@@ -182 +186 @@
     interpolation_points decides at which points interpolated
     quantities are to be computed whenever object is called.
-
-
-
     If None, return average value
     """
@@ -213 +214 @@
 
         if quantities is None or len(quantities) < 1:
-            msg = 'ERROR: File must contain at least one independent value'
+            msg = 'ERROR: At least one independent value must be specified'
             raise msg
 
@@ -254 +255 @@
                 if verbose: print msg
                 domain.starttime = self.starttime #Modifying model time
-                if verbose: print 'Domain starttime is now set to %f' %domain.starttime
-
-        #Read all data in and produce values for desired data points at each timestep
+                if verbose: print 'Domain starttime is now set to %f'\
+                   %domain.starttime
+
+        #Read all data in and produce values for desired data points at
+        #each timestep
         self.spatial_interpolation(interpolation_points, verbose = verbose)
         fid.close()
 
     def spatial_interpolation(self, interpolation_points, verbose = False):
-        """For each timestep in fid: read surface, interpolate to desired points
-        and store values for use when object is called.
+        """For each timestep in fid: read surface,
+        interpolate to desired points and store values for use when
+        object is called.
         """
 
@@ -322 +326 @@
             for i, t in enumerate(self.T):
                 #Interpolate quantities at this timestep
-                #print ' time step %d of %d' %(i, len(self.T))
+                if verbose: print ' time step %d of %d' %(i, len(self.T))
                 for name in self.quantities:
                     self.values[name][i, :] =\
@@ -332 +336 @@
                 print 'File_function statistics:'
                 print '  Name: %s' %self.filename
-                print '  Reference:'
+                print '  References:'
+                #print '    Datum ....' #FIXME
                 print '    Lower left corner: [%f, %f]'\
                       %(self.xllcorner, self.yllcorner)
-                print '    Start time (file):   %f' %self.starttime
+                print '    Start time:   %f' %self.starttime                
+                #print '    Start time (file):   %f' %self.starttime
                 #print '    Start time (domain): %f' %domain.starttime
                 print '  Extent:'
@@ -359 +365 @@
                     print '    %s at interpolation points in [%f, %f]'\
                           %(name, min(q), max(q))
-
-
-
-
                 print '------------------------------------------------'
         else:
@@ -368 +370 @@
             msg += 'a list of interpolation points'
             raise msg
+            #FIXME: Should not be an error.
+            #__call__ could return an average or in case of time series only
+            #no point s are needed
+
 
     def __repr__(self):