source: anuga_validation/okushiri_2005/extract_timeseries.py @ 3845

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

Work on Okushiri

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