Changeset 5582


Ignore:
Timestamp:
Jul 29, 2008, 10:51:48 AM (11 years ago)
Author:
ole
Message:

Comments and more testing trying to find issues with urs2sts.

Location:
anuga_core/source/anuga/shallow_water
Files:
3 edited

Legend:

Unmodified
Added
Removed
  • anuga_core/source/anuga/shallow_water/data_manager.py

    r5579 r5582  
    50115011        msg = 'When combining multiple sources must specify a weight for '+\
    50125012              'mux2 source file'
    5013         assert len(weights)== numSrc, msg
     5013        assert len(weights) == numSrc, msg
     5014       
     5015    if verbose:
     5016        print 'Weights used in urs2sts:', weights
    50145017
    50155018    # Check output filename   
  • anuga_core/source/anuga/shallow_water/test_data_manager.py

    r5556 r5582  
    61316131                                      va=va)
    61326132
    6133         weights=ones(1,Float)
     6133        weights=ones(1, Float)
    61346134        #ensure that files are indeed mux2 files
    61356135        times, latitudes, longitudes, elevation, stage, starttime =read_mux2_py([files[0]],weights)
     
    61886188                                      va=va)
    61896189
    6190         weights=ones(1,Float)
     6190        weights=ones(1, Float)
    61916191        #ensure that files are indeed mux2 files
    61926192        times, latitudes, longitudes, elevation, stage,starttime=read_mux2_py([files[0]],weights)
     
    62496249                                      va=va)
    62506250
    6251         weights=ones(1,Float)
     6251        weights=ones(1, Float)
    62526252        #ensure that files are indeed mux2 files
    62536253        times, latitudes, longitudes, elevation, stage,starttime=read_mux2_py([files[0]],weights)
     
    64486448        urs2sts([base_nameI, base_nameII],
    64496449                basename_out=base_nameI,
    6450                 weights=[1.0,1.0],
     6450                weights=[1.0, 1.0],
    64516451                mean_stage=tide,
    64526452                verbose=False)
     
    66046604                basename_out=base_nameI,
    66056605                ordering_filename=ordering_filename,
    6606                 weights=[1.0,1.0],
     6606                weights=[1.0, 1.0],
    66076607                mean_stage=tide,
    66086608                verbose=False)
     
    67826782                    basename_out=base_nameI,
    67836783                    ordering_filename=ordering_filename,
    6784                     weights=[1.0,1.0],
     6784                    weights=[1.0, 1.0],
    67856785                    mean_stage=tide,
    67866786                    verbose=False) 
     
    68146814       
    68156815        # Create non-trivial weights
    6816         weights = [0.8, 1.5]
     6816        #weights = [0.8, 1.5] # OK
     6817        #weights = [0.8, 10.5] # Fail (up to allclose tolerance)
     6818        #weights = [10.5, 10.5] # OK
     6819        #weights = [0.0, 10.5] # OK
     6820        #weights = [0.8, 0.] # OK               
     6821        #weights = [8, 0.1] # OK                       
     6822        #weights = [0.8, 10.0] # OK                               
     6823        #weights = [0.8, 10.6] # OK           
     6824        weights = [3.8, 7.6] # OK                   
     6825        #weights = [0.5, 0.5] # OK                           
     6826        #weights = [2., 2.] # OK                           
     6827        #weights = [0.0, 0.5] # OK                                         
     6828        #weights = [1.0, 1.0] # OK                                                 
     6829       
    68176830       
    68186831        # Create different timeseries starting and ending at different times
     
    68606873
    68616874
     6875                 
    68626876        #print
    68636877        #print 'using varying start and end time'
    6864         #print 'ha0', ha0
     6878        #print 'ha0', ha0[4]
    68656879        #print 'ua0', ua0
    68666880        #print 'va0', va0       
     
    68786892                                             
    68796893        # Create data to be written to second mux file       
    6880         ha1=zeros((n,time_step_count),Float)
     6894        ha1=ones((n,time_step_count),Float)
    68816895        ha1[0]=sin(times_ref)
    68826896        ha1[1]=2*sin(times_ref - 3)
     
    69096923        #print
    69106924        #print 'using varying start and end time'
    6911         #print 'ha1', ha1
     6925        #print 'ha1', ha1[4]
    69126926        #print 'ua1', ua1
    69136927        #print 'va1', va1       
     
    69386952           
    69396953
    6940                                                
    6941         # Call urs2sts with multiple mux files
    6942         urs2sts([base_nameI, base_nameII],
     6954        #------------------------------------------------------------
     6955        # Now read the mux files one by one with out weights and test
     6956       
     6957        # Call urs2sts with mux file #0
     6958        urs2sts([base_nameI],
    69436959                basename_out=base_nameI,
    69446960                ordering_filename=ordering_filename,
    6945                 weights=weights,
    69466961                mean_stage=tide,
    69476962                verbose=False)
     
    69796994                       
    69806995
     6996        # Check sts values for mux #0
     6997        stage0 = fid.variables['stage'][:].copy()
     6998        xmomentum0 = fid.variables['xmomentum'][:].copy()
     6999        ymomentum0 = fid.variables['ymomentum'][:].copy()
     7000        elevation0 = fid.variables['elevation'][:].copy()
     7001
     7002       
     7003        #print 'stage', stage0
     7004        #print 'xmomentum', xmomentum0
     7005        #print 'ymomentum', ymomentum0       
     7006        #print 'elevation', elevation0
     7007       
     7008        # The quantities stored in the .sts file should be the weighted sum of the
     7009        # quantities written to the mux2 files subject to the permutation vector.
     7010       
     7011        ha_ref = take(ha0, permutation)
     7012        ua_ref = take(ua0, permutation)       
     7013        va_ref = take(va0, permutation)               
     7014
     7015        gauge_depth_ref = take(gauge_depth, permutation)                         
     7016
     7017
     7018        #print
     7019        #print stage0
     7020        #print transpose(ha_ref)+tide - stage0
     7021
     7022       
     7023        assert allclose(transpose(ha_ref)+tide, stage0)  # Meters
     7024       
     7025        #Check the momentums - ua
     7026        #momentum = velocity*(stage-elevation)
     7027        # elevation = - depth
     7028        #momentum = velocity_ua *(stage+depth)
     7029
     7030        depth_ref = zeros((len(permutation), time_step_count), Float)
     7031        for i in range(len(permutation)):
     7032            depth_ref[i]=gauge_depth_ref[i]+tide+ha_ref[i]
     7033
     7034
     7035        # The xmomentum stored in the .sts file should be the sum of the ua
     7036        # in the two mux2 files multiplied by the depth.
     7037        assert allclose(transpose(ua_ref*depth_ref), xmomentum0)
     7038
     7039        #Check the momentums - va
     7040        #momentum = velocity*(stage-elevation)
     7041        # elevation = - depth
     7042        #momentum = velocity_va *(stage+depth)
     7043
     7044        # The ymomentum stored in the .sts file should be the sum of the va
     7045        # in the two mux2 files multiplied by the depth.
     7046       
     7047       
     7048        #print transpose(va_ref*depth_ref)
     7049        #print ymomentum
     7050        assert allclose(transpose(va_ref*depth_ref), ymomentum0)       
     7051
     7052        # check the elevation values.
     7053        # -ve since urs measures depth, sww meshers height,
     7054        assert allclose(-gauge_depth_ref, elevation0)
     7055
     7056        fid.close()
     7057        os.remove(sts_file)
     7058       
     7059       
     7060
     7061       
     7062        # Call urs2sts with mux file #1
     7063        urs2sts([base_nameII],
     7064                basename_out=base_nameI,
     7065                ordering_filename=ordering_filename,
     7066                mean_stage=tide,
     7067                verbose=False)
     7068
     7069        # Now read the sts file and check that values have been stored correctly.
     7070        sts_file = base_nameI + '.sts'
     7071        fid = NetCDFFile(sts_file)
     7072
     7073        # Make x and y absolute
     7074        x = fid.variables['x'][:]
     7075        y = fid.variables['y'][:]
     7076
     7077        geo_reference = Geo_reference(NetCDFObject=fid)
     7078        points = geo_reference.get_absolute(map(None, x, y))
     7079        points = ensure_numeric(points)
     7080
     7081        x = points[:,0]
     7082        y = points[:,1]
     7083
     7084        for i, index in enumerate(permutation):
     7085            # Check that STS points are stored in the correct order
     7086           
     7087            # Work out the UTM coordinates sts point i
     7088            zone, e, n = redfearn(lat_long_points[index][0],
     7089                                  lat_long_points[index][1])             
     7090
     7091            #print i, [x[i],y[i]], [e,n]
     7092            assert allclose([x[i],y[i]], [e,n])
     7093           
     7094                       
     7095        # Check the time vector
     7096        times = fid.variables['time'][:]
     7097        assert allclose(ensure_numeric(times),
     7098                        ensure_numeric(times_ref))
     7099                       
     7100
     7101        # Check sts values for mux #1
     7102        stage1 = fid.variables['stage'][:].copy()
     7103        xmomentum1 = fid.variables['xmomentum'][:].copy()
     7104        ymomentum1 = fid.variables['ymomentum'][:].copy()
     7105        elevation1 = fid.variables['elevation'][:].copy()
     7106
     7107       
     7108        #print 'stage', stage1
     7109        #print 'xmomentum', xmomentum1
     7110        #print 'ymomentum', ymomentum1       
     7111        #print 'elevation', elevation1
     7112       
     7113        # The quantities stored in the .sts file should be the weighted sum of the
     7114        # quantities written to the mux2 files subject to the permutation vector.
     7115       
     7116        ha_ref = take(ha1, permutation)
     7117        ua_ref = take(ua1, permutation)       
     7118        va_ref = take(va1, permutation)               
     7119
     7120        gauge_depth_ref = take(gauge_depth, permutation)                         
     7121
     7122
     7123        #print
     7124        #print stage1
     7125        #print transpose(ha_ref)+tide - stage1
     7126       
     7127
     7128        assert allclose(transpose(ha_ref)+tide, stage1)  # Meters
     7129        #import sys; sys.exit()
     7130
     7131        #Check the momentums - ua
     7132        #momentum = velocity*(stage-elevation)
     7133        # elevation = - depth
     7134        #momentum = velocity_ua *(stage+depth)
     7135
     7136        depth_ref = zeros((len(permutation), time_step_count), Float)
     7137        for i in range(len(permutation)):
     7138            depth_ref[i]=gauge_depth_ref[i]+tide+ha_ref[i]
     7139
     7140
     7141        # The xmomentum stored in the .sts file should be the sum of the ua
     7142        # in the two mux2 files multiplied by the depth.
     7143        assert allclose(transpose(ua_ref*depth_ref), xmomentum1)
     7144
     7145        #Check the momentums - va
     7146        #momentum = velocity*(stage-elevation)
     7147        # elevation = - depth
     7148        #momentum = velocity_va *(stage+depth)
     7149
     7150        # The ymomentum stored in the .sts file should be the sum of the va
     7151        # in the two mux2 files multiplied by the depth.
     7152       
     7153       
     7154        #print transpose(va_ref*depth_ref)
     7155        #print ymomentum
     7156        assert allclose(transpose(va_ref*depth_ref), ymomentum1)       
     7157
     7158        # check the elevation values.
     7159        # -ve since urs measures depth, sww meshers height,
     7160        assert allclose(-gauge_depth_ref, elevation1)
     7161
     7162        fid.close()
     7163        os.remove(sts_file)
     7164       
     7165       
     7166       
     7167        #----------------------
     7168        # Then read the mux files together and test
     7169       
     7170                                               
     7171        # Call urs2sts with multiple mux files
     7172        urs2sts([base_nameI, base_nameII],
     7173                basename_out=base_nameI,
     7174                ordering_filename=ordering_filename,
     7175                weights=weights,
     7176                mean_stage=tide,
     7177                verbose=False)
     7178
     7179        # Now read the sts file and check that values have been stored correctly.
     7180        sts_file = base_nameI + '.sts'
     7181        fid = NetCDFFile(sts_file)
     7182
     7183        # Make x and y absolute
     7184        x = fid.variables['x'][:]
     7185        y = fid.variables['y'][:]
     7186
     7187        geo_reference = Geo_reference(NetCDFObject=fid)
     7188        points = geo_reference.get_absolute(map(None, x, y))
     7189        points = ensure_numeric(points)
     7190
     7191        x = points[:,0]
     7192        y = points[:,1]
     7193
     7194        for i, index in enumerate(permutation):
     7195            # Check that STS points are stored in the correct order
     7196           
     7197            # Work out the UTM coordinates sts point i
     7198            zone, e, n = redfearn(lat_long_points[index][0],
     7199                                  lat_long_points[index][1])             
     7200
     7201            #print i, [x[i],y[i]], [e,n]
     7202            assert allclose([x[i],y[i]], [e,n])
     7203           
     7204                       
     7205        # Check the time vector
     7206        times = fid.variables['time'][:]
     7207        assert allclose(ensure_numeric(times),
     7208                        ensure_numeric(times_ref))
     7209                       
     7210
    69817211        # Check sts values
    6982         stage = fid.variables['stage'][:]
    6983         xmomentum = fid.variables['xmomentum'][:]
    6984         ymomentum = fid.variables['ymomentum'][:]
    6985         elevation = fid.variables['elevation'][:]
     7212        stage = fid.variables['stage'][:].copy()
     7213        xmomentum = fid.variables['xmomentum'][:].copy()
     7214        ymomentum = fid.variables['ymomentum'][:].copy()
     7215        elevation = fid.variables['elevation'][:].copy()
    69867216
    69877217       
     
    69987228        gauge_depth_ref = take(gauge_depth, permutation)                         
    69997229
    7000        
     7230
     7231        #print
     7232        #print stage
     7233        #print transpose(ha_ref)+tide - stage
     7234
    70017235        assert allclose(transpose(ha_ref)+tide, stage)  # Meters
    70027236
     
    70477281       
    70487282
     7283       
     7284        #---------------
     7285        # "Manually" add the timeseries up with weights and test
     7286        # Tide is discounted from individual results and added back in       
     7287        #
     7288
     7289        stage_man = weights[0]*(stage0-tide) + weights[1]*(stage1-tide) + tide
     7290        assert allclose(stage_man, stage)
     7291       
    70497292       
    70507293       
  • anuga_core/source/anuga/shallow_water/urs_ext.c

    r5560 r5582  
    125125char isdata(float x)
    126126{
    127     //char value;
    128127    if(x < NODATA + EPSILON && NODATA < x + EPSILON)
    129128    {
    130        return 0;
     129      return 0;
    131130    }
    132131    else
    133132    {
    134         return 1; 
     133      return 1; 
    135134    }
    136135}
     
    272271
    273272        /* Compute the size of the data block for this source */
    274         numData = getNumData(fros + i*(*total_number_of_stations), lros + i*(*total_number_of_stations), (*total_number_of_stations));
     273        numData = getNumData(fros + i*(*total_number_of_stations),
     274                             lros + i*(*total_number_of_stations),
     275                             (*total_number_of_stations));
    275276
    276277        /* Sanity check */
     
    289290    }
    290291
     292    // Store time resolution and number of timesteps   
     293    // These are the same for all stations, so
     294    // we take the first one.
    291295    *delta_t = (double)mytgs0[0].dt;
    292296    *number_of_time_steps = mytgs0[0].nt;
    293297
    294         free(mytgs);
     298    free(mytgs);
    295299
    296300    return 1;
Note: See TracChangeset for help on using the changeset viewer.