[2229] | 1 | """Read in tms file, interpolate at specified locations (ch 5,7,9) and compare |
---|
| 2 | """ |
---|
| 3 | |
---|
| 4 | import sys |
---|
| 5 | from os import sep |
---|
| 6 | from caching import cache |
---|
| 7 | |
---|
| 8 | sys.path.append('..'+sep+'..'+sep) |
---|
| 9 | |
---|
| 10 | |
---|
| 11 | gauges = [[0.000, 1.696]] #Boundary |
---|
| 12 | gauges += [[4.521, 1.196], [4.521, 1.696], [4.521, 2.196]] #Ch 5-7-9 |
---|
| 13 | |
---|
| 14 | gauge_names = ['Boundary', 'ch5', 'ch7', 'ch9'] |
---|
| 15 | |
---|
| 16 | |
---|
| 17 | finaltime = 22.5 |
---|
| 18 | timestep = 0.05 |
---|
| 19 | |
---|
| 20 | #Read reference data |
---|
| 21 | |
---|
| 22 | #Input wave |
---|
| 23 | filename = 'Benchmark_2_input.tms' |
---|
| 24 | print 'Reading', filename |
---|
| 25 | from Scientific.IO.NetCDF import NetCDFFile |
---|
| 26 | from Numeric import allclose |
---|
| 27 | fid = NetCDFFile(filename, 'r')# |
---|
| 28 | |
---|
| 29 | input_time = fid.variables['time'][:] |
---|
| 30 | input_stage = fid.variables['stage'][:] |
---|
| 31 | |
---|
| 32 | |
---|
| 33 | #gauges |
---|
| 34 | reference_time = [] |
---|
| 35 | ch5 = [] |
---|
| 36 | ch7 = [] |
---|
| 37 | ch9 = [] |
---|
| 38 | filename = 'output_ch5-7-9.txt' |
---|
| 39 | fid = open(filename) |
---|
| 40 | lines = fid.readlines() |
---|
| 41 | fid.close() |
---|
| 42 | for i, line in enumerate(lines[1:]): |
---|
| 43 | if i == len(input_time): break |
---|
| 44 | |
---|
| 45 | fields = line.split() |
---|
| 46 | |
---|
| 47 | reference_time.append(float(fields[0])) |
---|
| 48 | ch5.append(float(fields[1])/100) #cm2m |
---|
| 49 | ch7.append(float(fields[2])/100) #cm2m |
---|
| 50 | ch9.append(float(fields[3])/100) #cm2m |
---|
| 51 | |
---|
| 52 | |
---|
[2609] | 53 | from pyvolution.util import file_function |
---|
| 54 | from utilities.numerical_tools import ensure_numeric |
---|
[2229] | 55 | gauge_values = [ensure_numeric(input_stage), |
---|
| 56 | ensure_numeric(ch5), |
---|
| 57 | ensure_numeric(ch7), |
---|
| 58 | ensure_numeric(ch9)] #Reference values |
---|
| 59 | |
---|
| 60 | |
---|
| 61 | |
---|
| 62 | #Read model output |
---|
| 63 | #filename = 'output.sww' |
---|
| 64 | filename = 'lwru2.sww' |
---|
| 65 | #filename = 'lwru2_.05s.sww' |
---|
| 66 | |
---|
| 67 | #f = file_function(filename, |
---|
| 68 | # quantities = 'stage', |
---|
| 69 | # interpolation_points = gauges, |
---|
| 70 | # verbose = True) |
---|
| 71 | |
---|
| 72 | f = cache(file_function,filename, |
---|
| 73 | {'quantities': 'stage', |
---|
| 74 | 'interpolation_points': gauges, |
---|
| 75 | 'verbose': True}, |
---|
| 76 | #evaluate = True, |
---|
| 77 | dependencies = [filename], |
---|
| 78 | verbose = True) |
---|
| 79 | |
---|
| 80 | |
---|
| 81 | |
---|
| 82 | |
---|
| 83 | #Checks |
---|
| 84 | #print reference_time |
---|
| 85 | #print input_time |
---|
| 86 | assert reference_time[0] == 0.0 |
---|
| 87 | assert reference_time[-1] == finaltime |
---|
| 88 | assert allclose(reference_time, input_time) |
---|
| 89 | |
---|
| 90 | |
---|
| 91 | |
---|
| 92 | #Validation |
---|
| 93 | |
---|
| 94 | |
---|
| 95 | for k, name in enumerate(gauge_names): |
---|
| 96 | sqsum = 0 |
---|
| 97 | denom = 0 |
---|
| 98 | model = [] |
---|
| 99 | print 'Validating ' + name |
---|
| 100 | for i, t in enumerate(reference_time): |
---|
| 101 | ref = gauge_values[k][i] |
---|
| 102 | val = f(t, point_id = k)[0] |
---|
| 103 | model.append(val) |
---|
| 104 | |
---|
| 105 | sqsum += (ref - val)**2 |
---|
| 106 | denom += ref**2 |
---|
| 107 | |
---|
| 108 | print sqsum |
---|
| 109 | print sqsum/denom |
---|
| 110 | |
---|
| 111 | from pylab import * |
---|
| 112 | ion() |
---|
| 113 | hold(False) |
---|
| 114 | plot(reference_time, gauge_values[k], 'r.-', |
---|
| 115 | reference_time, model, 'k-') |
---|
| 116 | title('Gauge %s' %name) |
---|
| 117 | xlabel('time(s)') |
---|
| 118 | ylabel('stage (m)') |
---|
| 119 | legend(('Observed', 'Modelled'), shadow=True, loc='upper left') |
---|
| 120 | savefig(name, dpi = 300) |
---|
| 121 | |
---|
| 122 | raw_input('Next') |
---|
| 123 | show() |
---|
| 124 | |
---|
| 125 | |
---|
| 126 | |
---|
| 127 | #from pylab import * |
---|
| 128 | #plot(time, stage) |
---|
| 129 | #show() |
---|