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

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

Removed raw_input to allow code to run in the background

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