Changeset 3900


Ignore:
Timestamp:
Nov 1, 2006, 10:58:50 AM (16 years ago)
Author:
ole
Message:

Added time_thinning to Interpolation_function and improved diagnostics.

Location:
anuga_core/source/anuga
Files:
3 edited

Legend:

Unmodified
Added
Removed
  • anuga_core/source/anuga/abstract_2d_finite_volumes/util.py

    r3850 r3900  
    1313
    1414def file_function(filename,
    15                   domain = None,
    16                   quantities = None,
    17                   interpolation_points = None,
    18                   verbose = False,
    19                   use_cache = False):
     15                  domain=None,
     16                  quantities=None,
     17                  interpolation_points=None,
     18                  time_thinning=1,
     19                  verbose=False,
     20                  use_cache=False):
    2021    """Read time history of spatial data from NetCDF file and return
    2122    a callable object.
     
    5859
    5960
     61    # Build arguments and keyword arguments for use with caching or apply.
     62    args = (filename,)
     63   
     64    kwargs = {'domain': domain,
     65              'quantities': quantities,
     66              'interpolation_points': interpolation_points,
     67              'time_thinning': time_thinning,                   
     68              'verbose': verbose}
     69
     70
     71    # Call underlying engine with or without caching
    6072    if use_cache is True:
    6173        try:
     
    6678            raise msg
    6779
    68 
    6980        f = cache(_file_function,
    70                   filename,
    71                   {'domain': domain,
    72                    'quantities': quantities,
    73                    'interpolation_points': interpolation_points,
    74                    'verbose': verbose},
    75                   dependencies = [filename],
    76                   compression = False,
    77                   verbose = verbose)
    78         #FIXME (Ole): Pass cache arguments, such as compression, in some sort of
    79         #structure
     81                  args, kwargs,
     82                  dependencies=[filename],
     83                  compression=False,                 
     84                  verbose=verbose)
     85
     86    else:
     87        f = apply(_file_function,
     88                  args, kwargs)
     89
     90
     91    #FIXME (Ole): Pass cache arguments, such as compression, in some sort of
     92    #structure
    8093       
    81     else:
    82         f = _file_function(filename,
    83                            domain,
    84                            quantities,
    85                            interpolation_points,
    86                            verbose)           
    8794
    8895    return f
     
    9198
    9299def _file_function(filename,
    93                    domain = None,
    94                    quantities = None,
    95                    interpolation_points = None,
    96                    verbose = False):
     100                   domain=None,
     101                   quantities=None,
     102                   interpolation_points=None,
     103                   time_thinning=1,                                               
     104                   verbose=False):
    97105    """Internal function
    98106   
     
    133141        return get_netcdf_file_function(filename, domain, quantities,
    134142                                        interpolation_points,
     143                                        time_thinning=time_thinning,
    135144                                        verbose = verbose)
    136145    else:
     
    143152                             quantity_names=None,
    144153                             interpolation_points=None,
     154                             time_thinning=1,                             
    145155                             verbose = False):
    146156    """Read time history of spatial data from NetCDF sww file and
     
    238248    time = fid.variables['time'][:]   
    239249
     250    # Get time independent stuff
    240251    if spatial:
    241252        #Get origin
     
    296307       
    297308   
    298     #Produce values for desired data points at
    299     #each timestep
    300 
     309    # Produce values for desired data points at
     310    # each timestep for each quantity
    301311    quantities = {}
    302312    for i, name in enumerate(quantity_names):
     
    318328                                  triangles,
    319329                                  interpolation_points,
     330                                  time_thinning=time_thinning,
    320331                                  verbose=verbose)
    321332
  • anuga_core/source/anuga/fit_interpolate/interpolate.py

    r3896 r3900  
    405405    FIXME (Ole): Need to allow vertex coordinates and interpolation points to be
    406406    geospatial data objects
     407
     408    Time assumed to be relative to starttime   
    407409    """
    408410 
     
    415417                 triangles=None,
    416418                 interpolation_points=None,
     419                 time_thinning=1,
    417420                 verbose=False):
    418421        """Initialise object and build spatial interpolation if required
     422
     423        Time_thinning_number controls how many timesteps to use. Only timesteps with
     424        index%time_thinning_number == 0 will used, or in other words a value of 3, say,
     425        will cause the algorithm to use every third time step.
    419426        """
    420427
     
    458465            self.spatial = True           
    459466
     467        # Thin timesteps if needed
     468        # Note array() is used to make the thinned arrays contiguous in memory
     469        self.time = array(time[::time_thinning])         
     470        for name in quantity_names:
     471            if len(quantities[name].shape) == 2:
     472                quantities[name] = array(quantities[name][::time_thinning,:])
    460473             
    461474        # Save for use with statistics
     
    468481        self.vertex_coordinates = vertex_coordinates
    469482        self.interpolation_points = interpolation_points
    470         self.time = time[:]  # Time assumed to be relative to starttime
     483       
     484
    471485        self.index = 0    # Initial time index
    472486        self.precomputed_values = {}
     
    495509
    496510            # Build interpolator
    497             if verbose: print 'Building interpolation matrix'           
     511            if verbose:
     512                msg = 'Building interpolation matrix from source mesh '
     513                msg += '(%d vertices, %d triangles)' %(vertex_coordinates.shape[0],
     514                                                       triangles.shape[0])
     515                print msg
     516
     517               
    498518            interpol = Interpolate(vertex_coordinates,
    499519                                   triangles,
    500                                    #point_coordinates = \
    501                                    #self.interpolation_points,
    502                                    #alpha = 0,
    503520                                   verbose=verbose)
    504521
    505             if verbose: print 'Interpolate'
     522            if verbose:
     523                print 'Interpolating (%d interpolation points, %d timesteps).'\
     524                      %(self.interpolation_points.shape[0], self.time.shape[0]),
     525           
     526                if time_thinning > 1:
     527                    print 'Timesteps were thinned by a factor of %d' %time_thinning
     528                else:
     529                    print
     530
    506531            for i, t in enumerate(self.time):
    507532                # Interpolate quantities at this timestep
  • anuga_core/source/anuga/fit_interpolate/test_interpolate.py

    r3846 r3900  
    11241124        else:
    11251125            raise 'Should raise exception'
     1126
     1127
     1128
     1129    def test_interpolation_interface_with_time_thinning(self):
     1130        # Test spatio-temporal interpolation
     1131        # Test that spatio temporal function performs the correct
     1132        # interpolations in both time and space
     1133   
     1134        #Three timesteps
     1135        time = [1.0, 2.0, 4.0, 5.0, 7.0, 8.0, 9.0, 10.0]   
     1136
     1137        #Setup mesh used to represent fitted function
     1138        a = [0.0, 0.0]
     1139        b = [0.0, 2.0]
     1140        c = [2.0, 0.0]
     1141        d = [0.0, 4.0]
     1142        e = [2.0, 2.0]
     1143        f = [4.0, 0.0]
     1144
     1145        points = [a, b, c, d, e, f]
     1146        #bac, bce, ecf, dbe
     1147        triangles = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
     1148
     1149
     1150        #New datapoints where interpolated values are sought
     1151        interpolation_points = [[ 0.0, 0.0],
     1152                                [ 0.5, 0.5],
     1153                                [ 0.7, 0.7],
     1154                                [ 1.0, 0.5],
     1155                                [ 2.0, 0.4],
     1156                                [ 2.8, 1.2]]
     1157
     1158        #One quantity
     1159        Q = zeros( (8,6), Float )
     1160
     1161        #Linear in time and space
     1162        for i, t in enumerate(time):
     1163            Q[i, :] = t*linear_function(points)
     1164
     1165        # Check interpolation of one quantity using interpolaton points) using default
     1166        # time_thinning of 1
     1167        I = Interpolation_function(time, Q,
     1168                                   vertex_coordinates=points,
     1169                                   triangles=triangles,
     1170                                   interpolation_points=interpolation_points,
     1171                                   verbose=False)
     1172
     1173        answer = linear_function(interpolation_points)
     1174
     1175       
     1176        t = time[0]
     1177        for j in range(50): #t in [1, 6]
     1178            for id in range(len(interpolation_points)):
     1179                assert allclose(I(t, id), t*answer[id])
     1180            t += 0.1   
     1181
     1182
     1183        # Now check time_thinning
     1184        I = Interpolation_function(time, Q,
     1185                                   vertex_coordinates=points,
     1186                                   triangles=triangles,
     1187                                   interpolation_points=interpolation_points,
     1188                                   time_thinning=2,
     1189                                   verbose=False)
     1190
     1191
     1192        assert len(I.time) == 4
     1193        assert( allclose(I.time, [1.0, 4.0, 7.0, 9.0] ))   
     1194
     1195        answer = linear_function(interpolation_points)
     1196
     1197        t = time[0]
     1198        for j in range(50): #t in [1, 6]
     1199            for id in range(len(interpolation_points)):
     1200                assert allclose(I(t, id), t*answer[id])
     1201            t += 0.1   
     1202
     1203
    11261204
    11271205
Note: See TracChangeset for help on using the changeset viewer.