[3845] | 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. |
---|
[2229] | 5 | """ |
---|
| 6 | |
---|
[3913] | 7 | from Numeric import allclose, argmin, argmax |
---|
[3850] | 8 | from Scientific.IO.NetCDF import NetCDFFile |
---|
[3845] | 9 | |
---|
[3850] | 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 |
---|
[3845] | 13 | |
---|
[3850] | 14 | import project |
---|
[2229] | 15 | |
---|
[3850] | 16 | try: |
---|
| 17 | from pylab import ion, hold, plot, title, legend, xlabel, ylabel, savefig |
---|
| 18 | except: |
---|
| 19 | plotting = False |
---|
| 20 | else: |
---|
| 21 | plotting = True |
---|
[2229] | 22 | |
---|
[3861] | 23 | #plotting = False |
---|
[2229] | 24 | |
---|
[3850] | 25 | #------------------------- |
---|
| 26 | # Basic data |
---|
| 27 | #------------------------- |
---|
[3845] | 28 | |
---|
[3850] | 29 | finaltime = 22.5 |
---|
| 30 | timestep = 0.05 |
---|
[2229] | 31 | |
---|
[3850] | 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 |
---|
[2229] | 34 | gauge_names = ['Boundary', 'ch5', 'ch7', 'ch9'] |
---|
| 35 | |
---|
[3850] | 36 | validation_data = {} |
---|
| 37 | for 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] | 81 | expected_covariance = {'Boundary': 5.265850426352561254e-05, |
---|
| 82 | 'ch5': 1.171913883449135149e-04, |
---|
| 83 | 'ch7': 1.130834789656423632e-04, |
---|
| 84 | 'ch9': 1.257648236914514148e-04} |
---|
[3883] | 85 | |
---|
[3913] | 86 | expected_difference = {'Boundary': 9.188561500629542633e-04, |
---|
| 87 | 'ch5': 3.439740513194261117e-03, |
---|
| 88 | 'ch7': 2.868914369295409543e-03, |
---|
| 89 | 'ch9': 3.287722754612562841e-03} |
---|
[3883] | 90 | |
---|
[3913] | 91 | expected_maximum_diff = {'Boundary': 4.865498851439054029e-05, |
---|
| 92 | 'ch5': 1.808007760316858448e-03, |
---|
| 93 | 'ch7': 5.464055727564115506e-04, |
---|
| 94 | 'ch9': 1.856547605566957748e-03} |
---|
[3883] | 95 | |
---|
[3913] | 96 | expected_minimum_diff = {'Boundary': 2.293007809784416290e-04, |
---|
| 97 | 'ch5': 1.494336068527772621e-04, |
---|
| 98 | 'ch7': 4.515664278503157131e-03, |
---|
| 99 | 'ch9': 3.406785970601119290e-03} |
---|
[3883] | 100 | |
---|
[3913] | 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 | |
---|
[3850] | 151 | #------------------------- |
---|
| 152 | # Read validation dataa |
---|
| 153 | #------------------------- |
---|
[2229] | 154 | |
---|
[3850] | 155 | print 'Reading', project.boundary_filename |
---|
| 156 | fid = NetCDFFile(project.boundary_filename, 'r') |
---|
[2229] | 157 | input_time = fid.variables['time'][:] |
---|
[3850] | 158 | validation_data['Boundary'] = fid.variables['stage'][:] |
---|
[2229] | 159 | |
---|
| 160 | reference_time = [] |
---|
[3850] | 161 | fid = open(project.validation_filename) |
---|
[2229] | 162 | lines = fid.readlines() |
---|
| 163 | fid.close() |
---|
[3850] | 164 | |
---|
[2229] | 165 | for 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] | 177 | assert reference_time[0] == 0.0 |
---|
| 178 | assert reference_time[-1] == finaltime |
---|
| 179 | assert allclose(reference_time, input_time) |
---|
| 180 | |
---|
[3850] | 181 | for 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 | |
---|
| 189 | import sys |
---|
| 190 | if len(sys.argv) > 1: |
---|
| 191 | sww_filename = sys.argv[1] |
---|
| 192 | else: |
---|
| 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] | 197 | f = 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 | |
---|
| 208 | eps = get_machine_precision() |
---|
[2229] | 209 | for 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 | |
---|