source: anuga_validation/okushiri_2005/compare_timeseries.py @ 4227

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

Diagnostics

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