source: anuga_validation/okushiri_2005/extract_timeseries.py @ 3637

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

Work on Okushiri validation (towards recreating 18 Aug 2005 results)

File size: 2.7 KB
Line 
1"""Read in tms file, interpolate at specified locations (ch 5,7,9) and compare
2"""
3
4import sys
5from os import sep
6from anuga.caching import cache
7
8sys.path.append('..'+sep+'..'+sep)
9
10
11gauges = [[0.000, 1.696]]  #Boundary
12gauges += [[4.521, 1.196],  [4.521, 1.696],  [4.521, 2.196]] #Ch 5-7-9
13
14gauge_names = ['Boundary', 'ch5', 'ch7', 'ch9']
15
16
17finaltime = 22.5
18timestep = 0.05
19
20#Read reference data
21
22#Input wave
23filename = 'Benchmark_2_input.tms'
24print 'Reading', filename
25from Scientific.IO.NetCDF import NetCDFFile
26from Numeric import allclose
27fid = NetCDFFile(filename, 'r')#
28
29input_time = fid.variables['time'][:]
30input_stage = fid.variables['stage'][:]
31
32
33#gauges
34reference_time = []
35ch5 = []
36ch7 = []
37ch9 = []
38filename = 'output_ch5-7-9.txt'
39fid = open(filename)
40lines = fid.readlines()
41fid.close()
42for 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
53from anuga.abstract_2d_finite_volumes.util import file_function
54from anuga.utilities.numerical_tools import ensure_numeric
55gauge_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'
64filename = 'lwru2.sww'
65#filename = 'lwru2_.05s.sww'
66
67#f = file_function(filename,
68#                  quantities = 'stage',
69#                  interpolation_points = gauges,
70#                  verbose = True)
71
72f = 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
86assert reference_time[0] == 0.0
87assert reference_time[-1] == finaltime
88assert allclose(reference_time, input_time)
89
90
91
92#Validation
93
94
95for 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')
123show()
124
125
126
127#from pylab import *
128#plot(time, stage)
129#show()
Note: See TracBrowser for help on using the repository browser.