source: branches/numpy_anuga_validation/automated_validation_tests/okushiri_tank_validation/compare_timeseries_with_measures.py @ 6473

Last change on this file since 6473 was 6160, checked in by rwilson, 15 years ago

Change Numeric imports to general form - ready to change to NumPy?.

File size: 7.6 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 as png files.
5No plots are shown on screen.
6"""
7
8import Numeric as num
9from Scientific.IO.NetCDF import NetCDFFile
10
11from anuga.abstract_2d_finite_volumes.util import file_function
12from anuga.utilities.numerical_tools import\
13     ensure_numeric, cov, get_machine_precision
14
15import project
16import sys
17
18if sys.platform == 'win32':
19    # Windows has a problem when this module is run through
20    # os.system as done by validate_okushiri.
21    # See https://datamining.anu.edu.au/anuga/ticket/235
22
23    # If you want to see the plots from this validation,
24    # run this module by itself with this if clause removed.
25    plotting = False
26   
27else:
28    try:
29        from pylab import ion, hold, plot, title, legend
30        from pylab import xlabel, ylabel, savefig
31        hold(False)  # Check if this command can be issued
32    except:
33        print 'Could not import pylab'
34        plotting = False
35    else:
36        # Create plots as png files
37        plotting = True
38   
39
40#-------------------------
41# Basic data
42#-------------------------
43
44finaltime = 22.5
45timestep = 0.05
46
47gauge_locations = [[0.000, 1.696]] # Boundary gauge
48gauge_locations += [[4.521, 1.196],  [4.521, 1.696],  [4.521, 2.196]] #Ch 5-7-9
49gauge_names = ['Boundary', 'ch5', 'ch7', 'ch9']
50
51validation_data = {}
52for key in gauge_names:
53    validation_data[key] = []
54
55
56# Expected values
57expected_covariance = {'Boundary': 5.269569575007607815e-05, 
58                       'ch5': 1.166277999581819919e-04,
59                       'ch7': 1.127136457890861503e-04,
60                       'ch9': 1.250659477418482129e-04}
61
62expected_difference = {'Boundary': 8.350712673810733924e-04,
63                       'ch5': 3.405426180525532483e-03,
64                       'ch7': 2.852870417368218517e-03,
65                       'ch9': 3.248778982037564891e-03}
66
67expected_maximum = {'Boundary': 1.611749508386188523e-02, 
68                    'ch5': 3.551308418158714147e-02,
69                    'ch7': 3.858418457126511908e-02,
70                    'ch9': 4.317962986578308127e-02}
71
72expected_minimum = {'Boundary': -1.164547474575844919e-02, 
73                    'ch5': -8.664439185502026408e-03,
74                    'ch7': -2.726335488279797541e-03,
75                    'ch9': -5.977581218447349659e-03}
76
77expected_argmax = {'Boundary': 1.255000000000000071e+01, 
78                   'ch5': 1.839999999999999858e+01,
79                   'ch7': 1.700000000000000000e+01,
80                   'ch9': 1.685000000000000142e+01}
81
82expected_argmin = {'Boundary': 2.064999999999999858e+01, 
83                   'ch5': 1.459999999999999964e+01,
84                   'ch7': 1.230000000000000071e+01,
85                   'ch9': 1.315000000000000036e+01}
86
87#-------------------------
88# Read validation dataa
89#-------------------------
90
91print 'Reading', project.boundary_filename
92fid = NetCDFFile(project.boundary_filename, 'r')
93input_time = fid.variables['time'][:]
94validation_data['Boundary'] = fid.variables['stage'][:]
95
96reference_time = []
97fid = open(project.validation_filename)
98lines = fid.readlines()
99fid.close()
100
101for i, line in enumerate(lines[1:]):
102    if i == len(input_time): break
103   
104    fields = line.split()
105
106    reference_time.append(float(fields[0]))    # Record reference time
107    for j, key in enumerate(gauge_names[1:]):  # Omit boundary gauge
108        value = float(fields[1:][j])           # Omit time
109        validation_data[key].append(value/100) # Convert cm2m
110
111
112# Checks
113assert reference_time[0] == 0.0
114assert reference_time[-1] == finaltime
115assert num.allclose(reference_time, input_time)
116
117for key in gauge_names:
118    validation_data[key] = ensure_numeric(validation_data[key])
119
120#--------------------------------------------------
121# Read and interpolate model output
122#--------------------------------------------------
123
124import sys
125if len(sys.argv) > 1:
126    sww_filename = sys.argv[1]
127else:   
128    sww_filename = project.output_filename
129   
130f = file_function(sww_filename,
131                  quantities='stage',
132                  interpolation_points=gauge_locations,
133                  use_cache=True,
134                  verbose=True)
135
136
137def report_difference(name, computed_value, reference_value, rtol, atol):
138
139    if abs(reference_value) > 0:
140        msg = '%s (expected, computed):\n  (%.18e, %.18e):\n  Relative error=%.18e'\
141              %(name, reference_value, computed_value,
142                abs(reference_value-computed_value)/reference_value)
143        print msg
144       
145
146    msg = '  Absolute error=%.18e'\
147          %(abs(reference_value-computed_value))       
148    print msg
149
150   
151    #print 'Allclose:', allclose(reference_value, computed_value,
152    #                            rtol=rtol, atol=atol)
153    if plotting is False:
154        assert num.allclose(reference_value, computed_value,
155                            rtol=rtol, atol=atol), msg
156   
157
158
159#--------------------------------------------------
160# Compare model output to validation data
161#--------------------------------------------------
162
163
164#eps = get_machine_precision()
165
166# Windows tolerances
167rtol = 1.0e-2
168atol = 1.0e-3
169print 'Precisions used: rtol=%e, atol=%e' %(rtol, atol)
170
171for k, name in enumerate(gauge_names):
172
173    sqsum = 0
174    denom = 0
175    model = []
176    print 
177    print 'Validating ' + name
178    observed_timeseries = validation_data[name]
179    for i, t in enumerate(reference_time):
180        model.append(f(t, point_id=k)[0])
181
182    # Covariance measure   
183    res = cov(observed_timeseries, model)
184    report_difference('Covariance', res, expected_covariance[name], rtol, atol)
185     
186    # Difference measures   
187    res = sum(abs(observed_timeseries-model))/len(model)
188    report_difference('Accumulated difference', res,
189                      expected_difference[name], rtol, atol)   
190
191    # Extrema
192    res = max(model)
193    report_difference('Maximum', res, expected_maximum[name], rtol, atol)
194   
195    res = min(model)
196    report_difference('Minimum', res, expected_minimum[name], rtol, atol)   
197
198    # Locations of extrema
199    #i0 = argmax(observed_timeseries)
200    i1 = num.argmax(model)
201    res = reference_time[i1]
202    report_difference('Location of maximum', res, expected_argmax[name], rtol, atol)   
203   
204
205    if not name in ['ch7', 'ch9']:
206        # Minima of ch7 and ch9 are very flat and hard to pinpoint
207        i1 = num.argmin(model)
208        res = reference_time[i1]
209        report_difference('Location of minimum', res, expected_argmin[name],
210                          rtol, atol)       
211
212
213    if plotting is True:
214        #ion() # No plotting on screen
215        hold(False)
216   
217        plot(reference_time, validation_data[name], 'r-',
218             reference_time, model, 'k-')
219        title('Gauge %s' %name)
220        xlabel('time(s)')
221        ylabel('stage (m)')   
222        legend(('Observed', 'Modelled'), shadow=True, loc='upper left')
223        savefig(name, dpi = 300)       
224
225
226
227
228# Check max runup
229
230from anuga.shallow_water.data_manager import get_maximum_inundation_elevation
231from anuga.shallow_water.data_manager import get_maximum_inundation_location
232from anuga.utilities.polygon import is_inside_polygon
233
234q = get_maximum_inundation_elevation(sww_filename)
235loc = get_maximum_inundation_location(sww_filename)
236
237
238print 'Max runup elevation: ', q
239print 'Max runup elevation (scaled by 400): ', q*400
240print 'Max runup location:  ', loc
241
242
243from create_okushiri import gulleys
244assert is_inside_polygon(loc, gulleys)
245
246# FIXME more asserts here
247
248
249
250#msg = 'We got %f, should have been %f' %(q, q_max)
251#assert allclose(q, q_max, rtol=1.0/N), msg
252##print 'loc', loc, q
253#assert allclose(-loc[0]/2, q) # From topography formula
254
255print 'OK'
Note: See TracBrowser for help on using the repository browser.