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