Changeset 849


Ignore:
Timestamp:
Feb 8, 2005, 2:55:56 PM (20 years ago)
Author:
ole
Message:

Work towards spatio-temporal boundary
Also changed naming from 'z' to 'elevation' in data manager (kept 'z' for backwards compatibility, though)

Location:
inundation/ga/storm_surge/pyvolution
Files:
11 edited

Legend:

Unmodified
Added
Removed

  • inundation/ga/storm_surge/pyvolution/data_manager.py

    r848 r849  
    229229            fid.createVariable('x', self.precision, ('number_of_points',))
    230230            fid.createVariable('y', self.precision, ('number_of_points',))
    231             fid.createVariable('z', self.precision, ('number_of_points',))
     231            fid.createVariable('elevation', self.precision, ('number_of_points',))
     232           
     233            #FIXME: Backwards compatibility
     234            fid.createVariable('z', self.precision, ('number_of_points',))
     235            #################################
    232236       
    233237            fid.createVariable('volumes', Int, ('number_of_volumes',
    … …  
    272276        x = fid.variables['x']
    273277        y = fid.variables['y']
    274         z = fid.variables['z']
     278        z = fid.variables['elevation']
    275279       
    276280        volumes = fid.variables['volumes']
    … …  
    286290        y[:] = Y.astype(self.precision)
    287291        z[:] = Z.astype(self.precision)
     292       
     293        #FIXME: Backwards compatibility
     294        z = fid.variables['z']
     295        z[:] = Z.astype(self.precision)
     296        ################################
    288297
    289298        volumes[:] = V.astype(volumes.typecode())
    … …  
    403412            fid.createVariable('x', self.precision, ('number_of_points',))
    404413            fid.createVariable('y', self.precision, ('number_of_points',))
    405             #fid.createVariable('z', self.precision, ('number_of_points',))
     414
    406415       
    407416            fid.createVariable('volumes', Int, ('number_of_volumes',
    … …  
    648657    x = fid.variables['x']
    649658    y = fid.variables['y']
    650     z = fid.variables['z']       
     659    z = fid.variables['elevation']       
    651660    time = fid.variables['time']
    652661    stage = fid.variables['stage']

  • inundation/ga/storm_surge/pyvolution/general_mesh.py

    r715 r849  
    3737   
    3838        #creates two triangles: bac and bce
    39        
     39
     40    Other:
     41   
     42      In addition mesh computes an Nx6 array called vertex_coordinates.
     43      This structure is derived from coordinates and contains for each
     44      triangle the three x,y coordinates at the vertices.
     45                               
    4046
    4147        This is a cut down version of mesh from pyvolution\mesh.py

  • inundation/ga/storm_surge/pyvolution/generic_boundary_conditions.py

    r844 r849  
    175175        from Numeric import array
    176176        from config import time_format
    177         from util import File_function
    178        
    179         Boundary.__init__(self)
    180 
    181         self.F = File_function(filename, domain)
     177        from util import file_function
     178       
     179        Boundary.__init__(self)
     180
     181        self.F = file_function(filename, domain)
    182182        self.domain = domain
    183183
    … …  
    205205
    206206        d = len(domain.conserved_quantities)
    207         msg = 'Values specified in file must be a list or an array of length %d' %d
     207        msg = 'Values specified in file %s must be a list or an array of length %d' %(filename, d)
    208208        assert len(q) == d, msg
    209209

  • inundation/ga/storm_surge/pyvolution/interpolate_sww.py

    r754 r849  
    106106        try:
    107107            if quantity_name == 'height':
    108                 z = fid.variables['z'][:]
     108                z = fid.variables['elevation'][:]
    109109                stage = fid.variables['stage'][:,:]
    110110                quantity = stage - z  # 2D, using broadcasting

  • inundation/ga/storm_surge/pyvolution/mesh.py

    r655 r849  
    409409                break
    410410
    411             #Should this warning be an exception?
    412             #assert v is not None, msg
    413411
    414412               
    … …  
    453451        #for id, edge in self.boundary:
    454452        #    assert self.neighbours[id,edge] < 0
    455 
    456 
     453        #
     454        #NOTE (Ole): I reckon this was resolved late 2004?
     455        #
     456        #See domain.set_boundary
    457457       
    458458           

  • inundation/ga/storm_surge/pyvolution/shallow_water.py

    r845 r849  
    2222
    2323symbol    variable name    explanation
    24 z         elevation        elevation of bed on which flow is modelled
    25 h         height           water height above z
    26 w         stage            water level, w = z+h
    27 u                          speed in the x direction
    28 v                          speed in the y direction
    29 uh        xmomentum        momentum in the x direction
    30 vh        ymomentum        momentum in the y direction
    31 
    32 eta                        mannings friction coefficient
    33 nu                         wind stress coefficient
     24x         x                horizontal distance from origin [m]
     25y         y                vertical distance from origin [m]
     26z         elevation        elevation of bed on which flow is modelled [m]
     27h         height           water height above z [m]
     28w         stage            water level, w = z+h [m]
     29u                          speed in the x direction [m/s]
     30v                          speed in the y direction [m/s]
     31uh        xmomentum        momentum in the x direction [m^2/s]
     32vh        ymomentum        momentum in the y direction [m^2/s]
     33
     34eta                        mannings friction coefficient []
     35nu                         wind stress coefficient []
    3436
    3537The conserved quantities are w, uh, vh
    … …  
    111113        assert self.conserved_quantities[2] == 'ymomentum', msg
    112114       
    113 
    114         #Check that stages are >= bed elevation
    115         from Numeric import alltrue, greater_equal
    116        
    117         stage = self.quantities['stage']
    118         bed = self.quantities['elevation']       
    119 
    120115
    121116    def compute_fluxes(self):
    … …  
    766761
    767762       
    768 class Spatio_temporal_boundary(Boundary):
    769     """The spatio-temporal boundary, reads values for the conserved
    770     quantities from an sww NetCDF file, and returns interpolated values
    771     at the midpoints of each associated boundaty segment.
    772     Time dependency is interpolated linearly as in util.File_function.
    773      
    774     """
    775 
    776     pass
     763#class Spatio_temporal_boundary(Boundary):
     764#    """The spatio-temporal boundary, reads values for the conserved
     765#    quantities from an sww NetCDF file, and returns interpolated values
     766#    at the midpoints of each associated boundaty segment.
     767#    Time dependency is interpolated linearly as in util.File_function.#
     768#
     769#    Example:
     770#    Bf = Spatio_temporal_boundary('source_file.sww', domain)
     771#         
     772#    """
     773Spatio_temporal_boundary = File_boundary
     774       
     775
    777776       
    778777

  • inundation/ga/storm_surge/pyvolution/test_data_manager.py

    r839 r849  
    149149        x = fid.variables['x']
    150150        y = fid.variables['y']
    151         z = fid.variables['z']
     151        z = fid.variables['elevation']
    152152
    153153        volumes = fid.variables['volumes']       
    … …  
    196196        x = fid.variables['x']
    197197        y = fid.variables['y']
    198         z = fid.variables['z']
     198        z = fid.variables['elevation']
    199199
    200200        volumes = fid.variables['volumes']       
    … …  
    256256        x = fid.variables['x']
    257257        y = fid.variables['y']
    258         z = fid.variables['z']       
     258        z = fid.variables['elevation']       
    259259        time = fid.variables['time']
    260260        stage = fid.variables['stage']
    … …  
    314314        x = fid.variables['x']
    315315        y = fid.variables['y']
    316         z = fid.variables['z']
     316        z = fid.variables['elevation']
    317317        time = fid.variables['time']
    318318        stage = fid.variables['stage']       
    … …  
    374374        x = fid.variables['x']
    375375        y = fid.variables['y']
    376         z = fid.variables['z']
     376        z = fid.variables['elevation']
    377377        time = fid.variables['time']
    378378        stage = fid.variables['stage']       
    … …  
    459459        x = fid.variables['x']
    460460        y = fid.variables['y']
    461         z = fid.variables['z']
     461        z = fid.variables['elevation']
    462462
    463463        volumes = fid.variables['volumes']   

  • inundation/ga/storm_surge/pyvolution/test_interpolate_sww.py

    r773 r849  
    9898        x = fid.variables['x']
    9999        y = fid.variables['y']
    100         z = fid.variables['z']
     100        z = fid.variables['elevation']
    101101
    102102        volumes = fid.variables['volumes']   
    … …  
    218218        # look at error catching
    219219        try:
    220             interp = Interpolate_sww(sww.filename, 'z')
     220            interp = Interpolate_sww(sww.filename, 'elevation')
    221221        except KeyError:
    222222            pass

  • inundation/ga/storm_surge/pyvolution/test_shallow_water.py

    r845 r849  
    55
    66from config import g, epsilon
    7 from Numeric import allclose, array, zeros, ones, Float
     7from Numeric import allclose, array, zeros, ones, Float, take
    88from shallow_water import *
    99
    … …  
    949949        from math import pi, cos, sin, sqrt
    950950        from config import time_format
    951         from util import File_function
     951        from util import file_function
    952952        import time
    953953       
    … …  
    995995
    996996        #Setup wind stress
    997         F = File_function(filename)
     997        F = file_function(filename)
    998998        W = Wind_stress(F)
    999999        domain.forcing_terms = []
    … …  
    23802380       
    23812381        #Create sww file of simple propagation from left to right
    2382         #through rectangular domain: domain1
     2382        #through rectangular domain
    23832383       
    23842384        from mesh_factory import rectangular
    … …  
    23882388
    23892389        #Create shallow water domain
    2390         domain = Domain(points, vertices, boundary)
    2391         domain.smooth = False
    2392 
    2393         domain.visualise = False
    2394         domain.default_order = 2
    2395 
    2396         domain.set_datadir('.')
    2397         domain.store = True     
    2398         domain.set_name('test_spatio_temporal_boundary_' + str(time.time()))
    2399        
     2390        domain1 = Domain(points, vertices, boundary)
     2391       
     2392        from util import mean
     2393        domain1.reduction = mean
     2394        domain1.smooth = True #FIXME: This may be a requirement for the
     2395                             #interpolation process
     2396
     2397        #domain1.visualise = True
     2398        domain1.default_order = 2
     2399
     2400        domain1.store = True   
     2401        domain1.set_datadir('.')
     2402        domain1.set_name('spatio_temporal_boundary_source' + str(time.time()))
     2403
     2404               
    24002405        #Bed-slope and friction at vertices (and interpolated elsewhere)
    2401         domain.set_quantity('elevation', 0)
    2402         domain.set_quantity('friction', 0)     
     2406        domain1.set_quantity('elevation', 0)
     2407        domain1.set_quantity('friction', 0)     
    24032408
    24042409        # Boundary conditions
    2405         Br = Reflective_boundary(domain)
     2410        Br = Reflective_boundary(domain1)
    24062411        Bd = Dirichlet_boundary([0.3,0,0])     
    2407         Bt = Transmissive_boundary(domain)
    2408         domain.set_boundary({'left': Bd, 'right': Bt, 'top': Br, 'bottom': Br})
     2412        Bt = Transmissive_boundary(domain1)
     2413        domain1.set_boundary({'left': Bd, 'right': Bt, 'top': Br, 'bottom': Br})
    24092414
    24102415        #Initial condition
    2411         domain.set_quantity('stage', 0)
    2412         domain.check_integrity()
     2416        domain1.set_quantity('stage', 0)
     2417        domain1.check_integrity()
    24132418       
    24142419        #Evolution
    2415         for t in domain.evolve(yieldstep = 0.05, finaltime = 4):
    2416             import time
    2417             #time.sleep(0.01)
    2418             #domain.write_time()           
    2419 
    2420 
    2421        
    2422        
    2423        
    2424         #Create an triangle shaped domain, formed from the
    2425         #lower and right hand  boundaries and the sw-ne diagonal
    2426         #from domain 1. Called it domain2
     2420        for t in domain1.evolve(yieldstep = 0.05, finaltime = 4):
     2421            pass
     2422
     2423        cv1 = domain1.quantities['stage'].centroid_values
     2424       
     2425        #print take(cv1, (12,14,16))  #Right
     2426        #print take(cv1, (0,6,12))  #Bottom
     2427       
     2428       
     2429        #Create an triangle shaped domain (reusing coordinates from domain 1),
     2430        #formed from the lower and right hand  boundaries and
     2431        #the sw-ne diagonal
     2432        #from domain 1. Call it domain2
     2433       
     2434        points = [ [0,0], [1.0/3,0], [1.0/3,1.0/3],
     2435                   [2.0/3,0], [2.0/3,1.0/3], [2.0/3,2.0/3],
     2436                   [1,0], [1,1.0/3], [1,2.0/3], [1,1]]
     2437       
     2438        vertices = [ [1,2,0],
     2439                     [3,4,1], [2,1,4], [4,5,2],
     2440                     [6,7,3], [4,3,7], [7,8,4], [5,4,8], [8,9,5]]       
     2441
     2442        boundary = { (0,1):'bottom', (1,1):'bottom', (4,1): 'bottom',
     2443                     (4,2):'right', (6,2):'right', (8,2):'right',
     2444                     (0,0):'diagonal', (3,0):'diagonal', (8,0):'diagonal'}   
     2445                     
     2446        domain2 = Domain(points, vertices, boundary)
     2447       
     2448        domain2.reduction = domain1.reduction
     2449        domain2.smooth = False
     2450        domain2.visualise = True
     2451        domain2.default_order = 2
     2452       
     2453        #Bed-slope and friction at vertices (and interpolated elsewhere)
     2454        domain2.set_quantity('elevation', 0)
     2455        domain2.set_quantity('friction', 0)     
     2456
     2457        # Boundary conditions
     2458        Br = Reflective_boundary(domain2)
     2459        #Bd = Dirichlet_boundary([0.3,0,0])     
     2460        Bt = Transmissive_boundary(domain2)
     2461       
     2462        return
     2463       
     2464        Bf = Spatio_temporal_boundary(domain1.filename + '.' + domain1.format,
     2465                                      domain2)
     2466        domain2.set_boundary({'right':Bt, 'bottom':Br, 'diagonal':Bf})
     2467
     2468        #Initial condition
     2469        domain2.set_quantity('stage', 0)
     2470        domain2.check_integrity()       
     2471
     2472        #Evolution
     2473        for t in domain2.evolve(yieldstep = 0.05, finaltime = 4):
     2474            pass
     2475
    24272476       
    24282477        #Use output from domain1 as spatio-temporal boundary for domain2
    24292478        #and verify that results at right hand side are close. 
    24302479
    2431        
    2432        
     2480        cv2 = domain2.quantities['stage'].centroid_values
     2481       
     2482        print take(cv1, (12,14,16))  #Right
     2483        print take(cv2, (4,6,8))
     2484       
     2485        #print take(cv1, (0,6,12))  #Bottom     
     2486       
     2487        assert allclose( take(cv1, (12,14,16)), take(cv2, (4,6,8))) #RHS
     2488        assert allclose( take(cv1, (0,6,12)), take(cv2, (0,1,4))) #Bottom       
     2489       
     2490        #Cleanup
     2491        os.remove(domain1.filename + '.' + domain1.format)
    24332492   
    24342493#-------------------------------------------------------------
    24352494if __name__ == "__main__":
    2436     suite = unittest.makeSuite(TestCase,'test')
    2437     #suite = unittest.makeSuite(TestCase,'test_boundary_conditionsII')
     2495    #suite = unittest.makeSuite(TestCase,'test')
     2496    suite = unittest.makeSuite(TestCase,'test_spatio_temporal_boundary')
    24382497    runner = unittest.TextTestRunner()
    24392498    runner.run(suite)

  • inundation/ga/storm_surge/pyvolution/test_util.py

    r841 r849  
    166166        fid.close()
    167167
    168         F = File_function(filename)
     168        F = file_function(filename)
    169169
    170170        #Now try interpolation
    … …  
    229229
    230230        #Check that domain.starttime is updated if non-existing
    231         F = File_function(filename, domain)
     231        F = file_function(filename, domain)
    232232        assert allclose(domain.starttime, start)       
    233233
    234234        #Check that domain.starttime is updated if too early
    235235        domain.starttime = start - 1               
    236         F = File_function(filename, domain)
     236        F = file_function(filename, domain)
    237237        assert allclose(domain.starttime, start)
    238238
    239239        #Check that domain.starttime isn't updated if later
    240240        #domain.starttime = start + 1               
    241         #F = File_function(filename, domain)
     241        #F = file_function(filename, domain)
    242242        #assert allclose(domain.starttime, start+1)               
    243243
    … …  
    307307        delta = 23
    308308        domain.starttime = start + delta               
    309         F = File_function(filename, domain)
     309        F = file_function(filename, domain)
    310310        assert allclose(domain.starttime, start+delta)               
    311311

  • inundation/ga/storm_surge/pyvolution/util.py

    r844 r849  
    8080      return False 
    8181   
    82              
    83 class File_function:
     82
     83def file_function(filename, domain=None): 
     84    assert type(filename) == type(''),\
     85               'First argument to File_function must be a string'
     86
     87    try:
     88        fid = open(filename)
     89    except Exception, e:
     90        msg = 'File "%s" could not be opened: Error="%s"'\
     91                  %(filename, e)
     92        raise msg
     93
     94    line = fid.readline()
     95    fid.close()
     96
     97    #Choose format
     98    #FIXME: Maybe these can be merged later on
     99    if line[:3] == 'CDF':
     100        return File_function_NetCDF(filename, domain)
     101    else:
     102        return File_function_ASCII(filename, domain)   
     103           
     104             
     105class File_function_NetCDF:
     106    """Read time history of spatial data from NetCDF sww file and
     107    return a callable object f(t,x,y) 
     108    which will return interpolated values based on the input file.   
     109   
     110    x, y may be either scalars or vectors
     111   
     112    #FIXME: More about format, interpolation  and order of quantities
     113    """
     114
     115    def __init__(self, filename, domain=None):
     116        """Initialise callable object from a file.
     117        See docstring for class File_function_NetCDF
     118
     119        If domain is specified, model time (domain.starttime)
     120        will be checked and possibly modified
     121
     122        All times are assumed to be in UTC
     123        """
     124
     125        import time, calendar
     126        from Numeric import array
     127        from config import time_format
     128        from Scientific.IO.NetCDF import NetCDFFile     
     129
     130        #Open NetCDF file
     131        fid = NetCDFFile(filename, 'r')
     132   
     133        print fid.variables.keys()
     134   
     135        fields = line.split(',')
     136        msg = 'File %s must have the format date, value0 value1 value2 ...'
     137        msg += ' or date, x, y, value0 value1 value2 ...'
     138        mode = len(fields)
     139        assert mode in [2,4], msg
     140
     141        try:
     142            starttime = calendar.timegm(time.strptime(fields[0], time_format))
     143        except ValueError:   
     144            msg = 'First field in file %s must be' %filename
     145            msg += ' date-time with format %s.\n' %time_format
     146            msg += 'I got %s instead.' %fields[0]
     147            raise msg
     148
     149
     150        #Split values
     151        values = []
     152        for value in fields[mode-1].split():
     153            values.append(float(value))
     154
     155        q = array(values)
     156       
     157        msg = 'ERROR: File must contain at least one independent value'
     158        assert len(q.shape) == 1, msg
     159           
     160        self.number_of_values = len(q)
     161
     162        self.filename = filename
     163        self.starttime = starttime
     164        self.domain = domain
     165
     166        if domain is not None:
     167            if domain.starttime is None:
     168                domain.starttime = self.starttime
     169            else:
     170                msg = 'WARNING: Start time as specified in domain (%s)'\
     171                      %domain.starttime
     172                msg += ' is earlier than the starttime of file %s: %s.'\
     173                         %(self.filename, self.starttime)
     174                msg += 'Modifying starttime accordingly.'
     175                if self.starttime > domain.starttime:
     176                    #FIXME: Print depending on some verbosity setting
     177                    #print msg
     178                    domain.starttime = self.starttime #Modifying model time
     179
     180        if mode == 2:       
     181            self.read_times() #Now read all times in.
     182        else:
     183            raise 'x,y dependency not yet implemented'
     184
     185   
     186             
     187class File_function_ASCII:
     188    """Read time series from file and return a callable object:
     189    f(t,x,y)
     190    which will return interpolated values based on the input file.
     191
     192    The file format is assumed to be either two fields separated by a comma:
     193
     194        time [DD/MM/YY hh:mm:ss], value0 value1 value2 ...
     195
     196    or four comma separated fields
     197
     198        time [DD/MM/YY hh:mm:ss], x, y, value0 value1 value2 ...
     199
     200    In either case, the callable object will return a tuple of
     201    interpolated values, one each value stated in the file.
     202
     203   
     204    E.g
     205
     206      31/08/04 00:00:00, 1.328223 0 0
     207      31/08/04 00:15:00, 1.292912 0 0
     208
     209    will provide a time dependent function f(t,x=None,y=None),
     210    ignoring x and y, which returns three values per call.
     211
     212   
     213    NOTE: At this stage the function is assumed to depend on
     214    time only, i.e  no spatial dependency!!!!!
     215    When that is needed we can use the least_squares interpolation.
     216   
     217    #FIXME: This should work with netcdf (e.g. sww) and thus render the
     218    #spatio-temporal boundary condition in shallow water fully general
     219    """
     220
     221   
     222    def __init__(self, filename, domain=None):
     223        """Initialise callable object from a file.
     224        See docstring for class File_function
     225
     226        If domain is specified, model time (domain,starttime)
     227        will be checked and possibly modified
     228
     229        All times are assumed to be in UTC
     230        """
     231
     232        import time, calendar
     233        from Numeric import array
     234        from config import time_format
     235
     236        fid = open(filename)
     237        line = fid.readline()
     238        fid.close()
     239   
     240        fields = line.split(',')
     241        msg = 'File %s must have the format date, value0 value1 value2 ...'
     242        msg += ' or date, x, y, value0 value1 value2 ...'
     243        mode = len(fields)
     244        assert mode in [2,4], msg
     245
     246        try:
     247            starttime = calendar.timegm(time.strptime(fields[0], time_format))
     248        except ValueError:   
     249            msg = 'First field in file %s must be' %filename
     250            msg += ' date-time with format %s.\n' %time_format
     251            msg += 'I got %s instead.' %fields[0]
     252            raise msg
     253
     254
     255        #Split values
     256        values = []
     257        for value in fields[mode-1].split():
     258            values.append(float(value))
     259
     260        q = array(values)
     261       
     262        msg = 'ERROR: File must contain at least one independent value'
     263        assert len(q.shape) == 1, msg
     264           
     265        self.number_of_values = len(q)
     266
     267        self.filename = filename
     268        self.starttime = starttime
     269        self.domain = domain
     270
     271        if domain is not None:
     272            if domain.starttime is None:
     273                domain.starttime = self.starttime
     274            else:
     275                msg = 'WARNING: Start time as specified in domain (%s)'\
     276                      %domain.starttime
     277                msg += ' is earlier than the starttime of file %s: %s.'\
     278                         %(self.filename, self.starttime)
     279                msg += 'Modifying starttime accordingly.'
     280                if self.starttime > domain.starttime:
     281                    #FIXME: Print depending on some verbosity setting
     282                    #print msg
     283                    domain.starttime = self.starttime #Modifying model time
     284
     285        if mode == 2:       
     286            self.read_times() #Now read all times in.
     287        else:
     288            raise 'x,y dependency not yet implemented'
     289
     290       
     291    def read_times(self):
     292        """Read time and values
     293        """
     294        from Numeric import zeros, Float, alltrue
     295        from config import time_format
     296        import time, calendar
     297       
     298        fid = open(self.filename)       
     299        lines = fid.readlines()
     300        fid.close()
     301       
     302        N = len(lines)
     303        d = self.number_of_values
     304
     305        T = zeros(N, Float)       #Time
     306        Q = zeros((N, d), Float)  #Values
     307       
     308        for i, line in enumerate(lines):
     309            fields = line.split(',')
     310            realtime = calendar.timegm(time.strptime(fields[0], time_format))
     311
     312            T[i] = realtime - self.starttime
     313
     314            for j, value in enumerate(fields[1].split()):
     315                Q[i, j] = float(value)
     316
     317        msg = 'File %s must list time as a monotonuosly ' %self.filename
     318        msg += 'increasing sequence'
     319        assert alltrue( T[1:] - T[:-1] > 0 ), msg
     320
     321        self.T = T     #Time
     322        self.Q = Q     #Values
     323        self.index = 0 #Initial index
     324
     325       
     326    def __repr__(self):
     327        return 'File function'
     328
     329       
     330
     331    def __call__(self, t, x=None, y=None):
     332        """Evaluate f(t,x,y)
     333
     334        FIXME: x, y dependency not yet implemented except that
     335        result is a vector of same length as x and y
     336        FIXME: Naaaa
     337        """
     338
     339        from math import pi, cos, sin, sqrt
     340       
     341
     342        #Find time tau such that
     343        #
     344        # domain.starttime + t == self.starttime + tau
     345       
     346        if self.domain is not None:
     347            tau = self.domain.starttime-self.starttime+t
     348        else:
     349            tau = t
     350       
     351               
     352        msg = 'Time interval derived from file %s (%s:%s) does not match model time: %s'\
     353              %(self.filename, self.T[0], self.T[1], tau)
     354        if tau < self.T[0]: raise msg
     355        if tau > self.T[-1]: raise msg       
     356
     357        oldindex = self.index
     358
     359        #Find slot
     360        while tau > self.T[self.index]: self.index += 1
     361        while tau < self.T[self.index]: self.index -= 1
     362
     363        #t is now between index and index+1
     364        ratio = (tau - self.T[self.index])/\
     365                (self.T[self.index+1] - self.T[self.index])
     366
     367        #Compute interpolated values
     368        q = self.Q[self.index,:] +\
     369            ratio*(self.Q[self.index+1,:] - self.Q[self.index,:])
     370
     371        #Return vector of interpolated values
     372        if x == None and y == None:
     373            return q
     374        else:
     375            try:
     376                N = len(x)
     377            except:
     378                return q
     379            else:
     380                from Numeric import ones, Float
     381                #x is a vector - Create one constant column for each value
     382                N = len(x)
     383                assert len(y) == N, 'x and y must have same length'
     384             
     385                res = []
     386                for col in q:
     387                    res.append(col*ones(N, Float))
     388           
     389                return res
     390           
     391           
     392#FIXME: TEMP       
     393class File_function_copy:
    84394    """Read time series from file and return a callable object:
    85395    f(t,x,y)
Note: See TracChangeset for help on using the changeset viewer.