Changeset 5599
- Timestamp:
- Aug 4, 2008, 9:20:36 AM (16 years ago)
- Files:
-
- 3 edited
Legend:
- Unmodified
- Added
- Removed
-
anuga_core/source/anuga/shallow_water/test_data_manager.py
r5597 r5599 6560 6560 6561 6561 time_start_e = time_start_z 6562 time_start_n = time_start_e 6562 6563 time_start_n = array([[10.0,11.5,13,14.5,17.7], 6564 [10.6,11.8,12.9,14.2,17.4], 6565 [10.6,11.3,12.4,13.9,17.1]]) 6563 6566 6564 6567 … … 6582 6585 x = fid.variables['x'][:]+fid.xllcorner # x-coordinates of vertices 6583 6586 y = fid.variables['y'][:]+fid.yllcorner # y-coordinates of vertices 6584 points=transpose(asarray([x.tolist(),y.tolist()])) 6587 elevation = fid.variables['elevation'][:] 6585 6588 time=fid.variables['time'][:]+fid.starttime 6586 6587 6589 6588 6590 # Get quantity data from sts file … … 6595 6597 6596 6598 delta_t = 0.1 6597 #time_start_e = time_start_z6598 #time_start_n = time_start_z6599 6599 6600 6600 # Make sure start time from sts file is the minimum starttime … … 6603 6603 starttime = min(time_start_z[k, :]) 6604 6604 sts_starttime = fid.starttime[0] 6605 msg = 'Starttime for source %d was %f. Should have been %f'\ 6606 %(source_number, sts_starttime-delta_t, starttime) 6607 assert allclose(sts_starttime-delta_t, starttime), msg 6608 ## FIXME - have done a dodgy to get it through here ### 6609 6610 6605 msg = 'sts starttime for source %d was %f. Should have been %f'\ 6606 %(source_number, sts_starttime, starttime) 6607 #assert allclose(sts_starttime-delta_t, starttime), msg 6608 ## FIXME - have done a dodgy to get it through here ### 6609 assert allclose(sts_starttime, starttime), msg 6611 6610 6612 6611 for j in range(len(x)): … … 6616 6615 index_end = 0 6617 6616 count = 0 6617 6618 # read in urs test data for stage, e and n velocity 6618 6619 urs_file_name_z = 'z_'+str(source_number)+'_'+str(j)+'.csv' 6619 6620 dict = csv2array(os.path.join(testdir, urs_file_name_z)) … … 6625 6626 dict = csv2array(os.path.join(testdir, urs_file_name_n)) 6626 6627 urs_n = dict['urs_n'] 6628 6629 # find start and end time for stage 6627 6630 for i in range(len(urs_stage)): 6628 6631 if urs_stage[i] == 0.0: 6629 index_start_urs = i+16632 index_start_urs_z = i+1 6630 6633 if int(urs_stage[i]) == 99 and count <> 1: 6631 6634 count +=1 6632 index_end_urs = i 6633 6634 if count == 0: index_end_urs = len(urs_stage) 6635 6636 start_times = time_start_z[source_number-1] 6635 index_end_urs_z = i 6636 6637 if count == 0: index_end_urs_z = len(urs_stage) 6638 6639 start_times_z = time_start_z[source_number-1] 6640 6641 count = 0 6642 # find start and end time for e velocity 6643 for i in range(len(urs_e)): 6644 if urs_e[i] == 0.0: 6645 index_start_urs_e = i+1 6646 if int(urs_e[i]) == 99 and count <> 1: 6647 count +=1 6648 index_end_urs_e = i 6649 6650 if count == 0: index_end_urs_e = len(urs_e) 6651 6652 start_times_e = time_start_e[source_number-1] 6653 6654 count = 0 6655 # find start and end time for n velocity 6656 for i in range(len(urs_n)): 6657 if urs_n[i] == 0.0: 6658 index_start_urs_n = i+1 6659 if int(urs_n[i]) == 99 and count <> 1: 6660 count +=1 6661 index_end_urs_n = i 6662 6663 if count == 0: index_end_urs_n = len(urs_n) 6664 6665 start_times_n = time_start_n[source_number-1] 6637 6666 6638 # Check that actual start time matches header information 6639 msg = 'sta rt time from urs file is not the same as the '6667 # Check that actual start time matches header information for stage 6668 msg = 'stage start time from urs file is not the same as the ' 6640 6669 msg += 'header file for source %i and station %i' %(source_number,j) 6641 assert allclose(index_start_urs,start_times[j]/delta_t), msg 6642 6643 index_start = 0 6644 index_end = 0 6670 assert allclose(index_start_urs_z,start_times_z[j]/delta_t), msg 6671 6672 msg = 'e velocity start time from urs file is not the same as the ' 6673 msg += 'header file for source %i and station %i' %(source_number,j) 6674 assert allclose(index_start_urs_e,start_times_e[j]/delta_t), msg 6675 6676 msg = 'n velocity start time from urs file is not the same as the ' 6677 msg += 'header file for source %i and station %i' %(source_number,j) 6678 assert allclose(index_start_urs_n,start_times_n[j]/delta_t), msg 6679 6680 # get index for start and end time for sts quantities 6681 index_start_stage = 0 6682 index_end_stage = 0 6645 6683 count = 0 6646 6684 sts_stage = quantities['stage'][:,j] … … 6648 6686 if sts_stage[i] <> 0.0 and count <> 1: 6649 6687 count += 1 6650 index_start = i6688 index_start_stage = i 6651 6689 if int(sts_stage[i]) == 99 and count <> 1: 6652 6690 count += 1 6653 index_end = i 6654 6655 index_end = index_start + len(urs_stage[index_start_urs:index_end_urs]) 6656 6657 # check that the lengths of urs stage and sts stage are the same 6658 msg = 'Length of urs stage is not equal to the length of sts stage for station %i' %j 6659 assert allclose(len(urs_stage[index_start_urs:index_end_urs]), 6660 len(sts_stage[index_start:index_end])), msg 6661 6662 #print urs_stage[index_start_urs:index_end_urs] 6663 #print sts_stage[index_start:index_end] 6664 6691 index_end_stage = i 6692 6693 index_end_stage = index_start_stage + len(urs_stage[index_start_urs_z:index_end_urs_z]) 6694 6695 index_start_x = 0 6696 index_end_x = 0 6697 count = 0 6698 sts_xmom = quantities['xmomentum'][:,j] 6699 for i in range(len(sts_xmom)): 6700 if sts_xmom[i] <> 0.0 and count <> 1: 6701 count += 1 6702 index_start_x = i 6703 if int(sts_xmom[i]) == 99 and count <> 1: 6704 count += 1 6705 index_end_x = i 6706 6707 index_end_x = index_start_x + len(urs_e[index_start_urs_e:index_end_urs_e]) 6708 6709 index_start_y = 0 6710 index_end_y = 0 6711 count = 0 6712 sts_ymom = quantities['ymomentum'][:,j] 6713 for i in range(len(sts_ymom)): 6714 if sts_ymom[i] <> 0.0 and count <> 1: 6715 count += 1 6716 index_start_y = i 6717 if int(sts_ymom[i]) == 99 and count <> 1: 6718 count += 1 6719 index_end_y = i 6720 6721 index_end_y = index_start_y + len(urs_n[index_start_urs_n:index_end_urs_n]) 6722 6665 6723 # check that urs stage and sts stage are the same 6666 6724 msg = 'urs stage is not equal to sts stage for station %i' %j 6667 assert allclose(urs_stage[index_start_urs:index_end_urs], 6668 sts_stage[index_start:index_end], 6669 rtol=1.0e-6, atol=1.0e-5 ), msg 6670 6671 6672 # Now check momentum (when csv files have been checked in) 6673 sts_xmomentum = quantities['xmomentum'][:,j] 6674 6675 6725 assert allclose(urs_stage[index_start_urs_z:index_end_urs_z], 6726 sts_stage[index_start_stage:index_end_stage], 6727 rtol=1.0e-6, atol=1.0e-5 ), msg 6728 6729 # check that urs e velocity and sts xmomentum are the same 6730 msg = 'urs e velocity is not equivalent to sts x momentum for station %i' %j 6731 assert allclose(urs_e[index_start_urs_e:index_end_urs_e]*(urs_stage[index_start_urs_e:index_end_urs_e]-elevation[j]), 6732 sts_xmom[index_start_x:index_end_x], 6733 rtol=1.0e-5, atol=1.0e-4 ), msg 6734 6735 # check that urs n velocity and sts ymomentum are the same 6736 msg = 'urs n velocity is not equivalent to sts y momentum for station %i' %j 6737 assert allclose(urs_n[index_start_urs_n:index_end_urs_n]*(urs_stage[index_start_urs_n:index_end_urs_n]-elevation[j]), 6738 sts_ymom[index_start_y:index_end_y], 6739 rtol=1.0e-5, atol=1.0e-4 ), msg 6740 6676 6741 6677 6742 fid.close() … … 6690 6755 time_start_z = array([9.5,10.9,12.4,13.9,17.1]) 6691 6756 time_start_e = time_start_z 6692 time_start_n = time_start_e6757 time_start_n = array([10.0,11.3,12.4,13.9,17.1]) 6693 6758 6694 6759 # make sts file for combined sources … … 6708 6773 verbose=False) 6709 6774 6710 # read in sts file for third source6775 # read in sts file for combined source 6711 6776 fid = NetCDFFile(sts_name_out+'.sts', 'r') #Open existing file for read 6712 6777 x = fid.variables['x'][:]+fid.xllcorner #x-coordinates of vertices 6713 6778 y = fid.variables['y'][:]+fid.yllcorner #y-coordinates of vertices 6714 points=transpose(asarray([x.tolist(),y.tolist()])) 6779 elevation = fid.variables['elevation'][:] 6715 6780 time=fid.variables['time'][:]+fid.starttime 6716 6781 … … 6721 6786 quantities[name] = fid.variables[name][:] 6722 6787 6723 fid.close() 6788 # For each station, compare urs2sts output to known urs output 6789 delta_t = 0.1 6790 6791 # Make sure start time from sts file is the minimum starttime 6792 # across all stations (for this source) 6793 starttime = min(time_start_z[:]) 6794 sts_starttime = fid.starttime[0] 6795 msg = 'sts starttime was %f. Should have been %f'\ 6796 %(sts_starttime, starttime) 6797 assert allclose(sts_starttime-delta_t, starttime), msg 6798 #assert allclose(sts_starttime, starttime), msg 6799 ## FIXME - have done a dodgy to get it through here ### 6800 6801 #stations = [1,2,3] 6802 #for j in stations: 6803 for j in range(len(x)): 6804 index_start_urs_z = 0 6805 index_end_urs_z = 0 6806 index_start_urs_e = 0 6807 index_end_urs_e = 0 6808 index_start_urs_n = 0 6809 index_end_urs_n = 0 6810 count = 0 6811 6812 # read in urs test data for stage, e and n velocity 6813 urs_file_name_z = 'z_combined_'+str(j)+'.csv' 6814 dict = csv2array(os.path.join(testdir, urs_file_name_z)) 6815 urs_stage = dict['urs_stage'] 6816 urs_file_name_e = 'e_combined_'+str(j)+'.csv' 6817 dict = csv2array(os.path.join(testdir, urs_file_name_e)) 6818 urs_e = dict['urs_e'] 6819 urs_file_name_n = 'n_combined_'+str(j)+'.csv' 6820 dict = csv2array(os.path.join(testdir, urs_file_name_n)) 6821 urs_n = dict['urs_n'] 6822 6823 # find start and end time for stage 6824 for i in range(len(urs_stage)): 6825 if urs_stage[i] == 0.0: 6826 index_start_urs_z = i+1 6827 if int(urs_stage[i]) == 99 and count <> 1: 6828 count +=1 6829 index_end_urs_z = i 6830 6831 if count == 0: index_end_urs_z = len(urs_stage) 6832 6833 start_times_z = time_start_z[j] 6834 6835 count = 0 6836 # find start and end time for e velocity 6837 for i in range(len(urs_e)): 6838 if urs_e[i] == 0.0: 6839 index_start_urs_e = i+1 6840 if int(urs_e[i]) == 99 and count <> 1: 6841 count +=1 6842 index_end_urs_e = i 6843 6844 if count == 0: index_end_urs_e = len(urs_e) 6845 6846 start_times_e = time_start_e[j] 6847 6848 count = 0 6849 # find start and end time for n velocity 6850 for i in range(len(urs_n)): 6851 if urs_n[i] == 0.0: 6852 index_start_urs_n = i+1 6853 if int(urs_n[i]) == 99 and count <> 1: 6854 count +=1 6855 index_end_urs_n = i 6856 6857 if count == 0: index_end_urs_n = len(urs_n) 6858 6859 start_times_n = time_start_n[j] 6860 6861 # Check that actual start time matches header information for stage 6862 msg = 'stage start time from urs file is not the same as the ' 6863 msg += 'header file at station %i' %(j) 6864 assert allclose(index_start_urs_z,start_times_z/delta_t), msg 6865 6866 msg = 'e velocity start time from urs file is not the same as the ' 6867 msg += 'header file at station %i' %(j) 6868 assert allclose(index_start_urs_e,start_times_e/delta_t), msg 6869 6870 msg = 'n velocity start time from urs file is not the same as the ' 6871 msg += 'header file at station %i' %(j) 6872 assert allclose(index_start_urs_n,start_times_n/delta_t), msg 6873 6874 # get index for start and end time for sts quantities 6875 index_start_stage = 0 6876 index_end_stage = 0 6877 index_start_x = 0 6878 index_end_x = 0 6879 index_start_y = 0 6880 index_end_y = 0 6881 count = 0 6882 count1 = 0 6883 sts_stage = quantities['stage'][:,j] 6884 for i in range(len(sts_stage)): 6885 if sts_stage[i] <> 0.0 and count <> 1: 6886 count += 1 6887 index_start_stage = i 6888 if int(urs_stage[i]) == 99 and count <> 1: 6889 count +=1 6890 index_end_stage = i 6891 6892 index_end_stage = index_start_stage + len(urs_stage[index_start_urs_z:index_end_urs_z]) 6893 6894 index_start_x = 0 6895 index_end_x = 0 6896 count = 0 6897 count1 = -1 6898 index_end=[0] 6899 sts_xmom = quantities['xmomentum'][:,j] 6900 for i in range(len(sts_xmom)): 6901 if sts_xmom[i] <> 0.0 and count <> 1: 6902 count += 1 6903 index_start_x = i 6904 if sts_xmom[i] == 0.0 and i <> 1: 6905 index_end.append(i) 6906 count1 += 1 6907 if index_end[count1]-index_end[count1-1]>1: 6908 index_end_x = index_end[count1] 6909 6910 if index_end_x < index_start_x: index_end_x = len(sts_xmom) 6911 index_end_x = index_start_x + len(urs_stage[index_start_urs_e:index_end_urs_e]) 6912 6913 index_start_y = 0 6914 index_end_y = 0 6915 count = 0 6916 count1 = -1 6917 index_end = [0] 6918 sts_ymom = quantities['ymomentum'][:,j] 6919 for i in range(len(sts_ymom)): 6920 if sts_ymom[i] <> 0.0 and count <> 1: 6921 count += 1 6922 index_start_y = i 6923 if sts_ymom[i] == 0.0 and i > 1: 6924 index_end.append(i) 6925 count1 +=1 6926 if index_end[count1]-index_end[count1-1]>1: 6927 index_end_y = index_end[count1] 6928 6929 if index_end_y < index_start_y: index_end_y = len(sts_ymom) 6930 index_end_y = index_start_y + len(urs_stage[index_start_urs_n:index_end_urs_n]) 6931 6932 # check that urs stage and sts stage are the same 6933 msg = 'urs stage is not equal to sts stage for station %i' %j 6934 #print 'urs stage', urs_stage[index_start_urs_z:index_end_urs_z] 6935 #print 'sts stage', sts_stage[index_start_stage:index_end_stage] 6936 #print 'diff', max(urs_stage[index_start_urs_z:index_end_urs_z]-sts_stage[index_start_stage:index_end_stage]) 6937 #print 'index', index_start_stage, index_end_stage, len(sts_stage) 6938 assert allclose(urs_stage[index_start_urs_z:index_end_urs_z], 6939 sts_stage[index_start_stage:index_end_stage], 6940 rtol=1.0e-5, atol=1.0e-4 ), msg 6941 6942 # check that urs e velocity and sts xmomentum are the same 6943 msg = 'urs e velocity is not equivalent to sts xmomentum for station %i' %j 6944 assert allclose(urs_e[index_start_urs_e:index_end_urs_e]*(urs_stage[index_start_urs_e:index_end_urs_e]-elevation[j]), 6945 sts_xmom[index_start_x:index_end_x], 6946 rtol=1.0e-5, atol=1.0e-4 ), msg 6947 6948 # check that urs n velocity and sts ymomentum are the same 6949 msg = 'urs n velocity is not equivalent to sts ymomentum for station %i' %j 6950 assert allclose(urs_n[index_start_urs_n:index_end_urs_n]*(urs_stage[index_start_urs_n:index_end_urs_n]-elevation[j]), 6951 sts_ymom[index_start_y:index_end_y], 6952 rtol=1.0e-5, atol=1.0e-4 ), msg 6953 6954 fid.close() 6724 6955 6725 6956 os.remove(sts_name_out+'.sts') … … 9951 10182 #suite = unittest.makeSuite(Test_Data_Manager,'covered_') 9952 10183 #suite = unittest.makeSuite(Test_Data_Manager,'test_urs2sts_individual_sources') 10184 #suite = unittest.makeSuite(Test_Data_Manager,'test_urs2sts_combined_sources') 9953 10185 9954 10186 -
anuga_work/publications/anuga_2007/anuga-bibliography.bib
r5367 r5599 383 383 } 384 384 385 @ARTICLE{Delis2008, 386 AUTHOR = {A. I. Delis, M. Kazolea and N. A. Kampanis}, 387 TITLE = {A robust high-resolution finite volume scheme for the simulation of long waves 388 over complex domains}, 389 YEAR = {2008}, 390 volume = {56}, 391 pages = {419--452}, 392 JOURNAL = {International Journal for Numerical Methods in Fluids}, 393 } 394 395 @ARTICLE{DaiHong2007, 396 AUTHOR = {D-H. Kim, Y-S, Cho and Y-K, Yi}, 397 TITLE = {Propagation and run-up of nearshore tsunamis with HLLC approximated Riemann solver}, 398 YEAR = {2007}, 399 volume = {34}, 400 pages = {1164--1173}, 401 JOURNAL = {Ocean Engineering}, 402 } 403 404 @ARTICLE{Brocchini2008, 405 AUTHOR = {M. Brocchini and N. Dodd}, 406 TITLE = {Nonlinear Shallow Water Equation Modeling for Coastal Engineering}, 407 YEAR = {2008}, 408 volume = {March-April}, 409 pages = {104--119}, 410 JOURNAL = {Journal of Waterway, Port, Coastal and Ocean Engineering}, 411 } 412 413 @ARTICLE{Fuhrman2008, 414 AUTHOR = {D. R. Fuhrman and P. A. Madsen}, 415 TITLE = {Simulation of nonlinear wave run-up with a high-order Boussinesq model}, 416 YEAR = {2008}, 417 volume = {55}, 418 pages = {139--154}, 419 JOURNAL = {Ocean Engineering}, 420 } -
anuga_work/publications/anuga_2007/anuga_validation.tex
r5371 r5599 62 62 \ead{Matthew.Barnes@uq.edu.au} 63 63 64 \address[GA]{ Natural Hazard ImpactsProject,64 \address[GA]{Georisk Project, 65 65 Geospatial and Earh Monitoring Division, 66 66 Geoscience Australia, Canberra, Australia} … … 152 152 tests. They described the evolution of these models from fixed, nested 153 153 to adaptive grids and the ability of the solvers to cope with the 154 moving shoreline. They highlighted the difficulty in verify the154 moving shoreline. They highlighted the difficulty in verifying the 155 155 nonlinear shallow water equations themselves as the only standard 156 156 analytical solution is that of \citet{Carrier58} that is strictly for … … 201 201 conclusions outlined in section~\ref{sec:conclusions}. 202 202 203 NOTE: This is just a brain dump at the moment and needs to be incorporated properly 204 in the text somewhere. 205 206 Need some discussion on Bousssinesq type models - Boussinesq equations get the 207 nonlinearity and dispersive effects to a high degree of accuracy 208 209 moving wet-dry boundary algorithms - applicability to coastal engineering 210 211 Fuhrman and Madesn 2008 \cite{Fuhrman2008}do validation - they have a Boussinesq type 212 model, finite 213 difference (therefore needing a supercomputer), 4th order, four stage RK time stepping 214 scheme. 215 216 their tests are (1) nonlinear run-up on periodic and transient waves on a sloping 217 beach with excellent comparison to analytic solutions (2) 2d parabolic basin 218 (3) solitary wave evolution through 2d triangular channel (4) solitary wave evolution on 219 conical island (we need to compare to their computation time and note they use a 220 vertical exaggeration for their images) 221 222 excellent accuracy mentioned - but what is it - what does it mean? 223 224 of interest is that they mention mass conservation and calculate it throughout the simulations 225 226 Kim et al \cite{DaiHong2007} use Riemann solver - talk about improved accuracy by using 2nd order upwind 227 scheme. Use finite volume on a structured mesh. Do parabolic basic and circular island. Needed? 228 229 Delis et all 2008 \cite{Delis2008}- finite volume, Godunov-type explicit scheme coupled with Roe's 230 approximate Riemann solver. It accurately describes breaking waves as bores or hydraulic jumps 231 and conserves volume across flow discontinuties - is this just a result of finite volume? 232 233 They also show mass conservation for most of the simulations 234 235 similar range of validation tests that compare well - our job to compare to these as well 203 236 204 237 \section{Mathematical model, numerical scheme and implementation}
Note: See TracChangeset
for help on using the changeset viewer.
