 Timestamp:
 Apr 7, 2009, 4:00:51 PM (14 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

anuga_core/source/anuga/shallow_water/test_data_manager.py
r6711 r6752 6042 6042 quantities_init[2].append(va[i]) # South is negative in MUX 6043 6043 6044 file_handle, base_name = tempfile.mkstemp(" ")6044 file_handle, base_name = tempfile.mkstemp("write_mux2") 6045 6045 os.close(file_handle) 6046 6046 os.remove(base_name) … … 6103 6103 #print 'writing', time, point_i, q_time[time, point_i] 6104 6104 f.write(pack('f', q_time[time, point_i])) 6105 6106 6105 f.close() 6107 6106 … … 6163 6162 msg='incorrect gauge va time series returned' 6164 6163 assert num.allclose(yvelocity, va) 6164 6165 6165 6166 6166 6167 def test_urs2sts_read_mux2_pyII(self): … … 6396 6397 if j == 2: assert num.allclose(data[i][:parameters_index], va0[permutation[i], :]) 6397 6398 6399 self.delete_mux(filesI) 6398 6400 6399 6401 … … 6406 6408 wrong values Win32 6407 6409 6408 This test does not pass on Windows but test_read_mux_platform_problem1 does 6410 This test does not pass on Windows but test_read_mux_platform_problem1 6411 does 6409 6412 """ 6410 6413 … … 6421 6424 times_ref = num.arange(0, time_step_count*time_step, time_step) 6422 6425 6423 lat_long_points = [(21.5,114.5), (21,114.5), (21.5,115), (21.,115.), (22., 117.)] 6426 lat_long_points = [(21.5,114.5), (21,114.5), (21.5,115), 6427 (21.,115.), (22., 117.)] 6424 6428 n = len(lat_long_points) 6425 6429 … … 6458 6462 va1[1]=3*num.ones(time_step_count) 6459 6463 va1[3]=2*num.sin(times_ref0.71) 6460 6464 #print "va1[0]", va1[0] # The 8th element is what will go bad. 6461 6465 # Ensure data used to write mux file to be zero when gauges are 6462 6466 # not recording 6463 6467 for i in range(n): 6464 6468 # For each point 6465 for j in range(0, first_tstep[i]1) + range(last_tstep[i], time_step_count): 6469 for j in range(0, first_tstep[i]1) + range(last_tstep[i], 6470 time_step_count): 6466 6471 # For timesteps before and after recording range 6467 ha1[i][j] = ua1[i][j] = va1[i][j] = 0.0 6472 ha1[i][j] = ua1[i][j] = va1[i][j] = 0.0 6468 6473 6469 6474 … … 6473 6478 #print 'va', va1[0,:] 6474 6479 6475 # Write second mux file to be combined by urs2sts 6480 # Write second mux file to be combined by urs2sts 6476 6481 base_nameII, filesII = self.write_mux2(lat_long_points, 6477 6482 time_step_count, time_step, … … 6527 6532 if ha is None: 6528 6533 this_ha = e 6529 quantities_init[0].append(num.ones(time_step_count,num.Float)*this_ha) # HA 6534 quantities_init[0].append(num.ones(time_step_count, 6535 num.Float)*this_ha) # HA 6530 6536 else: 6531 6537 quantities_init[0].append(ha[i]) 6532 6538 if ua is None: 6533 6539 this_ua = n 6534 quantities_init[1].append(num.ones(time_step_count,num.Float)*this_ua) # UA 6540 quantities_init[1].append(num.ones(time_step_count, 6541 num.Float)*this_ua) # UA 6535 6542 else: 6536 6543 quantities_init[1].append(ua[i]) 6537 6544 if va is None: 6538 6545 this_va = e 6539 quantities_init[2].append(num.ones(time_step_count,num.Float)*this_va) # 6546 quantities_init[2].append(num.ones(time_step_count, 6547 num.Float)*this_va) # 6540 6548 else: 6541 6549 quantities_init[2].append(va[i]) … … 6601 6609 #print time, x, q_time[time, point_i] 6602 6610 if q == 'VA': x = x # South is positive in MUX 6611 #if q == 'VA': 6612 #print q+" q_time[%d, %d] = %f" %(time, point_i, 6613 # q_time[time, point_i]) 6614 # q_time[time, point_i] 6603 6615 assert abs(q_time[time, point_i]x)<epsilon 6604 6616 6605 6617 f.close() 6606 6607 6608 6609 6610 6611 6612 6613 6618 6614 6619 # Create ordering file 6615 6620 permutation = ensure_numeric([4,0,2]) 6616 6621 6617 _, ordering_filename = tempfile.mkstemp('') 6618 order_fid = open(ordering_filename, 'w') 6619 order_fid.write('index, longitude, latitude\n') 6620 for index in permutation: 6621 order_fid.write('%d, %f, %f\n' %(index, 6622 lat_long_points[index][1], 6623 lat_long_points[index][0])) 6624 order_fid.close() 6625 6626 6627 6622 # _, ordering_filename = tempfile.mkstemp('') 6623 # order_fid = open(ordering_filename, 'w') 6624 # order_fid.write('index, longitude, latitude\n') 6625 # for index in permutation: 6626 # order_fid.write('%d, %f, %f\n' %(index, 6627 # lat_long_points[index][1], 6628 # lat_long_points[index][0])) 6629 # order_fid.close() 6630 6628 6631 #  6629 6632 # Now read files back and check values … … 6637 6640 for j, file in enumerate(filesII): 6638 6641 # Read stage, u, v enumerated as j 6639 6640 6641 6642 #print 'Reading', j, file 6642 data = read_mux2(1, [file], weights, file_params, permutation, verbose) 6643 data = read_mux2(1, [file], weights, file_params, 6644 permutation, verbose) 6643 6645 6644 6646 #print 'Data received by Python' 6645 6647 #print data[1][8] 6646 6647 6648 6648 number_of_selected_stations = data.shape[0] 6649 6649 … … 6651 6651 parameters_index = data.shape[1]OFFSET 6652 6652 6653 quantity=num.zeros((number_of_selected_stations, parameters_index), num.Float) 6653 quantity=num.zeros((number_of_selected_stations, 6654 parameters_index), num.Float) 6654 6655 6655 6656 … … 6658 6659 #print i, parameters_index 6659 6660 #print quantity[i][:] 6660 6661 6662 6661 if j == 0: assert num.allclose(data[i][:parameters_index], ha1[permutation[i], :]) 6663 6662 if j == 1: assert num.allclose(data[i][:parameters_index], ua1[permutation[i], :]) … … 6668 6667 #print j, i 6669 6668 #print 'Input' 6670 #print 'u', ua1[permutation[i], 8] 6669 #print 'u', ua1[permutation[i], 8] 6671 6670 #print 'v', va1[permutation[i], 8] 6672 6671 6673 6672 #print 'Output' 6674 #print 'v ', data[i][:parameters_index][8] 6673 #print 'v ', data[i][:parameters_index][8] 6674 6675 # South is positive in MUX 6676 #print "data[i][:parameters_index]", data[i][:parameters_index] 6677 #print "va1[permutation[i], :]", va1[permutation[i], :] 6678 assert num.allclose(data[i][:parameters_index], va1[permutation[i], :]) 6679 6680 self.delete_mux(filesII) 6681 6682 def test_read_mux_platform_problem3(self): 6683 6684 # This is to test a situation where read_mux returned 6685 # wrong values Win32 6686 6687 6688 from urs_ext import read_mux2 6689 6690 from anuga.config import single_precision as epsilon 6691 6692 verbose = False 6675 6693 6676 # South is positive in MUX 6694 tide = 1.5 6695 time_step_count = 10 6696 time_step = 0.02 6697 6698 ''' 6699 Win results 6700 time_step = 0.2000001 6701 This is OK 6702 ''' 6703 6704 ''' 6705 Win results 6706 time_step = 0.20000001 6707 6708 ====================================================================== 6709 ERROR: test_read_mux_platform_problem3 (__main__.Test_Data_Manager) 6710  6711 Traceback (most recent call last): 6712 File "test_data_manager.py", line 6718, in test_read_mux_platform_problem3 6713 ha1[0]=num.sin(times_ref) 6714 ValueError: matrices are not aligned for copy 6715 6716 ''' 6717 6718 ''' 6719 Win results 6720 time_step = 0.200000001 6721 FAIL 6722 assert num.allclose(data[i][:parameters_index], 6723 va1[permutation[i], :]) 6724 ''' 6725 times_ref = num.arange(0, time_step_count*time_step, time_step) 6726 #print "times_ref", times_ref 6727 6728 lat_long_points = [(21.5,114.5), (21,114.5), (21.5,115), 6729 (21.,115.), (22., 117.)] 6730 stations = len(lat_long_points) 6731 6732 # Create different timeseries starting and ending at different times 6733 first_tstep=num.ones(stations,num.Int) 6734 first_tstep[0]+=2 # Point 0 starts at 2 6735 first_tstep[1]+=4 # Point 1 starts at 4 6736 first_tstep[2]+=3 # Point 2 starts at 3 6737 6738 last_tstep=(time_step_count)*num.ones(stations,num.Int) 6739 last_tstep[0]=1 # Point 0 ends 1 step early 6740 last_tstep[1]=2 # Point 1 ends 2 steps early 6741 last_tstep[4]=3 # Point 4 ends 3 steps early 6742 6743 # Create varying elevation data (positive values for seafloor) 6744 gauge_depth=20*num.ones(stations,num.Float) 6745 for i in range(stations): 6746 gauge_depth[i] += i**2 6747 6748 # Create data to be written to second mux file 6749 ha1=num.ones((stations,time_step_count),num.Float) 6750 ha1[0]=num.sin(times_ref) 6751 ha1[1]=2*num.sin(times_ref  3) 6752 ha1[2]=5*num.sin(4*times_ref) 6753 ha1[3]=num.sin(times_ref) 6754 ha1[4]=num.sin(2*times_ref0.7) 6755 6756 ua1=num.zeros((stations,time_step_count),num.Float) 6757 ua1[0]=3*num.cos(times_ref) 6758 ua1[1]=2*num.sin(times_ref0.7) 6759 ua1[2]=num.arange(3*time_step_count,4*time_step_count) 6760 ua1[4]=2*num.ones(time_step_count) 6761 6762 va1=num.zeros((stations,time_step_count),num.Float) 6763 va1[0]=2*num.cos(times_ref0.87) 6764 va1[1]=3*num.ones(time_step_count) 6765 va1[3]=2*num.sin(times_ref0.71) 6766 #print "va1[0]", va1[0] # The 8th element is what will go bad. 6767 # Ensure data used to write mux file to be zero when gauges are 6768 # not recording 6769 for i in range(stations): 6770 # For each point 6771 for j in range(0, first_tstep[i]1) + range(last_tstep[i], 6772 time_step_count): 6773 # For timesteps before and after recording range 6774 ha1[i][j] = ua1[i][j] = va1[i][j] = 0.0 6775 6776 6777 #print 'Second station to be written to MUX' 6778 #print 'ha', ha1[0,:] 6779 #print 'ua', ua1[0,:] 6780 #print 'va', va1[0,:] 6781 6782 # Write second mux file to be combined by urs2sts 6783 base_nameII, filesII = self.write_mux2(lat_long_points, 6784 time_step_count, time_step, 6785 first_tstep, last_tstep, 6786 depth=gauge_depth, 6787 ha=ha1, 6788 ua=ua1, 6789 va=va1) 6790 #print "filesII", filesII 6791 6792 6793 6794 6795 # Read mux file back and verify it's correcness 6796 6797 #################################################### 6798 # FIXME (Ole): This is where the test should 6799 # verify that the MUX files are correct. 6800 6801 #JJ: It appears as though 6802 #that certain quantities are not being stored with enough precision 6803 #inn muxfile or more likely that they are being cast into a 6804 #lower precision when read in using read_mux2 Time step and q_time 6805 # are equal but only to approx 1e7 6806 #################################################### 6807 6808 #define information as it should be stored in mus2 files 6809 points_num=len(lat_long_points) 6810 depth=gauge_depth 6811 ha=ha1 6812 ua=ua1 6813 va=va1 6814 6815 quantities = ['HA','UA','VA'] 6816 mux_names = [WAVEHEIGHT_MUX2_LABEL, 6817 EAST_VELOCITY_MUX2_LABEL, 6818 NORTH_VELOCITY_MUX2_LABEL] 6819 quantities_init = [[],[],[]] 6820 latlondeps = [] 6821 #irrelevant header information 6822 ig=ilon=ilat=0 6823 mcolat=mcolon=centerlat=centerlon=offset=az=baz=id=0.0 6824 # urs binary is latitude fastest 6825 for i,point in enumerate(lat_long_points): 6826 lat = point[0] 6827 lon = point[1] 6828 _ , e, n = redfearn(lat, lon) 6829 if depth is None: 6830 this_depth = n 6831 else: 6832 this_depth = depth[i] 6833 latlondeps.append([lat, lon, this_depth]) 6834 6835 if ha is None: 6836 this_ha = e 6837 quantities_init[0].append(num.ones(time_step_count, 6838 num.Float)*this_ha) # HA 6839 else: 6840 quantities_init[0].append(ha[i]) 6841 if ua is None: 6842 this_ua = n 6843 quantities_init[1].append(num.ones(time_step_count, 6844 num.Float)*this_ua) # UA 6845 else: 6846 quantities_init[1].append(ua[i]) 6847 if va is None: 6848 this_va = e 6849 quantities_init[2].append(num.ones(time_step_count, 6850 num.Float)*this_va) # 6851 else: 6852 quantities_init[2].append(va[i]) 6853 6854 for i, q in enumerate(quantities): 6855 #print 6856 #print i, q 6857 6858 q_time = num.zeros((time_step_count, points_num), num.Float) 6859 quantities_init[i] = ensure_numeric(quantities_init[i]) 6860 for time in range(time_step_count): 6861 #print i, q, time, quantities_init[i][:,time] 6862 q_time[time,:] = quantities_init[i][:,time] 6863 #print i, q, time, q_time[time, :] 6864 6865 6866 filename = base_nameII + mux_names[i] 6867 f = open(filename, 'rb') 6868 if self.verbose: print 'Reading' + filename 6869 assert abs(points_numunpack('i',f.read(4))[0])<epsilon 6870 #write mux 2 header 6871 for latlondep in latlondeps: 6872 assert abs(latlondep[0]unpack('f',f.read(4))[0])<epsilon 6873 assert abs(latlondep[1]unpack('f',f.read(4))[0])<epsilon 6874 assert abs(mcolatunpack('f',f.read(4))[0])<epsilon 6875 assert abs(mcolonunpack('f',f.read(4))[0])<epsilon 6876 assert abs(igunpack('i',f.read(4))[0])<epsilon 6877 assert abs(ilonunpack('i',f.read(4))[0])<epsilon 6878 assert abs(ilatunpack('i',f.read(4))[0])<epsilon 6879 assert abs(latlondep[2]unpack('f',f.read(4))[0])<epsilon 6880 assert abs(centerlatunpack('f',f.read(4))[0])<epsilon 6881 assert abs(centerlonunpack('f',f.read(4))[0])<epsilon 6882 assert abs(offsetunpack('f',f.read(4))[0])<epsilon 6883 assert abs(azunpack('f',f.read(4))[0])<epsilon 6884 assert abs(bazunpack('f',f.read(4))[0])<epsilon 6885 6886 x = unpack('f', f.read(4))[0] 6887 #print time_step 6888 #print x 6889 assert abs(time_stepx)<epsilon 6890 assert abs(time_step_countunpack('i',f.read(4))[0])<epsilon 6891 for j in range(4): # identifier 6892 assert abs(idunpack('i',f.read(4))[0])<epsilon 6893 6894 #first_tstep=1 6895 #last_tstep=time_step_count 6896 for i,latlondep in enumerate(latlondeps): 6897 assert abs(first_tstep[i]unpack('i',f.read(4))[0])<epsilon 6898 for i,latlondep in enumerate(latlondeps): 6899 assert abs(last_tstep[i]unpack('i',f.read(4))[0])<epsilon 6900 6901 # Find when first station starts recording 6902 min_tstep = min(first_tstep) 6903 # Find when all stations have stopped recording 6904 max_tstep = max(last_tstep) 6905 6906 #for time in range(time_step_count): 6907 for time in range(min_tstep1,max_tstep): 6908 assert abs(time*time_stepunpack('f',f.read(4))[0])<epsilon 6909 for point_i in range(points_num): 6910 if time+1>=first_tstep[point_i] and time+1<=last_tstep[point_i]: 6911 x = unpack('f',f.read(4))[0] 6912 #print time, x, q_time[time, point_i] 6913 if q == 'VA': x = x # South is positive in MUX 6914 #print q+" q_time[%d, %d] = %f" %(time, point_i, 6915 #q_time[time, point_i]) 6916 assert abs(q_time[time, point_i]x)<epsilon 6917 6918 f.close() 6919 6920 permutation = ensure_numeric([4,0,2]) 6921 6922 # Create ordering file 6923 # _, ordering_filename = tempfile.mkstemp('') 6924 # order_fid = open(ordering_filename, 'w') 6925 # order_fid.write('index, longitude, latitude\n') 6926 # for index in permutation: 6927 # order_fid.write('%d, %f, %f\n' %(index, 6928 # lat_long_points[index][1], 6929 # lat_long_points[index][0])) 6930 # order_fid.close() 6931 6932 #  6933 # Now read files back and check values 6934 weights = ensure_numeric([1.0]) 6935 6936 # For each quantity read the associated list of source mux2 file with 6937 # extention associated with that quantity 6938 file_params=1*num.ones(3,num.Float) # [nsta,dt,nt] 6939 OFFSET = 5 6940 6941 for j, file in enumerate(filesII): 6942 # Read stage, u, v enumerated as j 6943 #print 'Reading', j, file 6944 #print "file", file 6945 data = read_mux2(1, [file], weights, file_params, 6946 permutation, verbose) 6947 #print str(j) + "data", data 6948 6949 #print 'Data received by Python' 6950 #print data[1][8] 6951 number_of_selected_stations = data.shape[0] 6952 #print "number_of_selected_stations", number_of_selected_stations 6953 #print "stations", stations 6954 6955 # Index where data ends and parameters begin 6956 parameters_index = data.shape[1]OFFSET 6957 6958 for i in range(number_of_selected_stations): 6959 6960 #print i, parameters_index 6961 if j == 0: 6962 assert num.allclose(data[i][:parameters_index], 6963 ha1[permutation[i], :]) 6964 6965 if j == 1: assert num.allclose(data[i][:parameters_index], ua1[permutation[i], :]) 6966 if j == 2: 6677 6967 assert num.allclose(data[i][:parameters_index], va1[permutation[i], :]) 6678 6968 6969 self.delete_mux(filesII) 6679 6970 6680 6971 def test_urs2sts0(self): … … 6857 7148 assert num.allclose([x[i],y[i]], [e,n]) 6858 7149 assert zone==1 7150 7151 self.delete_mux(files) 6859 7152 6860 7153 … … 6921 7214 assert num.allclose([x[i],y[i]], [e,n]) 6922 7215 assert zone==geo_reference.zone 7216 7217 self.delete_mux(files) 6923 7218 6924 7219 … … 7651 7946 raise Exception, msg 7652 7947 7948 7949 self.delete_mux(filesI) 7950 self.delete_mux(filesII) 7653 7951 7654 7952 … … 7785 8083 7786 8084 7787 # Write second mux file to be combined by urs2sts 8085 # Write second mux file to be combined by urs2sts 7788 8086 base_nameII, filesII = self.write_mux2(lat_long_points, 7789 8087 time_step_count, time_step, … … 8033 8331 os.remove(sts_file) 8034 8332 8035 8036 8037 8333 # 8038 8334 # Then read the mux files together and test … … 8140 8436 8141 8437 fid.close() 8142 self.delete_mux(filesI)8143 self.delete_mux(filesII)8144 8438 os.remove(sts_file) 8145 8146 8147 8439 8148 8440 # … … 8153 8445 stage_man = weights[0]*(stage0tide) + weights[1]*(stage1tide) + tide 8154 8446 assert num.allclose(stage_man, stage) 8155 8156 8157 8158 8159 8160 8447 8161 8448 … … 8284 8571 8285 8572 assert num.allclose(domain_fbound.quantities['xmomentum'].vertex_values, 8286 domain_drchlt.quantities['xmomentum'].vertex_values) 8573 domain_drchlt.quantities['xmomentum'].vertex_values) 8287 8574 8288 8575 assert num.allclose(domain_fbound.quantities['ymomentum'].vertex_values, 8289 domain_drchlt.quantities['ymomentum'].vertex_values) 8576 domain_drchlt.quantities['ymomentum'].vertex_values) 8290 8577 8291 8578 8292 8579 os.remove(sts_file+'.sts') 8293 8580 os.remove(meshname) 8294 8295 8296 8297 8581 8298 8582 … … 8415 8699 8416 8700 8417 8418 8419 8420 8421 8422 8423 8424 8701 domain_drchlt = Domain(meshname) 8425 8702 domain_drchlt.set_quantity('stage', tide) … … 8443 8720 8444 8721 assert num.allclose(domain_fbound.quantities['xmomentum'].vertex_values, 8445 domain_drchlt.quantities['xmomentum'].vertex_values) 8722 domain_drchlt.quantities['xmomentum'].vertex_values) 8446 8723 8447 8724 assert num.allclose(domain_fbound.quantities['ymomentum'].vertex_values, 8448 domain_drchlt.quantities['ymomentum'].vertex_values) 8449 8725 domain_drchlt.quantities['ymomentum'].vertex_values) 8450 8726 8451 8727 os.remove(sts_file+'.sts') 8452 8728 os.remove(meshname) 8453 8454 8455 8456 8729 8457 8730 … … 11365 11638 #suite = unittest.makeSuite(Test_Data_Manager,'test_urs2sts_ordering_different_sources') 11366 11639 11367 # FIXME (Ole): This is the test that fails under Windows 11368 #suite = unittest.makeSuite(Test_Data_Manager,'test_read_mux_platform_problem2') 11640 #suite = unittest.makeSuite(Test_Data_Manager,'test_read_mux_platform_problem3') 11369 11641 #suite = unittest.makeSuite(Test_Data_Manager,'test_file_boundary_stsIV') 11370 11642
Note: See TracChangeset
for help on using the changeset viewer.