source: anuga_validation/automated_validation_tests/okushiri_tank_validation/validate_okushiri.py @ 3825

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

Work on validation

File size: 6.2 KB
Line 
1"""Automated validation of ANUGA.
2
3This script runs the benchmark on a regular grid and verifies that
4the predicted timeseries at selected gauges are close enought to the observed
5timeseries provided from the experiment.
6
7See anuga_validation/okushiri_2005 for scripts that produce simulations at
8greater resolution.
9
10"""
11
12# Module imports
13from Numeric import array, zeros, Float, allclose
14
15from anuga.utilities.numerical_tools import ensure_numeric
16from anuga.shallow_water import Domain
17from anuga.shallow_water import Reflective_boundary
18from anuga.shallow_water import Transmissive_Momentum_Set_Stage_boundary
19from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular_cross
20from anuga.abstract_2d_finite_volumes.util import file_function
21
22import project
23
24try:
25    import pylab
26except:
27    plotting = False
28else:
29    plotting = True
30   
31
32finaltime = 22.5
33timestep = 0.05
34
35#-------------------------
36# Create Domain
37#-------------------------
38#N = 50
39#N = 50
40#points, vertices, boundary = rectangular_cross(N, N/5*3,
41#                                               len1=5.448, len2=3.402)
42#domain = Domain(points, vertices, boundary)
43
44
45print 'Creating domain from', project.mesh_filename
46
47domain = Domain(project.mesh_filename, use_cache=True, verbose=True)
48
49
50domain.set_name('okushiri_automated_validation_varmesh')
51domain.set_default_order(2)
52domain.set_minimum_storable_height(0.001)
53domain.set_maximum_allowed_speed(0) # The default in August 2005
54
55
56# Set old (pre Sep 2006) defaults for limiters
57domain.beta_w      = 0.9
58domain.beta_w_dry  = 0.9
59domain.beta_uh     = 0.9
60domain.beta_uh_dry = 0.9
61domain.beta_vh     = 0.9
62domain.beta_vh_dry = 0.9
63
64domain.check_integrity()
65
66print domain.statistics()
67
68#-------------------------
69# Initial Conditions
70#-------------------------
71domain.set_quantity('friction', 0.0)
72domain.set_quantity('stage', 0.0)
73domain.set_quantity('elevation',
74                    filename=project.bathymetry_filename[:-4] + '.pts',
75                    alpha=0.02,                   
76                    verbose=True,
77                    use_cache=True)
78
79#-------------------------
80# Boundary Conditions
81#-------------------------
82Br = Reflective_boundary(domain)
83
84function = file_function(project.boundary_filename[:-4] + '.tms',
85                         domain,
86                         verbose=True)
87
88Bts = Transmissive_Momentum_Set_Stage_boundary(domain, function)
89
90domain.set_boundary({'wave': Bts, 'wall': Br})
91
92
93
94
95#-------------------------
96# Set up validation data
97#-------------------------
98gauge_locations = [[0.000, 1.696]]  #Boundary
99gauge_locations += [[4.521, 1.196], [4.521, 1.696], [4.521, 2.196]] #Ch 5-7-9
100
101#Boundary input wave
102filename = project.boundary_filename[:-4] + '.tms'
103print 'Reading', filename
104from Scientific.IO.NetCDF import NetCDFFile
105from Numeric import allclose
106fid = NetCDFFile(filename, 'r')#
107
108input_time = fid.variables['time'][:]
109input_stage = fid.variables['stage'][:]
110
111
112# reference gauge timeseries
113reference_time = []
114ch5 = []
115ch7 = []
116ch9 = []
117fid = open('output_ch5-7-9.txt')
118lines = fid.readlines()
119fid.close()
120#for i, line in enumerate(lines[1:]):
121for line in lines[1:]:
122    fields = line.split()
123
124    reference_time.append(float(fields[0]))
125    ch5.append(float(fields[1])/100)   #cm2m
126    ch7.append(float(fields[2])/100)   #cm2m
127    ch9.append(float(fields[3])/100)   #cm2m
128
129    if fields[0] == str(finaltime): break   
130   
131
132assert reference_time[0] == 0.0
133assert reference_time[-1] == finaltime
134
135reference_gauge_values = [ensure_numeric(input_stage),
136                          ensure_numeric(ch5),
137                          ensure_numeric(ch7),
138                          ensure_numeric(ch9)] 
139
140
141
142predicted_gauge_values = [[], [], [], []]
143
144
145#-------------------------
146# Evolve through time
147#-------------------------
148import time
149t0 = time.time()
150
151w_max = 0
152for j, t in enumerate(domain.evolve(yieldstep = timestep,
153                                    finaltime = finaltime)):
154    domain.write_time()
155
156    # Record gauge values
157    stage = domain.get_quantity('stage')
158    for i, (x, y) in enumerate(gauge_locations):
159        w = stage.get_values(interpolation_points=[[x,y]])[0]
160        predicted_gauge_values[i].append(w)
161
162        #print 'difference', x, y,\
163        #      predicted_gauge_values[i][j] - reference_gauge_values[i][j]
164
165
166    w = domain.get_maximum_inundation_elevation()
167    x, y = domain.get_maximum_inundation_location()
168    print '  Coastline elevation = %.2f at (x,y)=(%.2f, %.2f)' %(w, x, y)
169    print
170   
171    if w > w_max:
172        w_max = w
173        x_max = x
174        y_max = y
175
176
177print '**********************************************'
178print 'Max coastline elevation = %.2f at (%.2f, %.2f)' %(w_max, x_max, y_max)
179print '**********************************************'
180       
181   
182
183#-------------------------
184# Validate
185#-------------------------
186for i, (x, y) in enumerate(gauge_locations):       
187    predicted_gauge_values[i] = ensure_numeric(predicted_gauge_values[i])
188    predicted_gauge_values[i] = predicted_gauge_values[i][:451]   
189
190    print predicted_gauge_values[i].shape
191    print reference_gauge_values[i].shape
192    #print predicted_gauge_values.shape
193    #print reference_gauge_values.shape   
194   
195    sqsum = sum((predicted_gauge_values[i].flat - \
196                 reference_gauge_values[i].flat)**2)
197    denom = sum(reference_gauge_values[i].flat**2)
198
199    print i, sqsum, sqsum/denom
200
201    msg = 'Rms error is %f' %(sqsum/denom)
202    print msg
203   
204    #assert sqsum/denom < 0.01, msg
205   
206
207    if plotting is True:
208       
209        from pylab import *
210
211        ion()
212        hold(False)
213        plot(reference_time, reference_gauge_values[i], 'r.-',
214             reference_time, predicted_gauge_values[i], 'k-')
215
216        title('Gauge %d (%f,%f)' %(i,x,y))
217        xlabel('time(s)')
218        ylabel('stage (m)')   
219        legend(('Observed', 'Modelled'), shadow=True, loc='upper left')
220        savefig('Gauge_%d.png' %i, dpi = 300)
221       
222        #raw_input('Next')
223       
224
225if plotting is True:
226    show()
227   
228
229print 'That took %.2f seconds' %(time.time()-t0)
Note: See TracBrowser for help on using the repository browser.