Changeset 1664


Ignore:
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:
2 edited

Legend:

Unmodified
Added
Removed
  • inundation/ga/storm_surge/pyvolution/test_util.py

    r1486 r1664  
    275275
    276276       
     277    def test_spatio_temporal_file_function(self):
     278        """Test that spatio temporal file function performs the correct
     279        interpolations in both time and space
     280        NetCDF version (x,y,t dependency)       
     281        """
     282        import time
     283
     284        #Create sww file of simple propagation from left to right
     285        #through rectangular domain
     286        from shallow_water import Domain, Dirichlet_boundary
     287        from mesh_factory import rectangular
     288        from Numeric import take, concatenate, reshape
     289
     290        #Create basic mesh and shallow water domain
     291        points, vertices, boundary = rectangular(3, 3)
     292        domain1 = Domain(points, vertices, boundary)
     293
     294        from util import mean
     295        domain1.reduction = mean
     296        domain1.smooth = True #NOTE: Mimic sww output where each vertex has
     297                              # only one value.
     298
     299        domain1.default_order = 2
     300        domain1.store = True
     301        domain1.set_datadir('.')
     302        domain1.set_name('spatio_temporal_boundary_source_%d' %(id(self)))
     303        domain1.quantities_to_be_stored = ['stage', 'xmomentum', 'ymomentum']
     304
     305        #Bed-slope, friction and IC at vertices (and interpolated elsewhere)
     306        domain1.set_quantity('elevation', 0)
     307        domain1.set_quantity('friction', 0)
     308        domain1.set_quantity('stage', 0)
     309
     310        # Boundary conditions
     311        B0 = Dirichlet_boundary([0,0,0])
     312        B6 = Dirichlet_boundary([0.6,0,0])
     313        domain1.set_boundary({'left': B6, 'top': B6, 'right': B0, 'bottom': B0})
     314        domain1.check_integrity()
     315
     316        finaltime = 8
     317        #Evolution
     318        for t in domain1.evolve(yieldstep = 0.1, finaltime = finaltime):
     319            pass
     320            #domain1.write_time()
     321
     322
     323        #Now read data from sww and check
     324        from Scientific.IO.NetCDF import NetCDFFile
     325        filename = domain1.get_name() + '.' + domain1.format
     326        fid = NetCDFFile(filename)
     327
     328        x = fid.variables['x'][:]
     329        y = fid.variables['y'][:]
     330        stage = fid.variables['stage'][:]
     331        xmomentum = fid.variables['xmomentum'][:]
     332        ymomentum = fid.variables['ymomentum'][:]
     333        time = fid.variables['time'][:]
     334
     335        #Take stage vertex values at last timestep on diagonal
     336        #Diagonal is identified by vertices: 0, 5, 10, 15
     337
     338        timestep = len(time)-1 #Last timestep
     339        d_stage = reshape(take(stage[timestep, :], [0,5,10,15]), (4,1))
     340        d_uh = reshape(take(xmomentum[timestep, :], [0,5,10,15]), (4,1))
     341        d_vh = reshape(take(ymomentum[timestep, :], [0,5,10,15]), (4,1))
     342        D = concatenate( (d_stage, d_uh, d_vh), axis=1)
     343
     344        #Reference interpolated values at midpoints on diagonal at
     345        #this timestep are
     346        r0 = (D[0] + D[1])/2
     347        r1 = (D[1] + D[2])/2
     348        r2 = (D[2] + D[3])/2
     349
     350        #And the midpoints are found now
     351        Dx = take(reshape(x, (16,1)), [0,5,10,15])
     352        Dy = take(reshape(y, (16,1)), [0,5,10,15])
     353
     354        diag = concatenate( (Dx, Dy), axis=1)
     355        d_midpoints = (diag[1:] + diag[:-1])/2
     356
     357        #Let us see if the file function can find the correct
     358        #values at the midpoints at the last timestep:
     359        f = file_function(filename, domain1,
     360                          interpolation_points = d_midpoints)
     361
     362        q = f(timestep/10., point_id=0); assert allclose(r0, q)
     363        q = f(timestep/10., point_id=1); assert allclose(r1, q)
     364        q = f(timestep/10., point_id=2); assert allclose(r2, q)
     365
     366
     367        ##################
     368        #Now do the same for the first timestep
     369
     370        timestep = 0 #First timestep
     371        d_stage = reshape(take(stage[timestep, :], [0,5,10,15]), (4,1))
     372        d_uh = reshape(take(xmomentum[timestep, :], [0,5,10,15]), (4,1))
     373        d_vh = reshape(take(ymomentum[timestep, :], [0,5,10,15]), (4,1))
     374        D = concatenate( (d_stage, d_uh, d_vh), axis=1)
     375
     376        #Reference interpolated values at midpoints on diagonal at
     377        #this timestep are
     378        r0 = (D[0] + D[1])/2
     379        r1 = (D[1] + D[2])/2
     380        r2 = (D[2] + D[3])/2
     381
     382        #Let us see if the file function can find the correct
     383        #values
     384        q = f(0, point_id=0); assert allclose(r0, q)
     385        q = f(0, point_id=1); assert allclose(r1, q)
     386        q = f(0, point_id=2); assert allclose(r2, q)
     387
     388
     389        ##################
     390        #Now do it again for a timestep in the middle
     391
     392        timestep = 33
     393        d_stage = reshape(take(stage[timestep, :], [0,5,10,15]), (4,1))
     394        d_uh = reshape(take(xmomentum[timestep, :], [0,5,10,15]), (4,1))
     395        d_vh = reshape(take(ymomentum[timestep, :], [0,5,10,15]), (4,1))
     396        D = concatenate( (d_stage, d_uh, d_vh), axis=1)
     397
     398        #Reference interpolated values at midpoints on diagonal at
     399        #this timestep are
     400        r0 = (D[0] + D[1])/2
     401        r1 = (D[1] + D[2])/2
     402        r2 = (D[2] + D[3])/2
     403
     404        q = f(timestep/10., point_id=0); assert allclose(r0, q)
     405        q = f(timestep/10., point_id=1); assert allclose(r1, q)
     406        q = f(timestep/10., point_id=2); assert allclose(r2, q)
     407
     408
     409        ##################
     410        #Now check temporal interpolation
     411        #Halfway between timestep 15 and 16
     412
     413        timestep = 15
     414        d_stage = reshape(take(stage[timestep, :], [0,5,10,15]), (4,1))
     415        d_uh = reshape(take(xmomentum[timestep, :], [0,5,10,15]), (4,1))
     416        d_vh = reshape(take(ymomentum[timestep, :], [0,5,10,15]), (4,1))
     417        D = concatenate( (d_stage, d_uh, d_vh), axis=1)
     418
     419        #Reference interpolated values at midpoints on diagonal at
     420        #this timestep are
     421        r0_0 = (D[0] + D[1])/2
     422        r1_0 = (D[1] + D[2])/2
     423        r2_0 = (D[2] + D[3])/2
     424
     425        #
     426        timestep = 16
     427        d_stage = reshape(take(stage[timestep, :], [0,5,10,15]), (4,1))
     428        d_uh = reshape(take(xmomentum[timestep, :], [0,5,10,15]), (4,1))
     429        d_vh = reshape(take(ymomentum[timestep, :], [0,5,10,15]), (4,1))
     430        D = concatenate( (d_stage, d_uh, d_vh), axis=1)
     431
     432        #Reference interpolated values at midpoints on diagonal at
     433        #this timestep are
     434        r0_1 = (D[0] + D[1])/2
     435        r1_1 = (D[1] + D[2])/2
     436        r2_1 = (D[2] + D[3])/2
     437
     438        # The reference values are
     439        r0 = (r0_0 + r0_1)/2
     440        r1 = (r1_0 + r1_1)/2
     441        r2 = (r2_0 + r2_1)/2
     442
     443        q = f((timestep - 0.5)/10., point_id=0); assert allclose(r0, q)
     444        q = f((timestep - 0.5)/10., point_id=1); assert allclose(r1, q)
     445        q = f((timestep - 0.5)/10., point_id=2); assert allclose(r2, q)
     446
     447        ##################
     448        #Finally check interpolation 2 thirds of the way
     449        #between timestep 15 and 16
     450
     451        # The reference values are
     452        r0 = (r0_0 + 2*r0_1)/3
     453        r1 = (r1_0 + 2*r1_1)/3
     454        r2 = (r2_0 + 2*r2_1)/3
     455
     456        #And the file function gives
     457        q = f((timestep - 1.0/3)/10., point_id=0); assert allclose(r0, q)
     458        q = f((timestep - 1.0/3)/10., point_id=1); assert allclose(r1, q)
     459        q = f((timestep - 1.0/3)/10., point_id=2); assert allclose(r2, q)
     460
     461        fid.close()
     462        import os
     463        os.remove(filename)
     464
    277465       
    278466
     
    294482        #   ymomentum = x**2 + y**2 * sin(t*pi/600)
    295483
    296         #Nice test that may render some of the others redundant.
     484        #NOTE: Nice test that may render some of the others redundant.
    297485
    298486        import os, time
     
    354542
    355543
    356         #Set domain.starttime to too early
     544        #Deliberately set domain.starttime to too early
    357545        domain.starttime = start - 1
    358546
     
    591779        os.remove(filename)
    592780
    593 
    594     def test_spatio_temporal_file_function(self):
    595         """Test that spatio temporal file function performs the correct
    596         interpolations in both time and space
    597         """
    598         import time
    599 
    600         #Create sww file of simple propagation from left to right
    601         #through rectangular domain
    602         from shallow_water import Domain, Dirichlet_boundary
    603         from mesh_factory import rectangular
    604         from Numeric import take, concatenate, reshape
    605 
    606         #Create basic mesh and shallow water domain
    607         points, vertices, boundary = rectangular(3, 3)
    608         domain1 = Domain(points, vertices, boundary)
    609 
    610         from util import mean
    611         domain1.reduction = mean
    612         domain1.smooth = True #NOTE: Mimic sww output where each vertex has
    613                               # only one value.
    614 
    615         domain1.default_order = 2
    616         domain1.store = True
    617         domain1.set_datadir('.')
    618         domain1.set_name('spatio_temporal_boundary_source_%d' %(id(self)))
    619         domain1.quantities_to_be_stored = ['stage', 'xmomentum', 'ymomentum']
    620 
    621         #Bed-slope, friction and IC at vertices (and interpolated elsewhere)
    622         domain1.set_quantity('elevation', 0)
    623         domain1.set_quantity('friction', 0)
    624         domain1.set_quantity('stage', 0)
    625 
    626         # Boundary conditions
    627         B0 = Dirichlet_boundary([0,0,0])
    628         B6 = Dirichlet_boundary([0.6,0,0])
    629         domain1.set_boundary({'left': B6, 'top': B6, 'right': B0, 'bottom': B0})
    630         domain1.check_integrity()
    631 
    632         finaltime = 8
    633         #Evolution
    634         for t in domain1.evolve(yieldstep = 0.1, finaltime = finaltime):
    635             pass
    636             #domain1.write_time()
    637 
    638 
    639         #Now read data from sww and check
    640         from Scientific.IO.NetCDF import NetCDFFile
    641         filename = domain1.get_name() + '.' + domain1.format
    642         fid = NetCDFFile(filename)
    643 
    644         x = fid.variables['x'][:]
    645         y = fid.variables['y'][:]
    646         stage = fid.variables['stage'][:]
    647         xmomentum = fid.variables['xmomentum'][:]
    648         ymomentum = fid.variables['ymomentum'][:]
    649         time = fid.variables['time'][:]
    650 
    651         #Take stage vertex values at last timestep on diagonal
    652         #Diagonal is identified by vertices: 0, 5, 10, 15
    653 
    654         timestep = len(time)-1 #Last timestep
    655         d_stage = reshape(take(stage[timestep, :], [0,5,10,15]), (4,1))
    656         d_uh = reshape(take(xmomentum[timestep, :], [0,5,10,15]), (4,1))
    657         d_vh = reshape(take(ymomentum[timestep, :], [0,5,10,15]), (4,1))
    658         D = concatenate( (d_stage, d_uh, d_vh), axis=1)
    659 
    660         #Reference interpolated values at midpoints on diagonal at
    661         #this timestep are
    662         r0 = (D[0] + D[1])/2
    663         r1 = (D[1] + D[2])/2
    664         r2 = (D[2] + D[3])/2
    665 
    666         #And the midpoints are found now
    667         Dx = take(reshape(x, (16,1)), [0,5,10,15])
    668         Dy = take(reshape(y, (16,1)), [0,5,10,15])
    669 
    670         diag = concatenate( (Dx, Dy), axis=1)
    671         d_midpoints = (diag[1:] + diag[:-1])/2
    672 
    673         #Let us see if the file function can find the correct
    674         #values at the midpoints at the last timestep:
    675         f = file_function(filename, domain1,
    676                           interpolation_points = d_midpoints)
    677 
    678         q = f(timestep/10., point_id=0); assert allclose(r0, q)
    679         q = f(timestep/10., point_id=1); assert allclose(r1, q)
    680         q = f(timestep/10., point_id=2); assert allclose(r2, q)
    681 
    682 
    683         ##################
    684         #Now do the same for the first timestep
    685 
    686         timestep = 0 #First timestep
    687         d_stage = reshape(take(stage[timestep, :], [0,5,10,15]), (4,1))
    688         d_uh = reshape(take(xmomentum[timestep, :], [0,5,10,15]), (4,1))
    689         d_vh = reshape(take(ymomentum[timestep, :], [0,5,10,15]), (4,1))
    690         D = concatenate( (d_stage, d_uh, d_vh), axis=1)
    691 
    692         #Reference interpolated values at midpoints on diagonal at
    693         #this timestep are
    694         r0 = (D[0] + D[1])/2
    695         r1 = (D[1] + D[2])/2
    696         r2 = (D[2] + D[3])/2
    697 
    698         #Let us see if the file function can find the correct
    699         #values
    700         q = f(0, point_id=0); assert allclose(r0, q)
    701         q = f(0, point_id=1); assert allclose(r1, q)
    702         q = f(0, point_id=2); assert allclose(r2, q)
    703 
    704 
    705         ##################
    706         #Now do it again for a timestep in the middle
    707 
    708         timestep = 33
    709         d_stage = reshape(take(stage[timestep, :], [0,5,10,15]), (4,1))
    710         d_uh = reshape(take(xmomentum[timestep, :], [0,5,10,15]), (4,1))
    711         d_vh = reshape(take(ymomentum[timestep, :], [0,5,10,15]), (4,1))
    712         D = concatenate( (d_stage, d_uh, d_vh), axis=1)
    713 
    714         #Reference interpolated values at midpoints on diagonal at
    715         #this timestep are
    716         r0 = (D[0] + D[1])/2
    717         r1 = (D[1] + D[2])/2
    718         r2 = (D[2] + D[3])/2
    719 
    720         q = f(timestep/10., point_id=0); assert allclose(r0, q)
    721         q = f(timestep/10., point_id=1); assert allclose(r1, q)
    722         q = f(timestep/10., point_id=2); assert allclose(r2, q)
    723 
    724 
    725         ##################
    726         #Now check temporal interpolation
    727         #Halfway between timestep 15 and 16
    728 
    729         timestep = 15
    730         d_stage = reshape(take(stage[timestep, :], [0,5,10,15]), (4,1))
    731         d_uh = reshape(take(xmomentum[timestep, :], [0,5,10,15]), (4,1))
    732         d_vh = reshape(take(ymomentum[timestep, :], [0,5,10,15]), (4,1))
    733         D = concatenate( (d_stage, d_uh, d_vh), axis=1)
    734 
    735         #Reference interpolated values at midpoints on diagonal at
    736         #this timestep are
    737         r0_0 = (D[0] + D[1])/2
    738         r1_0 = (D[1] + D[2])/2
    739         r2_0 = (D[2] + D[3])/2
    740 
    741         #
    742         timestep = 16
    743         d_stage = reshape(take(stage[timestep, :], [0,5,10,15]), (4,1))
    744         d_uh = reshape(take(xmomentum[timestep, :], [0,5,10,15]), (4,1))
    745         d_vh = reshape(take(ymomentum[timestep, :], [0,5,10,15]), (4,1))
    746         D = concatenate( (d_stage, d_uh, d_vh), axis=1)
    747 
    748         #Reference interpolated values at midpoints on diagonal at
    749         #this timestep are
    750         r0_1 = (D[0] + D[1])/2
    751         r1_1 = (D[1] + D[2])/2
    752         r2_1 = (D[2] + D[3])/2
    753 
    754         # The reference values are
    755         r0 = (r0_0 + r0_1)/2
    756         r1 = (r1_0 + r1_1)/2
    757         r2 = (r2_0 + r2_1)/2
    758 
    759         q = f((timestep - 0.5)/10., point_id=0); assert allclose(r0, q)
    760         q = f((timestep - 0.5)/10., point_id=1); assert allclose(r1, q)
    761         q = f((timestep - 0.5)/10., point_id=2); assert allclose(r2, q)
    762 
    763         ##################
    764         #Finally check interpolation 2 thirds of the way
    765         #between timestep 15 and 16
    766 
    767         # The reference values are
    768         r0 = (r0_0 + 2*r0_1)/3
    769         r1 = (r1_0 + 2*r1_1)/3
    770         r2 = (r2_0 + 2*r2_1)/3
    771 
    772         #And the file function gives
    773         q = f((timestep - 1.0/3)/10., point_id=0); assert allclose(r0, q)
    774         q = f((timestep - 1.0/3)/10., point_id=1); assert allclose(r1, q)
    775         q = f((timestep - 1.0/3)/10., point_id=2); assert allclose(r2, q)
    776 
    777         fid.close()
    778         import os
    779         os.remove(filename)
    780781
    781782
  • inundation/ga/storm_surge/pyvolution/util.py

    r1658 r1664  
    1717    v2 = v[1]/l
    1818
     19    theta = acos(v1)
    1920    try:
    2021        theta = acos(v1)
     
    2526        # why is it happening?
    2627        # is it catching something we should avoid?
     28        #FIXME (Ole): When did this happen? We need a unit test to
     29        #reveal this condition or otherwise remove the hack
     30       
    2731        s = 1e-6
    2832        if (v1+s >  1.0) and (v1-s < 1.0) :
     
    154158    #Choose format
    155159    #FIXME: Maybe these can be merged later on
    156     #FIXME: We shoulp probably not waste enrgy on ASCII files -
     160    #FIXME: We shoulp probably not waste energy on ASCII files -
    157161    #the different formats are a pain. E.g. the time field.
    158162    #
     
    182186    interpolation_points decides at which points interpolated
    183187    quantities are to be computed whenever object is called.
    184 
    185 
    186 
    187188    If None, return average value
    188189    """
     
    213214
    214215        if quantities is None or len(quantities) < 1:
    215             msg = 'ERROR: File must contain at least one independent value'
     216            msg = 'ERROR: At least one independent value must be specified'
    216217            raise msg
    217218
     
    254255                if verbose: print msg
    255256                domain.starttime = self.starttime #Modifying model time
    256                 if verbose: print 'Domain starttime is now set to %f' %domain.starttime
    257 
    258         #Read all data in and produce values for desired data points at each timestep
     257                if verbose: print 'Domain starttime is now set to %f'\
     258                   %domain.starttime
     259
     260        #Read all data in and produce values for desired data points at
     261        #each timestep
    259262        self.spatial_interpolation(interpolation_points, verbose = verbose)
    260263        fid.close()
    261264
    262265    def spatial_interpolation(self, interpolation_points, verbose = False):
    263         """For each timestep in fid: read surface, interpolate to desired points
    264         and store values for use when object is called.
     266        """For each timestep in fid: read surface,
     267        interpolate to desired points and store values for use when
     268        object is called.
    265269        """
    266270
     
    322326            for i, t in enumerate(self.T):
    323327                #Interpolate quantities at this timestep
    324                 #print ' time step %d of %d' %(i, len(self.T))
     328                if verbose: print ' time step %d of %d' %(i, len(self.T))
    325329                for name in self.quantities:
    326330                    self.values[name][i, :] =\
     
    332336                print 'File_function statistics:'
    333337                print '  Name: %s' %self.filename
    334                 print '  Reference:'
     338                print '  References:'
     339                #print '    Datum ....' #FIXME
    335340                print '    Lower left corner: [%f, %f]'\
    336341                      %(self.xllcorner, self.yllcorner)
    337                 print '    Start time (file):   %f' %self.starttime
     342                print '    Start time:   %f' %self.starttime               
     343                #print '    Start time (file):   %f' %self.starttime
    338344                #print '    Start time (domain): %f' %domain.starttime
    339345                print '  Extent:'
     
    359365                    print '    %s at interpolation points in [%f, %f]'\
    360366                          %(name, min(q), max(q))
    361 
    362 
    363 
    364 
    365367                print '------------------------------------------------'
    366368        else:
     
    368370            msg += 'a list of interpolation points'
    369371            raise msg
     372            #FIXME: Should not be an error.
     373            #__call__ could return an average or in case of time series only
     374            #no point s are needed
     375
    370376
    371377    def __repr__(self):
Note: See TracChangeset for help on using the changeset viewer.