source: anuga_validation/okushiri_2005/compare_timeseries.py @ 3915

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

Okushiri demo

File size: 3.9 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
7from Numeric import allclose, argmin, argmax
8from Scientific.IO.NetCDF import NetCDFFile
9
10from anuga.abstract_2d_finite_volumes.util import file_function
11from anuga.utilities.numerical_tools import\
12     ensure_numeric, cov, get_machine_precision
13
14import project
15
16try:
17    from pylab import ion, hold, plot, title, legend, xlabel, ylabel, savefig
18except:
19    plotting = False
20else:
21    plotting = True
22
23#plotting = False
24
25#-------------------------
26# Basic data
27#-------------------------
28
29finaltime = 22.5
30timestep = 0.05
31
32gauge_locations = [[0.000, 1.696]] # Boundary gauge
33gauge_locations += [[4.521, 1.196],  [4.521, 1.696],  [4.521, 2.196]] #Ch 5-7-9
34gauge_names = ['Boundary', 'ch5', 'ch7', 'ch9']
35
36validation_data = {}
37for key in gauge_names:
38    validation_data[key] = []
39
40
41#-------------------------
42# Read validation dataa
43#-------------------------
44
45print 'Reading', project.boundary_filename
46fid = NetCDFFile(project.boundary_filename, 'r')
47input_time = fid.variables['time'][:]
48validation_data['Boundary'] = fid.variables['stage'][:]
49
50reference_time = []
51fid = open(project.validation_filename)
52lines = fid.readlines()
53fid.close()
54
55for i, line in enumerate(lines[1:]):
56    if i == len(input_time): break
57   
58    fields = line.split()
59
60    reference_time.append(float(fields[0]))    # Record reference time
61    for j, key in enumerate(gauge_names[1:]):  # Omit boundary gauge
62        value = float(fields[1:][j])           # Omit time
63        validation_data[key].append(value/100) # Convert cm2m
64
65
66# Checks
67assert reference_time[0] == 0.0
68assert reference_time[-1] == finaltime
69assert allclose(reference_time, input_time)
70
71for key in gauge_names:
72    validation_data[key] = ensure_numeric(validation_data[key])
73
74
75#--------------------------------------------------
76# Read and interpolate model output
77#--------------------------------------------------
78
79import sys
80if len(sys.argv) > 1:
81    sww_filename = sys.argv[1]
82else:   
83    sww_filename = project.output_filename
84   
85f = file_function(sww_filename,
86                  quantities='stage',
87                  interpolation_points=gauge_locations,
88                  use_cache=True,
89                  verbose=True)
90
91
92#--------------------------------------------------
93# Compare model output to validation data
94#--------------------------------------------------
95
96eps = get_machine_precision()
97for k, name in enumerate(gauge_names):
98    sqsum = 0
99    denom = 0
100    model = []
101    print 
102    print 'Validating ' + name
103    observed_timeseries = validation_data[name]
104    for i, t in enumerate(reference_time):
105        model.append(f(t, point_id=k)[0])
106
107
108    # Covariance measures   
109    res = cov(observed_timeseries, model)     
110    print 'Covariance = %.18e' %res
111
112    # Difference measures   
113    res = sum(abs(observed_timeseries-model))/len(model)     
114    print 'Accumulated difference = %.18e' %res
115
116
117    # Extrema
118    res = abs(max(observed_timeseries)-max(model))
119    print 'Difference in maxima = %.18e' %res
120
121   
122    res = abs(min(observed_timeseries)-min(model))
123    print 'Difference in minima = %.18e' %res
124
125    # Locations of extrema
126    i0 = argmax(observed_timeseries)
127    i1 = argmax(model)
128    res = abs(reference_time[i1] - reference_time[i0])
129    print 'Timelag between maxima = %.18e' %res
130   
131
132    i0 = argmin(observed_timeseries)
133    i1 = argmin(model)
134    res = abs(reference_time[i1] - reference_time[i0])
135    print 'Timelag between minima = %.18e' %res
136
137
138
139
140    if plotting is True:
141        ion()
142        hold(False)
143   
144        plot(reference_time, validation_data[name], 'r-',
145             reference_time, model, 'k-')
146        title('Gauge %s' %name)
147        xlabel('time(s)')
148        ylabel('stage (m)')   
149        legend(('Observed', 'Modelled'), shadow=True, loc='upper left')
150        savefig(name, dpi = 300)       
151
152        raw_input('Next')
153
154
155
Note: See TracBrowser for help on using the repository browser.