source: anuga_validation/okushiri_2005/compare_timeseries.py @ 3913

Last change on this file since 3913 was 3913, checked in by ole, 17 years ago

Got similar number of triangles as Aug 2005 okushiri validation

File size: 11.2 KB
RevLine 
[3845]1"""Verify that simulation produced by ANUGA compares to published
2validation timeseries ch5, ch7 and ch9 as well as the boundary timeseries.
3
4RMS norm is printed and plots are produced.
[2229]5"""
6
[3913]7from Numeric import allclose, argmin, argmax
[3850]8from Scientific.IO.NetCDF import NetCDFFile
[3845]9
[3850]10from anuga.abstract_2d_finite_volumes.util import file_function
11from anuga.utilities.numerical_tools import\
12     ensure_numeric, cov, get_machine_precision
[3845]13
[3850]14import project
[2229]15
[3850]16try:
17    from pylab import ion, hold, plot, title, legend, xlabel, ylabel, savefig
18except:
19    plotting = False
20else:
21    plotting = True
[2229]22
[3861]23#plotting = False
[2229]24
[3850]25#-------------------------
26# Basic data
27#-------------------------
[3845]28
[3850]29finaltime = 22.5
30timestep = 0.05
[2229]31
[3850]32gauge_locations = [[0.000, 1.696]] # Boundary gauge
33gauge_locations += [[4.521, 1.196],  [4.521, 1.696],  [4.521, 2.196]] #Ch 5-7-9
[2229]34gauge_names = ['Boundary', 'ch5', 'ch7', 'ch9']
35
[3850]36validation_data = {}
37for key in gauge_names:
38    validation_data[key] = []
[2229]39
[3850]40   
[3861]41#expected_covariances = {'Boundary': 5.288392008865989e-05,
42#                        'ch5': 1.166748190444681e-04,
43#                        'ch7': 1.121816242516758e-04,
44#                        'ch9': 1.249543278366778e-04}
[2229]45
[3854]46# old limiters
[3878]47#expected_covariances = {'Boundary': 5.288601162783020386e-05,
48#                        'ch5': 1.167001054284431472e-04,
49#                        'ch7': 1.121474766904651861e-04,
50#                        'ch9': 1.249244820847215335e-04}
51#
52#expected_differences = {'Boundary': 8.361144081847830638e-04,
53#                        'ch5': 3.423673831653336816e-03,
54#                        'ch7': 2.799962153549145211e-03,
55#                        'ch9': 3.198560464876740433e-03}
[2229]56
[3854]57
[3883]58#expected_covariances = {'Boundary': 5.288392008865989237e-05,
59#                        'ch5': 1.166748190444680592e-04,
60#                        'ch7': 1.121816242516757850e-04,
61#                        'ch9': 1.249543278366777640e-04}
62#
63#expected_differences = {'Boundary':  8.373150808730501615e-04,
64#                        'ch5': 3.425914311580337875e-03,
65#                        'ch7': 2.802327594773105189e-03,
66#                        'ch9': 3.198733498646373370e-03}
[3861]67
68
[3913]69#expected_covariances = {'Boundary': 5.299487474489660856e-05,
70#                        'ch5': 1.169650160375980370e-04,
71#                        'ch7': 1.123947836450141360e-04,
72#                        'ch9': 1.248329330436513066e-04}
73#
74#expected_differences = {'Boundary': 8.269932106586290379e-04,
75#                        'ch5': 3.411769502897898498e-03,
76#                        'ch7': 2.777024544282339583e-03,
77#                        'ch9': 3.183530359567784788e-03}
78#
[3878]79
80
[3913]81expected_covariance = {'Boundary': 5.265850426352561254e-05, 
82                       'ch5': 1.171913883449135149e-04,
83                       'ch7': 1.130834789656423632e-04,
84                       'ch9': 1.257648236914514148e-04}
[3883]85
[3913]86expected_difference = {'Boundary': 9.188561500629542633e-04,
87                       'ch5': 3.439740513194261117e-03,
88                       'ch7': 2.868914369295409543e-03,
89                       'ch9': 3.287722754612562841e-03}
[3883]90
[3913]91expected_maximum_diff = {'Boundary': 4.865498851439054029e-05,
92                         'ch5': 1.808007760316858448e-03,
93                         'ch7': 5.464055727564115506e-04,
94                         'ch9': 1.856547605566957748e-03}
[3883]95
[3913]96expected_minimum_diff = {'Boundary': 2.293007809784416290e-04,
97                         'ch5': 1.494336068527772621e-04,
98                         'ch7': 4.515664278503157131e-03,
99                         'ch9': 3.406785970601119290e-03}
[3883]100
[3913]101expected_argmax_timelag = {'Boundary': 3.000000000000007105e-01,
102                           'ch5': 9.999999999999786837e-02,
103                           'ch7': 0.000000000000000000e+00,
104                           'ch9': 0.000000000000000000e+00}
105
106expected_argmin_timelag = {'Boundary': 3.000000000000007105e-01,
107                           'ch5': 8.499999999999996447e-01,
108                           'ch7': 4.500000000000010658e-01,
109                           'ch9': 9.499999999999992895e-01}
110
111
112# Results from lwru2_variable_mesh.sww
113#
114#Validating Boundary
115#Covariance = 5.265850426352561254e-05
116#Accumulated difference = 9.188561500629542633e-04
117#Difference in maxima = 4.865498851439054029e-05
118#Difference in minima = 2.293007809784416290e-04
119#Timelag between maxima = 3.000000000000007105e-01
120#Timelag between minima = 3.000000000000007105e-01
121#Next
122#
123#Validating ch5
124#Covariance = 1.171913883449135149e-04
125#Accumulated difference = 3.439740513194261117e-03
126#Difference in maxima = 1.808007760316858448e-03
127#Difference in minima = 1.494336068527772621e-04
128#Timelag between maxima = 9.999999999999786837e-02
129#Timelag between minima = 8.499999999999996447e-01
130#Next#
131#
132#Validating ch7
133#Covariance = 1.130834789656423632e-04
134#Accumulated difference = 2.868914369295409543e-03
135#Difference in maxima = 5.464055727564115506e-04
136#Difference in minima = 4.515664278503157131e-03
137#Timelag between maxima = 0.000000000000000000e+00
138#Timelag between minima = 4.500000000000010658e-01
139#
140#Next
141#Validating ch9
142#Covariance = 1.257648236914514148e-04
143#Accumulated difference = 3.287722754612562841e-03
144#Difference in maxima = 1.856547605566957748e-03
145#Difference in minima = 3.406785970601119290e-03
146#Timelag between maxima = 0.000000000000000000e+00
147#Timelag between minima = 9.499999999999992895e-01
148
149
150
[3850]151#-------------------------
152# Read validation dataa
153#-------------------------
[2229]154
[3850]155print 'Reading', project.boundary_filename
156fid = NetCDFFile(project.boundary_filename, 'r')
[2229]157input_time = fid.variables['time'][:]
[3850]158validation_data['Boundary'] = fid.variables['stage'][:]
[2229]159
160reference_time = []
[3850]161fid = open(project.validation_filename)
[2229]162lines = fid.readlines()
163fid.close()
[3850]164
[2229]165for i, line in enumerate(lines[1:]):
166    if i == len(input_time): break
[3850]167   
[2229]168    fields = line.split()
169
[3850]170    reference_time.append(float(fields[0]))    # Record reference time
171    for j, key in enumerate(gauge_names[1:]):  # Omit boundary gauge
172        value = float(fields[1:][j])           # Omit time
173        validation_data[key].append(value/100) # Convert cm2m
[2229]174
175
[3850]176# Checks
[2229]177assert reference_time[0] == 0.0
178assert reference_time[-1] == finaltime
179assert allclose(reference_time, input_time)
180
[3850]181for key in gauge_names:
182    validation_data[key] = ensure_numeric(validation_data[key])
[2229]183
184
[3850]185#--------------------------------------------------
186# Read and interpolate model output
187#--------------------------------------------------
[3878]188
189import sys
190if len(sys.argv) > 1:
191    sww_filename = sys.argv[1]
192else:   
193    sww_filename = project.output_filename
194   
[3854]195#f = file_function('okushiri_new_limiters.sww',   #The best so far
[3850]196#f = file_function('okushiri_as2005_with_mxspd=0.1.sww',
[3878]197f = file_function(sww_filename,
[3850]198                  quantities='stage',
199                  interpolation_points=gauge_locations,
200                  use_cache=True,
201                  verbose=True)
[2229]202
203
[3850]204#--------------------------------------------------
205# Compare model output to validation data
206#--------------------------------------------------
207
208eps = get_machine_precision()
[2229]209for k, name in enumerate(gauge_names):
210    sqsum = 0
211    denom = 0
212    model = []
[3850]213    print 
[2229]214    print 'Validating ' + name
[3850]215    observed_timeseries = validation_data[name]
[2229]216    for i, t in enumerate(reference_time):
[3850]217        model.append(f(t, point_id=k)[0])
[2229]218
[3913]219
[3861]220    # Covariance measures   
[3850]221    res = cov(observed_timeseries, model)     
222    print 'Covariance = %.18e' %res
[2229]223
[3913]224    if res < expected_covariance[name]-eps:
225        print '  Result is better than expected by: %.18e'\
226              %(res-expected_covariance[name])
227        print '  Expect = %.18e' %expected_covariance[name]   
228    elif res > expected_covariance[name]+eps:
229        print '  FAIL: Result is worse than expected by: %.18e'\
230                                %(res-expected_covariance[name])
231        print '  Expect = %.18e' %expected_covariance[name]
[3850]232    else:
233        pass
[3913]234   
[2229]235
[3861]236    # Difference measures   
237    res = sum(abs(observed_timeseries-model))/len(model)     
238    print 'Accumulated difference = %.18e' %res
[2229]239
[3913]240    if res < expected_difference[name]-eps:
241        print '  Result is better than expected by: %.18e'\
242              %(res-expected_difference[name])
243        print '  Expect = %.18e' %expected_difference[name]   
244    elif res > expected_difference[name]+eps:
245        print '  FAIL: Result is worse than expected by: %.18e'\
246                                %(res-expected_difference[name])
247        print '  Expect = %.18e' %expected_difference[name]
[3861]248    else:
249        pass
250
251
[3913]252    # Extrema
253    res = abs(max(observed_timeseries)-max(model))
254    print 'Difference in maxima = %.18e' %res
[3861]255
[3913]256    if res < expected_maximum_diff[name]-eps:
257        print '  Result is better than expected by: %.18e'\
258              %(res-expected_maximum_diff[name])
259        print '  Expect = %.18e' %expected_maximum_diff[name]   
260    elif res > expected_maximum_diff[name]+eps:
261        print '  FAIL: Result is worse than expected by: %.18e'\
262                                %(res-expected_maximum_diff[name])
263        print '  Expect = %.18e' %expected_maximum_diff[name]
264    else:
265        pass
266   
[3861]267
[3913]268   
269    res = abs(min(observed_timeseries)-min(model))
270    print 'Difference in minima = %.18e' %res
271
272    if res < expected_minimum_diff[name]-eps:
273        print '  Result is better than expected by: %.18e'\
274              %(res-expected_minimum_diff[name])
275        print '  Expect = %.18e' %expected_minimum_diff[name]   
276    elif res > expected_minimum_diff[name]+eps:
277        print '  FAIL: Result is worse than expected by: %.18e'\
278                                %(res-expected_minimum_diff[name])
279        print '  Expect = %.18e' %expected_minimum_diff[name]
280    else:
281        pass
282   
283
284    # Locations of extrema
285    i0 = argmax(observed_timeseries)
286    i1 = argmax(model)
287    res = abs(reference_time[i1] - reference_time[i0])
288    print 'Timelag between maxima = %.18e' %res
289
290    if res < expected_argmax_timelag[name]-eps:
291        print '  Result is better than expected by: %.18e'\
292              %(res-expected_argmax_timelag[name])
293        print '  Expect = %.18e' %expected_argmax_timelag[name]   
294    elif res > expected_argmax_timelag[name]+eps:
295        print '  FAIL: Result is worse than expected by: %.18e'\
296                                %(res-expected_argmax_timelag[name])
297        print '  Expect = %.18e' %expected_argmax_timelag[name]
298    else:
299        pass
300   
301
302    i0 = argmin(observed_timeseries)
303    i1 = argmin(model)
304    res = abs(reference_time[i1] - reference_time[i0])
305    print 'Timelag between minima = %.18e' %res
306    if res < expected_argmin_timelag[name]-eps:
307        print '  Result is better than expected by: %.18e'\
308              %(res-expected_argmin_timelag[name])
309        print '  Expect = %.18e' %expected_argmin_timelag[name]   
310    elif res > expected_argmin_timelag[name]+eps:
311        print '  FAIL: Result is worse than expected by: %.18e'\
312                                %(res-expected_argmin_timelag[name])
313        print '  Expect = %.18e' %expected_argmin_timelag[name]
314    else:
315        pass
316
317
318
319
[3850]320    if plotting is True:
321        ion()
322        hold(False)
323   
[3913]324        plot(reference_time, validation_data[name], 'r-',
[3850]325             reference_time, model, 'k-')
326        title('Gauge %s' %name)
327        xlabel('time(s)')
328        ylabel('stage (m)')   
329        legend(('Observed', 'Modelled'), shadow=True, loc='upper left')
330        savefig(name, dpi = 300)       
[2229]331
[3850]332        raw_input('Next')
[2229]333
334
[3850]335
Note: See TracBrowser for help on using the repository browser.