source: anuga_validation/automated_validation_tests/okushiri_tank_validation/compare_timeseries_with_measures.py @ 3917

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

More automated validation

File size: 13.8 KB
Line 
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.
5"""
6
7from Numeric import allclose, argmin, argmax
8from Scientific.IO.NetCDF import NetCDFFile
9
10from anuga.abstract_2d_finite_volumes.util import file_function
11from anuga.utilities.numerical_tools import\
12     ensure_numeric, cov, get_machine_precision
13
14import project
15
16try:
17    from pylab import ion, hold, plot, title, legend, xlabel, ylabel, savefig
18except:
19    plotting = False
20else:
21    plotting = True
22
23#plotting = False
24
25#-------------------------
26# Basic data
27#-------------------------
28
29finaltime = 22.5
30timestep = 0.05
31
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
34gauge_names = ['Boundary', 'ch5', 'ch7', 'ch9']
35
36validation_data = {}
37for key in gauge_names:
38    validation_data[key] = []
39
40
41# Results from lwru2_variable_mesh.sww
42#
43#Validating Boundary
44#Covariance = 5.265850426352561254e-05
45#Accumulated difference = 9.188561500629542633e-04
46#Difference in maxima = 4.865498851439054029e-05
47#Difference in minima = 2.293007809784416290e-04
48#Timelag between maxima = 3.000000000000007105e-01
49#Timelag between minima = 3.000000000000007105e-01
50#Next
51#
52#Validating ch5
53#Covariance = 1.171913883449135149e-04
54#Accumulated difference = 3.439740513194261117e-03
55#Difference in maxima = 1.808007760316858448e-03
56#Difference in minima = 1.494336068527772621e-04
57#Timelag between maxima = 9.999999999999786837e-02
58#Timelag between minima = 8.499999999999996447e-01
59#Next#
60#
61#Validating ch7
62#Covariance = 1.130834789656423632e-04
63#Accumulated difference = 2.868914369295409543e-03
64#Difference in maxima = 5.464055727564115506e-04
65#Difference in minima = 4.515664278503157131e-03
66#Timelag between maxima = 0.000000000000000000e+00
67#Timelag between minima = 4.500000000000010658e-01
68#
69#Next
70#Validating ch9
71#Covariance = 1.257648236914514148e-04
72#Accumulated difference = 3.287722754612562841e-03
73#Difference in maxima = 1.856547605566957748e-03
74#Difference in minima = 3.406785970601119290e-03
75#Timelag between maxima = 0.000000000000000000e+00
76#Timelag between minima = 9.499999999999992895e-01
77
78
79
80
81
82# New validation mesh
83#Validating Boundary
84#Covariance = 5.294418463841652839e-05
85#Accumulated difference = 8.222031767747334978e-04
86#Difference in maxima = 2.074639250243537347e-05
87#Difference in minima = 2.594698240608837858e-04
88#Timelag between maxima = 2.500000000000000000e-01
89#Timelag between minima = 2.500000000000000000e-01
90#Next
91#
92#Validating ch5
93#Covariance = 1.169700827457795071e-04
94#Accumulated difference = 3.390769873961110153e-03
95#Difference in maxima = 1.443103996215826246e-03
96#Difference in minima = 6.360676810760659133e-04
97#Timelag between maxima = 4.999999999999715783e-02
98#Timelag between minima = 7.999999999999989342e-01
99#Next
100#
101#Validating ch7
102#Covariance = 1.131522443178145860e-04
103#Accumulated difference = 2.815084054773397005e-03
104#Difference in maxima = 8.380973795966278894e-05
105#Difference in minima = 4.604184345996111850e-03
106#Timelag between maxima = 0.000000000000000000e+00
107#Timelag between minima = 2.050000000000000711e+00
108#
109#Next
110#Validating ch9
111#Covariance = 1.251932626774203753e-04
112#Accumulated difference = 3.208757435748204462e-03
113#Difference in maxima = 2.081659579598611753e-03
114#Difference in minima = 3.428854197019200710e-03
115#Timelag between maxima = 0.000000000000000000e+00
116#Timelag between minima = 1.000000000000000000e+00
117
118
119#expected_covariance = {'Boundary': 5.294418463841652839e-05,
120#                       'ch5': 1.169700827457795071e-04,
121#                       'ch7': 1.131522443178145860e-04,
122#                       'ch9': 1.251932626774203753e-04}
123#
124#expected_difference = {'Boundary': 8.222031767747334978e-04,
125#                       'ch5': 3.390769873961110153e-03,
126#                       'ch7': 2.815084054773397005e-03,
127#                       'ch9': 3.208757435748204462e-03}
128#
129#expected_maximum_diff = {'Boundary': 2.074639250243537347e-05,
130#                         'ch5': 1.443103996215826246e-03,
131#                         'ch7': 8.380973795966278894e-05,
132#                         'ch9': 2.081659579598611753e-03}
133#
134#expected_minimum_diff = {'Boundary': 2.594698240608837858e-04,
135#                         'ch5': 6.360676810760659133e-04,
136#                         'ch7': 4.604184345996111850e-03,
137#                         'ch9': 3.428854197019200710e-03}
138#
139#expected_argmax_timelag = {'Boundary': 2.500000000000000000e-01,
140#                           'ch5': 4.999999999999715783e-02,
141#                           'ch7': 0.000000000000000000e+00,
142#                           'ch9': 0.000000000000000000e+00}
143#
144#expected_argmin_timelag = {'Boundary': 2.500000000000000000e-01,
145#                           'ch5': 7.999999999999989342e-01,
146#                           'ch7': 2.050000000000000711e+00,
147#                           'ch9': 1.000000000000000000e+00}
148
149
150
151
152# New validation mesh - smaller regions
153#Validating Boundary
154#Covariance = 5.269749049737429159e-05
155#Accumulated difference = 8.337114775060755880e-04
156#Difference in maxima = 7.166795837686146253e-05
157#Difference in minima = 1.335627610594206094e-04
158#Timelag between maxima = 3.000000000000007105e-01
159#Timelag between minima = 2.500000000000000000e-01
160#Next
161#
162#Validating ch5
163#Covariance = 1.163749077455013324e-04
164#Accumulated difference = 3.385911115917152742e-03
165#Difference in maxima = 1.326867779957524585e-03
166#Difference in minima = 9.334970576141445042e-04
167#Timelag between maxima = 4.999999999999715783e-02
168#Timelag between minima = 7.999999999999989342e-01
169#Next
170#
171#Validating ch7
172#Covariance = 1.126985405540967749e-04
173#Accumulated difference = 2.857119198201266575e-03
174#Difference in maxima = 3.582696745840288632e-04
175#Difference in minima = 4.523664511720201960e-03
176#Timelag between maxima = 0.000000000000000000e+00
177#Timelag between minima = 2.150000000000000355e+00
178#Next
179#
180#Validating ch9
181#Covariance = 1.250280909578863607e-04
182#Accumulated difference = 3.247257369996859894e-03
183#Difference in maxima = 2.212919173501755321e-03
184#Difference in minima = 3.427581218447350343e-03
185#Timelag between maxima = 0.000000000000000000e+00
186#Timelag between minima = 1.150000000000000355e+00
187
188expected_covariance = {'Boundary': 5.269749049737429159e-05, 
189                       'ch5': 1.163749077455013324e-04,
190                       'ch7': 1.126985405540967749e-04,
191                       'ch9': 1.250280909578863607e-04}
192
193expected_difference = {'Boundary': 8.337114775060755880e-04,
194                       'ch5': 3.385911115917152742e-03,
195                       'ch7': 2.857119198201266575e-03,
196                       'ch9': 3.247257369996859894e-03}
197
198expected_maximum_diff = {'Boundary': 7.166795837686146253e-05,
199                         'ch5': 1.326867779957524585e-03,
200                         'ch7': 3.582696745840288632e-04,
201                         'ch9': 2.212919173501755321e-03}
202
203expected_minimum_diff = {'Boundary': 1.335627610594206094e-04,
204                         'ch5': 9.334970576141445042e-04,
205                         'ch7': 4.523664511720201960e-03,
206                         'ch9': 3.427581218447350343e-03}
207
208expected_argmax_timelag = {'Boundary': 3.000000000000007105e-01,
209                           'ch5': 4.999999999999715783e-02,
210                           'ch7': 0.000000000000000000e+00,
211                           'ch9': 0.000000000000000000e+00}
212
213expected_argmin_timelag = {'Boundary': 2.500000000000000000e-01,
214                           'ch5': 7.999999999999989342e-01,
215                           'ch7': 2.150000000000000355e+00,
216                           'ch9': 1.150000000000000355e+00}
217
218
219
220#-------------------------
221# Read validation dataa
222#-------------------------
223
224print 'Reading', project.boundary_filename
225fid = NetCDFFile(project.boundary_filename, 'r')
226input_time = fid.variables['time'][:]
227validation_data['Boundary'] = fid.variables['stage'][:]
228
229reference_time = []
230fid = open(project.validation_filename)
231lines = fid.readlines()
232fid.close()
233
234for i, line in enumerate(lines[1:]):
235    if i == len(input_time): break
236   
237    fields = line.split()
238
239    reference_time.append(float(fields[0]))    # Record reference time
240    for j, key in enumerate(gauge_names[1:]):  # Omit boundary gauge
241        value = float(fields[1:][j])           # Omit time
242        validation_data[key].append(value/100) # Convert cm2m
243
244
245# Checks
246assert reference_time[0] == 0.0
247assert reference_time[-1] == finaltime
248assert allclose(reference_time, input_time)
249
250for key in gauge_names:
251    validation_data[key] = ensure_numeric(validation_data[key])
252
253
254#--------------------------------------------------
255# Read and interpolate model output
256#--------------------------------------------------
257
258import sys
259if len(sys.argv) > 1:
260    sww_filename = sys.argv[1]
261else:   
262    sww_filename = project.output_filename
263   
264f = file_function(sww_filename,
265                  quantities='stage',
266                  interpolation_points=gauge_locations,
267                  use_cache=True,
268                  verbose=True)
269
270
271#--------------------------------------------------
272# Compare model output to validation data
273#--------------------------------------------------
274
275eps = get_machine_precision()
276for k, name in enumerate(gauge_names):
277    sqsum = 0
278    denom = 0
279    model = []
280    print 
281    print 'Validating ' + name
282    observed_timeseries = validation_data[name]
283    for i, t in enumerate(reference_time):
284        model.append(f(t, point_id=k)[0])
285
286
287    # Covariance measures   
288    res = cov(observed_timeseries, model)     
289    print 'Covariance = %.18e' %res
290
291    if res < expected_covariance[name]-eps:
292        print '  Result is better than expected by: %.18e'\
293              %(res-expected_covariance[name])
294        print '  Expect = %.18e' %expected_covariance[name]   
295    elif res > expected_covariance[name]+eps:
296        print '  FAIL: Result is worse than expected by: %.18e'\
297                                %(res-expected_covariance[name])
298        print '  Expect = %.18e' %expected_covariance[name]
299    else:
300        pass
301
302    assert res-eps < expected_covariance[name] < res+eps
303   
304   
305
306    # Difference measures   
307    res = sum(abs(observed_timeseries-model))/len(model)     
308    print 'Accumulated difference = %.18e' %res
309
310    if res < expected_difference[name]-eps:
311        print '  Result is better than expected by: %.18e'\
312              %(res-expected_difference[name])
313        print '  Expect = %.18e' %expected_difference[name]   
314    elif res > expected_difference[name]+eps:
315        print '  FAIL: Result is worse than expected by: %.18e'\
316                                %(res-expected_difference[name])
317        print '  Expect = %.18e' %expected_difference[name]
318    else:
319        pass
320
321
322    assert res-eps < expected_difference[name] < res+eps
323
324    # Extrema
325    res = abs(max(observed_timeseries)-max(model))
326    print 'Difference in maxima = %.18e' %res
327
328    if res < expected_maximum_diff[name]-eps:
329        print '  Result is better than expected by: %.18e'\
330              %(res-expected_maximum_diff[name])
331        print '  Expect = %.18e' %expected_maximum_diff[name]   
332    elif res > expected_maximum_diff[name]+eps:
333        print '  FAIL: Result is worse than expected by: %.18e'\
334                                %(res-expected_maximum_diff[name])
335        print '  Expect = %.18e' %expected_maximum_diff[name]
336    else:
337        pass
338
339    assert res-eps < expected_maximum_diff[name] < res+eps   
340
341   
342    res = abs(min(observed_timeseries)-min(model))
343    print 'Difference in minima = %.18e' %res
344
345    if res < expected_minimum_diff[name]-eps:
346        print '  Result is better than expected by: %.18e'\
347              %(res-expected_minimum_diff[name])
348        print '  Expect = %.18e' %expected_minimum_diff[name]   
349    elif res > expected_minimum_diff[name]+eps:
350        print '  FAIL: Result is worse than expected by: %.18e'\
351                                %(res-expected_minimum_diff[name])
352        print '  Expect = %.18e' %expected_minimum_diff[name]
353    else:
354        pass
355
356    assert res-eps < expected_minimum_diff[name] < res+eps
357   
358
359    # Locations of extrema
360    i0 = argmax(observed_timeseries)
361    i1 = argmax(model)
362    res = abs(reference_time[i1] - reference_time[i0])
363    print 'Timelag between maxima = %.18e' %res
364
365    if res < expected_argmax_timelag[name]-eps:
366        print '  Result is better than expected by: %.18e'\
367              %(res-expected_argmax_timelag[name])
368        print '  Expect = %.18e' %expected_argmax_timelag[name]   
369    elif res > expected_argmax_timelag[name]+eps:
370        print '  FAIL: Result is worse than expected by: %.18e'\
371                                %(res-expected_argmax_timelag[name])
372        print '  Expect = %.18e' %expected_argmax_timelag[name]
373    else:
374        pass
375   
376    assert res-eps < expected_argmax_timelag[name] < res+eps
377
378    i0 = argmin(observed_timeseries)
379    i1 = argmin(model)
380    res = abs(reference_time[i1] - reference_time[i0])
381    print 'Timelag between minima = %.18e' %res
382    if res < expected_argmin_timelag[name]-eps:
383        print '  Result is better than expected by: %.18e'\
384              %(res-expected_argmin_timelag[name])
385        print '  Expect = %.18e' %expected_argmin_timelag[name]   
386    elif res > expected_argmin_timelag[name]+eps:
387        print '  FAIL: Result is worse than expected by: %.18e'\
388                                %(res-expected_argmin_timelag[name])
389        print '  Expect = %.18e' %expected_argmin_timelag[name]
390    else:
391        pass
392
393    assert res-eps < expected_argmin_timelag[name] < res+eps
394
395
396    if plotting is True:
397        ion()
398        hold(False)
399   
400        plot(reference_time, validation_data[name], 'r-',
401             reference_time, model, 'k-')
402        title('Gauge %s' %name)
403        xlabel('time(s)')
404        ylabel('stage (m)')   
405        legend(('Observed', 'Modelled'), shadow=True, loc='upper left')
406        savefig(name, dpi = 300)       
407
408        raw_input('Next')
409
410
411
Note: See TracBrowser for help on using the repository browser.