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

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

Back-merge from Numeric trunk.

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 == 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 num.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# Check max runup
94#--------------------------------------------------
95
96from anuga.shallow_water.data_manager import get_maximum_inundation_elevation
97from anuga.shallow_water.data_manager import get_maximum_inundation_location
98from anuga.utilities.polygon import is_inside_polygon
99
100q = get_maximum_inundation_elevation(sww_filename)
101loc = get_maximum_inundation_location(sww_filename)
102
103print 'Max runup elevation: ', q
104print 'Max runup elevation (scaled by 400): ', q*400
105print 'Max runup location:  ', loc
106
107#from create_okushiri import gulleys
108#assert is_inside_polygon(loc, gulleys)
109
110
111
112#--------------------------------------------------
113# Compare model output to validation data
114#--------------------------------------------------
115
116eps = get_machine_precision()
117for k, name in enumerate(gauge_names):
118    sqsum = 0
119    denom = 0
120    model = []
121    print 
122    print 'Validating ' + name
123    observed_timeseries = validation_data[name]
124    for i, t in enumerate(reference_time):
125        model.append(f(t, point_id=k)[0])
126
127
128    # Covariance measures   
129    res = cov(observed_timeseries, model)     
130    print 'Covariance = %.18e' %res
131
132    # Difference measures   
133    res = sum(abs(observed_timeseries-model))/len(model)     
134    print 'Accumulated difference = %.18e' %res
135
136
137    # Extrema
138    res = abs(max(observed_timeseries)-max(model))
139    print 'Difference in maxima = %.18e' %res
140
141   
142    res = abs(min(observed_timeseries)-min(model))
143    print 'Difference in minima = %.18e' %res
144
145    # Locations of extrema
146    i0 = num.argmax(observed_timeseries)
147    i1 = num.argmax(model)
148    res = abs(reference_time[i1] - reference_time[i0])
149    print 'Timelag between maxima = %.18e' %res
150   
151
152    i0 = num.argmin(observed_timeseries)
153    i1 = num.argmin(model)
154    res = abs(reference_time[i1] - reference_time[i0])
155    print 'Timelag between minima = %.18e' %res
156
157
158
159
160    if plotting is True:
161        ion()
162        hold(False)
163   
164        plot(reference_time, validation_data[name], 'r-',
165             reference_time, model, 'k-')
166        title('Gauge %s' %name)
167        xlabel('time(s)')
168        ylabel('stage (m)')   
169        legend(('Observed', 'Modelled'), shadow=True, loc='upper left')
170        savefig(name, dpi = 300)       
171
172        raw_input('Next')
173
174
175
Note: See TracBrowser for help on using the repository browser.