1 | """Verify that simulation produced by ANUGA compares to published |
---|
2 | validation timeseries ch5, ch7 and ch9 as well as the boundary timeseries. |
---|
3 | |
---|
4 | RMS norm is printed and plots are produced. |
---|
5 | """ |
---|
6 | |
---|
7 | from Numeric import allclose, argmin, argmax |
---|
8 | from Scientific.IO.NetCDF import NetCDFFile |
---|
9 | |
---|
10 | from anuga.abstract_2d_finite_volumes.util import file_function |
---|
11 | from anuga.utilities.numerical_tools import\ |
---|
12 | ensure_numeric, cov, get_machine_precision |
---|
13 | |
---|
14 | import project |
---|
15 | |
---|
16 | try: |
---|
17 | from pylab import ion, hold, plot, title, legend, xlabel, ylabel, savefig |
---|
18 | except: |
---|
19 | plotting = False |
---|
20 | else: |
---|
21 | plotting = True |
---|
22 | |
---|
23 | #plotting = False |
---|
24 | |
---|
25 | #------------------------- |
---|
26 | # Basic data |
---|
27 | #------------------------- |
---|
28 | |
---|
29 | finaltime = 22.5 |
---|
30 | timestep = 0.05 |
---|
31 | |
---|
32 | gauge_locations = [[0.000, 1.696]] # Boundary gauge |
---|
33 | gauge_locations += [[4.521, 1.196], [4.521, 1.696], [4.521, 2.196]] #Ch 5-7-9 |
---|
34 | gauge_names = ['Boundary', 'ch5', 'ch7', 'ch9'] |
---|
35 | |
---|
36 | validation_data = {} |
---|
37 | for key in gauge_names: |
---|
38 | validation_data[key] = [] |
---|
39 | |
---|
40 | |
---|
41 | #expected_covariances = {'Boundary': 5.288392008865989e-05, |
---|
42 | # 'ch5': 1.166748190444681e-04, |
---|
43 | # 'ch7': 1.121816242516758e-04, |
---|
44 | # 'ch9': 1.249543278366778e-04} |
---|
45 | |
---|
46 | # old limiters |
---|
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} |
---|
56 | |
---|
57 | |
---|
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} |
---|
67 | |
---|
68 | |
---|
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 | # |
---|
79 | |
---|
80 | |
---|
81 | expected_covariance = {'Boundary': 5.265850426352561254e-05, |
---|
82 | 'ch5': 1.171913883449135149e-04, |
---|
83 | 'ch7': 1.130834789656423632e-04, |
---|
84 | 'ch9': 1.257648236914514148e-04} |
---|
85 | |
---|
86 | expected_difference = {'Boundary': 9.188561500629542633e-04, |
---|
87 | 'ch5': 3.439740513194261117e-03, |
---|
88 | 'ch7': 2.868914369295409543e-03, |
---|
89 | 'ch9': 3.287722754612562841e-03} |
---|
90 | |
---|
91 | expected_maximum_diff = {'Boundary': 4.865498851439054029e-05, |
---|
92 | 'ch5': 1.808007760316858448e-03, |
---|
93 | 'ch7': 5.464055727564115506e-04, |
---|
94 | 'ch9': 1.856547605566957748e-03} |
---|
95 | |
---|
96 | expected_minimum_diff = {'Boundary': 2.293007809784416290e-04, |
---|
97 | 'ch5': 1.494336068527772621e-04, |
---|
98 | 'ch7': 4.515664278503157131e-03, |
---|
99 | 'ch9': 3.406785970601119290e-03} |
---|
100 | |
---|
101 | expected_argmax_timelag = {'Boundary': 3.000000000000007105e-01, |
---|
102 | 'ch5': 9.999999999999786837e-02, |
---|
103 | 'ch7': 0.000000000000000000e+00, |
---|
104 | 'ch9': 0.000000000000000000e+00} |
---|
105 | |
---|
106 | expected_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 | |
---|
151 | #------------------------- |
---|
152 | # Read validation dataa |
---|
153 | #------------------------- |
---|
154 | |
---|
155 | print 'Reading', project.boundary_filename |
---|
156 | fid = NetCDFFile(project.boundary_filename, 'r') |
---|
157 | input_time = fid.variables['time'][:] |
---|
158 | validation_data['Boundary'] = fid.variables['stage'][:] |
---|
159 | |
---|
160 | reference_time = [] |
---|
161 | fid = open(project.validation_filename) |
---|
162 | lines = fid.readlines() |
---|
163 | fid.close() |
---|
164 | |
---|
165 | for i, line in enumerate(lines[1:]): |
---|
166 | if i == len(input_time): break |
---|
167 | |
---|
168 | fields = line.split() |
---|
169 | |
---|
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 |
---|
174 | |
---|
175 | |
---|
176 | # Checks |
---|
177 | assert reference_time[0] == 0.0 |
---|
178 | assert reference_time[-1] == finaltime |
---|
179 | assert allclose(reference_time, input_time) |
---|
180 | |
---|
181 | for key in gauge_names: |
---|
182 | validation_data[key] = ensure_numeric(validation_data[key]) |
---|
183 | |
---|
184 | |
---|
185 | #-------------------------------------------------- |
---|
186 | # Read and interpolate model output |
---|
187 | #-------------------------------------------------- |
---|
188 | |
---|
189 | import sys |
---|
190 | if len(sys.argv) > 1: |
---|
191 | sww_filename = sys.argv[1] |
---|
192 | else: |
---|
193 | sww_filename = project.output_filename |
---|
194 | |
---|
195 | #f = file_function('okushiri_new_limiters.sww', #The best so far |
---|
196 | #f = file_function('okushiri_as2005_with_mxspd=0.1.sww', |
---|
197 | f = file_function(sww_filename, |
---|
198 | quantities='stage', |
---|
199 | interpolation_points=gauge_locations, |
---|
200 | use_cache=True, |
---|
201 | verbose=True) |
---|
202 | |
---|
203 | |
---|
204 | #-------------------------------------------------- |
---|
205 | # Compare model output to validation data |
---|
206 | #-------------------------------------------------- |
---|
207 | |
---|
208 | eps = get_machine_precision() |
---|
209 | for k, name in enumerate(gauge_names): |
---|
210 | sqsum = 0 |
---|
211 | denom = 0 |
---|
212 | model = [] |
---|
213 | print |
---|
214 | print 'Validating ' + name |
---|
215 | observed_timeseries = validation_data[name] |
---|
216 | for i, t in enumerate(reference_time): |
---|
217 | model.append(f(t, point_id=k)[0]) |
---|
218 | |
---|
219 | |
---|
220 | # Covariance measures |
---|
221 | res = cov(observed_timeseries, model) |
---|
222 | print 'Covariance = %.18e' %res |
---|
223 | |
---|
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] |
---|
232 | else: |
---|
233 | pass |
---|
234 | |
---|
235 | |
---|
236 | # Difference measures |
---|
237 | res = sum(abs(observed_timeseries-model))/len(model) |
---|
238 | print 'Accumulated difference = %.18e' %res |
---|
239 | |
---|
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] |
---|
248 | else: |
---|
249 | pass |
---|
250 | |
---|
251 | |
---|
252 | # Extrema |
---|
253 | res = abs(max(observed_timeseries)-max(model)) |
---|
254 | print 'Difference in maxima = %.18e' %res |
---|
255 | |
---|
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 | |
---|
267 | |
---|
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 | |
---|
320 | if plotting is True: |
---|
321 | ion() |
---|
322 | hold(False) |
---|
323 | |
---|
324 | plot(reference_time, validation_data[name], 'r-', |
---|
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) |
---|
331 | |
---|
332 | raw_input('Next') |
---|
333 | |
---|
334 | |
---|
335 | |
---|