source: branches/numpy_anuga_validation/okushiri_2005/compare_timeseries.py @ 7184

Last change on this file since 7184 was 7184, checked in by rwilson, 15 years ago

Updated all licence files to use new checksum.

File size: 4.6 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
7import numpy as num
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
58    import numpy as num== len(input_time): break
59   
60    fields = line.split()
61
62    reference_time.append(float(fields[0]))    # Record reference time
63    for j, key in enumerate(gauge_names[1:]):  # Omit boundary gauge
64        value = float(fields[1:][j])           # Omit time
65        validation_data[key].append(value/100) # Convert cm2m
66
67
68# Checks
69assert reference_time[0] == 0.0
70assert reference_time[-1] == finaltime
71assert num.allclose(reference_time, input_time)
72
73for key in gauge_names:
74    validation_data[key] = ensure_numeric(validation_data[key])
75
76
77#--------------------------------------------------
78# Read and interpolate model output
79#--------------------------------------------------
80
81import sys
82if len(sys.argv) > 1:
83    sww_filename = sys.argv[1]
84else:   
85    sww_filename = project.output_filename
86   
87f = file_function(sww_filename,
88                  quantities='stage',
89                  interpolation_points=gauge_locations,
90                  use_cache=True,
91                  verbose=True)
92
93#--------------------------------------------------
94# Check max runup
95#--------------------------------------------------
96
97from anuga.shallow_water.data_manager import get_maximum_inundation_elevation
98from anuga.shallow_water.data_manager import get_maximum_inundation_location
99from anuga.utilities.polygon import is_inside_polygon
100
101q = get_maximum_inundation_elevation(sww_filename)
102loc = get_maximum_inundation_location(sww_filename)
103
104print 'Max runup elevation: ', q
105print 'Max runup elevation (scaled by 400): ', q*400
106print 'Max runup location:  ', loc
107
108#from create_okushiri import gulleys
109#assert is_inside_polygon(loc, gulleys)
110
111
112
113#--------------------------------------------------
114# Compare model output to validation data
115#--------------------------------------------------
116
117eps = get_machine_precision()
118for k, name in enumerate(gauge_names):
119    sqsum = 0
120    denom = 0
121    model = []
122    print 
123    print 'Validating ' + name
124    observed_timeseries = validation_data[name]
125    for i, t in enumerate(reference_time):
126        model.append(f(t, point_id=k)[0])
127
128
129    # Covariance measures   
130    res = cov(observed_timeseries, model)     
131    print 'Covariance = %.18e' %res
132
133    # Difference measures   
134    res = sum(abs(observed_timeseries-model))/len(model)     
135    print 'Accumulated difference = %.18e' %res
136
137
138    # Extrema
139    res = abs(max(observed_timeseries)-max(model))
140    print 'Difference in maxima = %.18e' %res
141
142   
143    res = abs(min(observed_timeseries)-min(model))
144    print 'Difference in minima = %.18e' %res
145
146    # Locations of extrema
147    i0 = num.argmax(observed_timeseries)
148    i1 = num.argmax(model)
149    res = abs(reference_time[i1] - reference_time[i0])
150    print 'Timelag between maxima = %.18e' %res
151   
152
153    i0 = num.argmin(observed_timeseries)
154    i1 = num.argmin(model)
155    res = abs(reference_time[i1] - reference_time[i0])
156    print 'Timelag between minima = %.18e' %res
157
158
159
160
161    if plotting is True:
162        ion()
163        hold(False)
164   
165        plot(reference_time, validation_data[name], 'r-',
166             reference_time, model, 'k-')
167        title('Gauge %s' %name)
168        xlabel('time(s)')
169        ylabel('stage (m)')   
170        legend(('Observed', 'Modelled'), shadow=True, loc='upper left')
171        savefig(name, dpi = 300)       
172
173        raw_input('Next')
174
175
176
Note: See TracBrowser for help on using the repository browser.