1 | import numpy as num |
---|
2 | import unittest |
---|
3 | import tempfile |
---|
4 | import os |
---|
5 | import sys |
---|
6 | from Scientific.IO.NetCDF import NetCDFFile |
---|
7 | |
---|
8 | from anuga.utilities.system_tools import get_pathname_from_package |
---|
9 | from anuga.coordinate_transforms.geo_reference import Geo_reference |
---|
10 | from anuga.coordinate_transforms.redfearn import redfearn |
---|
11 | from anuga.utilities.numerical_tools import ensure_numeric |
---|
12 | from anuga.config import netcdf_mode_r, netcdf_mode_w, netcdf_mode_a |
---|
13 | from anuga.file.sts import create_sts_boundary |
---|
14 | from anuga.file.csv_file import load_csv_as_dict, load_csv_as_array |
---|
15 | |
---|
16 | from anuga.shallow_water.shallow_water_domain import Domain |
---|
17 | |
---|
18 | # boundary functions |
---|
19 | from anuga.shallow_water.boundaries import Reflective_boundary, \ |
---|
20 | Field_boundary, Transmissive_momentum_set_stage_boundary, \ |
---|
21 | Transmissive_stage_zero_momentum_boundary |
---|
22 | from anuga.abstract_2d_finite_volumes.generic_boundary_conditions\ |
---|
23 | import Transmissive_boundary, Dirichlet_boundary, \ |
---|
24 | Time_boundary, File_boundary, AWI_boundary |
---|
25 | |
---|
26 | from anuga.pmesh.mesh_interface import create_mesh_from_regions |
---|
27 | |
---|
28 | from urs2sts import urs2sts |
---|
29 | |
---|
30 | # Allow us to use helper methods from this test. |
---|
31 | from anuga.file.test_mux import Test_Mux |
---|
32 | |
---|
33 | class Test_Urs2Sts(Test_Mux): |
---|
34 | """ A suite of tests to test urs2sts file conversion functions. |
---|
35 | These tests are quite coarse-grained: converting a file |
---|
36 | and checking that its headers and some of its contents |
---|
37 | are correct. |
---|
38 | """ |
---|
39 | |
---|
40 | def test_urs2sts0(self): |
---|
41 | """ |
---|
42 | Test single source |
---|
43 | """ |
---|
44 | tide=0 |
---|
45 | time_step_count = 3 |
---|
46 | time_step = 2 |
---|
47 | lat_long_points =[(-21.5,114.5),(-21,114.5),(-21.5,115), (-21.,115.)] |
---|
48 | n=len(lat_long_points) |
---|
49 | first_tstep=num.ones(n,num.int) |
---|
50 | first_tstep[0]+=1 |
---|
51 | first_tstep[2]+=1 |
---|
52 | last_tstep=(time_step_count)*num.ones(n,num.int) |
---|
53 | last_tstep[0]-=1 |
---|
54 | |
---|
55 | gauge_depth=20*num.ones(n,num.float) |
---|
56 | ha=2*num.ones((n,time_step_count),num.float) |
---|
57 | ha[0]=num.arange(0,time_step_count) |
---|
58 | ha[1]=num.arange(time_step_count,2*time_step_count) |
---|
59 | ha[2]=num.arange(2*time_step_count,3*time_step_count) |
---|
60 | ha[3]=num.arange(3*time_step_count,4*time_step_count) |
---|
61 | ua=5*num.ones((n,time_step_count),num.float) |
---|
62 | va=-10*num.ones((n,time_step_count),num.float) |
---|
63 | |
---|
64 | base_name, files = self.write_mux2(lat_long_points, |
---|
65 | time_step_count, time_step, |
---|
66 | first_tstep, last_tstep, |
---|
67 | depth=gauge_depth, |
---|
68 | ha=ha, |
---|
69 | ua=ua, |
---|
70 | va=va) |
---|
71 | |
---|
72 | sts_file = base_name + '.sts' |
---|
73 | |
---|
74 | urs2sts(base_name, |
---|
75 | basename_out=sts_file, |
---|
76 | mean_stage=tide,verbose=False) |
---|
77 | |
---|
78 | |
---|
79 | #Let's interigate the sww file |
---|
80 | # Note, the sww info is not gridded. It is point data. |
---|
81 | fid = NetCDFFile(sts_file) |
---|
82 | |
---|
83 | # Make x and y absolute |
---|
84 | x = fid.variables['x'][:] |
---|
85 | y = fid.variables['y'][:] |
---|
86 | |
---|
87 | geo_reference = Geo_reference(NetCDFObject=fid) |
---|
88 | points = geo_reference.get_absolute(map(None, x, y)) |
---|
89 | points = ensure_numeric(points) |
---|
90 | |
---|
91 | x = points[:,0] |
---|
92 | y = points[:,1] |
---|
93 | |
---|
94 | #Check that first coordinate is correctly represented |
---|
95 | #Work out the UTM coordinates for first point |
---|
96 | for i in range(4): |
---|
97 | zone, e, n = redfearn(lat_long_points[i][0], lat_long_points[i][1]) |
---|
98 | assert num.allclose([x[i],y[i]], [e,n]) |
---|
99 | |
---|
100 | #Check the time vector |
---|
101 | times = fid.variables['time'][:] |
---|
102 | |
---|
103 | times_actual = [] |
---|
104 | for i in range(time_step_count): |
---|
105 | times_actual.append(time_step * i) |
---|
106 | |
---|
107 | assert num.allclose(ensure_numeric(times), |
---|
108 | ensure_numeric(times_actual)) |
---|
109 | |
---|
110 | #Check first value |
---|
111 | stage = fid.variables['stage'][:] |
---|
112 | xmomentum = fid.variables['xmomentum'][:] |
---|
113 | ymomentum = fid.variables['ymomentum'][:] |
---|
114 | elevation = fid.variables['elevation'][:] |
---|
115 | |
---|
116 | # Set original data used to write mux file to be zero when gauges are |
---|
117 | #not recdoring |
---|
118 | ha[0][0]=0.0 |
---|
119 | ha[0][time_step_count-1]=0.0; |
---|
120 | ha[2][0]=0.0; |
---|
121 | ua[0][0]=0.0 |
---|
122 | ua[0][time_step_count-1]=0.0; |
---|
123 | ua[2][0]=0.0; |
---|
124 | va[0][0]=0.0 |
---|
125 | va[0][time_step_count-1]=0.0; |
---|
126 | va[2][0]=0.0; |
---|
127 | |
---|
128 | assert num.allclose(num.transpose(ha),stage) #Meters |
---|
129 | |
---|
130 | #Check the momentums - ua |
---|
131 | #momentum = velocity*(stage-elevation) |
---|
132 | # elevation = - depth |
---|
133 | #momentum = velocity_ua *(stage+depth) |
---|
134 | |
---|
135 | depth=num.zeros((len(lat_long_points),time_step_count),num.float) |
---|
136 | for i in range(len(lat_long_points)): |
---|
137 | depth[i]=gauge_depth[i]+tide+ha[i] |
---|
138 | assert num.allclose(num.transpose(ua*depth),xmomentum) |
---|
139 | |
---|
140 | #Check the momentums - va |
---|
141 | #momentum = velocity*(stage-elevation) |
---|
142 | # elevation = - depth |
---|
143 | #momentum = velocity_va *(stage+depth) |
---|
144 | |
---|
145 | assert num.allclose(num.transpose(va*depth),ymomentum) |
---|
146 | |
---|
147 | # check the elevation values. |
---|
148 | # -ve since urs measures depth, sww meshers height, |
---|
149 | assert num.allclose(-elevation, gauge_depth) #Meters |
---|
150 | |
---|
151 | fid.close() |
---|
152 | self.delete_mux(files) |
---|
153 | os.remove(sts_file) |
---|
154 | |
---|
155 | def test_urs2sts_nonstandard_meridian(self): |
---|
156 | """ |
---|
157 | Test single source using the meridian from zone 50 as a nonstandard meridian |
---|
158 | """ |
---|
159 | tide=0 |
---|
160 | time_step_count = 3 |
---|
161 | time_step = 2 |
---|
162 | lat_long_points =[(-21.,114.5),(-21.,113.5),(-21.,114.), (-21.,115.)] |
---|
163 | n=len(lat_long_points) |
---|
164 | first_tstep=num.ones(n,num.int) |
---|
165 | first_tstep[0]+=1 |
---|
166 | first_tstep[2]+=1 |
---|
167 | last_tstep=(time_step_count)*num.ones(n,num.int) |
---|
168 | last_tstep[0]-=1 |
---|
169 | |
---|
170 | gauge_depth=20*num.ones(n,num.float) |
---|
171 | ha=2*num.ones((n,time_step_count),num.float) |
---|
172 | ha[0]=num.arange(0,time_step_count) |
---|
173 | ha[1]=num.arange(time_step_count,2*time_step_count) |
---|
174 | ha[2]=num.arange(2*time_step_count,3*time_step_count) |
---|
175 | ha[3]=num.arange(3*time_step_count,4*time_step_count) |
---|
176 | ua=5*num.ones((n,time_step_count),num.float) |
---|
177 | va=-10*num.ones((n,time_step_count),num.float) |
---|
178 | |
---|
179 | base_name, files = self.write_mux2(lat_long_points, |
---|
180 | time_step_count, time_step, |
---|
181 | first_tstep, last_tstep, |
---|
182 | depth=gauge_depth, |
---|
183 | ha=ha, |
---|
184 | ua=ua, |
---|
185 | va=va) |
---|
186 | |
---|
187 | sts_file = base_name + '.sts' |
---|
188 | |
---|
189 | urs2sts(base_name, |
---|
190 | basename_out=sts_file, |
---|
191 | central_meridian=123, |
---|
192 | mean_stage=tide, |
---|
193 | verbose=False) |
---|
194 | |
---|
195 | #Let's interigate the sww file |
---|
196 | # Note, the sww info is not gridded. It is point data. |
---|
197 | fid = NetCDFFile(sts_file) |
---|
198 | |
---|
199 | # Make x and y absolute |
---|
200 | x = fid.variables['x'][:] |
---|
201 | y = fid.variables['y'][:] |
---|
202 | |
---|
203 | geo_reference = Geo_reference(NetCDFObject=fid) |
---|
204 | points = geo_reference.get_absolute(map(None, x, y)) |
---|
205 | points = ensure_numeric(points) |
---|
206 | |
---|
207 | x = points[:,0] |
---|
208 | y = points[:,1] |
---|
209 | |
---|
210 | # Check that all coordinate are correctly represented |
---|
211 | # Using the non standard projection (50) |
---|
212 | for i in range(4): |
---|
213 | zone, e, n = redfearn(lat_long_points[i][0], |
---|
214 | lat_long_points[i][1], |
---|
215 | central_meridian=123) |
---|
216 | assert num.allclose([x[i],y[i]], [e,n]) |
---|
217 | assert zone==-1 |
---|
218 | self.delete_mux(files) |
---|
219 | |
---|
220 | |
---|
221 | |
---|
222 | def test_urs2sts_individual_sources(self): |
---|
223 | """Test that individual sources compare to actual urs output |
---|
224 | Test that the first recording time is the smallest |
---|
225 | over waveheight, easting and northing velocity |
---|
226 | """ |
---|
227 | |
---|
228 | # Get path where this test is run |
---|
229 | path = get_pathname_from_package('anuga.shallow_water') |
---|
230 | |
---|
231 | testdir = os.path.join(path, 'urs_test_data') |
---|
232 | ordering_filename=os.path.join(testdir, 'thinned_bound_order_test.txt') |
---|
233 | |
---|
234 | sources = ['1-z.grd','2-z.grd','3-z.grd'] |
---|
235 | |
---|
236 | # Start times by source and station taken manually from urs header files |
---|
237 | time_start_z = num.array([[10.0,11.5,13,14.5,17.7], |
---|
238 | [9.8,11.2,12.7,14.2,17.4], |
---|
239 | [9.5,10.9,12.4,13.9,17.1]]) |
---|
240 | |
---|
241 | time_start_e = time_start_n = time_start_z |
---|
242 | |
---|
243 | # time step in urs output |
---|
244 | delta_t = 0.1 |
---|
245 | |
---|
246 | # Make sts file for each source |
---|
247 | for k, source_filename in enumerate(sources): |
---|
248 | source_number = k + 1 # Source numbering starts at 1 |
---|
249 | |
---|
250 | urs_filenames = os.path.join(testdir, source_filename) |
---|
251 | weights = [1.] |
---|
252 | sts_name_out = 'test' |
---|
253 | |
---|
254 | urs2sts(urs_filenames, |
---|
255 | basename_out=sts_name_out, |
---|
256 | ordering_filename=ordering_filename, |
---|
257 | weights=weights, |
---|
258 | mean_stage=0.0, |
---|
259 | verbose=False) |
---|
260 | |
---|
261 | # Read in sts file for this source file |
---|
262 | fid = NetCDFFile(sts_name_out+'.sts', netcdf_mode_r) # Open existing file for read |
---|
263 | x = fid.variables['x'][:]+fid.xllcorner # x-coordinates of vertices |
---|
264 | y = fid.variables['y'][:]+fid.yllcorner # y-coordinates of vertices |
---|
265 | elevation = fid.variables['elevation'][:] |
---|
266 | time=fid.variables['time'][:]+fid.starttime |
---|
267 | |
---|
268 | # Get quantity data from sts file |
---|
269 | quantity_names=['stage','xmomentum','ymomentum'] |
---|
270 | quantities = {} |
---|
271 | for i, name in enumerate(quantity_names): |
---|
272 | quantities[name] = fid.variables[name][:] |
---|
273 | |
---|
274 | # Make sure start time from sts file is the minimum starttime |
---|
275 | # across all stations (for this source) |
---|
276 | #print k, time_start_z[k,:] |
---|
277 | starttime = min(time_start_z[k, :]) |
---|
278 | sts_starttime = fid.starttime[0] |
---|
279 | msg = 'sts starttime for source %d was %f. Should have been %f'\ |
---|
280 | %(source_number, sts_starttime, starttime) |
---|
281 | assert num.allclose(sts_starttime, starttime), msg |
---|
282 | |
---|
283 | # For each station, compare urs2sts output to known urs output |
---|
284 | for j in range(len(x)): |
---|
285 | index_start_urs = 0 |
---|
286 | index_end_urs = 0 |
---|
287 | index_start = 0 |
---|
288 | index_end = 0 |
---|
289 | count = 0 |
---|
290 | |
---|
291 | # read in urs test data for stage, e and n velocity |
---|
292 | urs_file_name_z = 'z_'+str(source_number)+'_'+str(j)+'.csv' |
---|
293 | dict = load_csv_as_array(os.path.join(testdir, urs_file_name_z)) |
---|
294 | urs_stage = dict['urs_stage'] |
---|
295 | urs_file_name_e = 'e_'+str(source_number)+'_'+str(j)+'.csv' |
---|
296 | dict = load_csv_as_array(os.path.join(testdir, urs_file_name_e)) |
---|
297 | urs_e = dict['urs_e'] |
---|
298 | urs_file_name_n = 'n_'+str(source_number)+'_'+str(j)+'.csv' |
---|
299 | dict = load_csv_as_array(os.path.join(testdir, urs_file_name_n)) |
---|
300 | urs_n = dict['urs_n'] |
---|
301 | |
---|
302 | # find start and end time for stage |
---|
303 | for i in range(len(urs_stage)): |
---|
304 | if urs_stage[i] == 0.0: |
---|
305 | index_start_urs_z = i+1 |
---|
306 | if int(urs_stage[i]) == 99 and count <> 1: |
---|
307 | count +=1 |
---|
308 | index_end_urs_z = i |
---|
309 | |
---|
310 | if count == 0: index_end_urs_z = len(urs_stage) |
---|
311 | |
---|
312 | start_times_z = time_start_z[source_number-1] |
---|
313 | |
---|
314 | # start times for easting velocities should be the same as stage |
---|
315 | start_times_e = time_start_e[source_number-1] |
---|
316 | index_start_urs_e = index_start_urs_z |
---|
317 | index_end_urs_e = index_end_urs_z |
---|
318 | |
---|
319 | # start times for northing velocities should be the same as stage |
---|
320 | start_times_n = time_start_n[source_number-1] |
---|
321 | index_start_urs_n = index_start_urs_z |
---|
322 | index_end_urs_n = index_end_urs_z |
---|
323 | |
---|
324 | # Check that actual start time matches header information for stage |
---|
325 | msg = 'stage start time from urs file is not the same as the ' |
---|
326 | msg += 'header file for source %i and station %i' %(source_number,j) |
---|
327 | assert num.allclose(index_start_urs_z,start_times_z[j]/delta_t), msg |
---|
328 | |
---|
329 | msg = 'e velocity start time from urs file is not the same as the ' |
---|
330 | msg += 'header file for source %i and station %i' %(source_number,j) |
---|
331 | assert num.allclose(index_start_urs_e,start_times_e[j]/delta_t), msg |
---|
332 | |
---|
333 | msg = 'n velocity start time from urs file is not the same as the ' |
---|
334 | msg += 'header file for source %i and station %i' %(source_number,j) |
---|
335 | assert num.allclose(index_start_urs_n,start_times_n[j]/delta_t), msg |
---|
336 | |
---|
337 | # get index for start and end time for sts quantities |
---|
338 | index_start_stage = 0 |
---|
339 | index_end_stage = 0 |
---|
340 | count = 0 |
---|
341 | sts_stage = quantities['stage'][:,j] |
---|
342 | for i in range(len(sts_stage)): |
---|
343 | if sts_stage[i] <> 0.0 and count <> 1: |
---|
344 | count += 1 |
---|
345 | index_start_stage = i |
---|
346 | if int(sts_stage[i]) == 99 and count <> 1: |
---|
347 | count += 1 |
---|
348 | index_end_stage = i |
---|
349 | |
---|
350 | index_end_stage = index_start_stage + len(urs_stage[index_start_urs_z:index_end_urs_z]) |
---|
351 | |
---|
352 | sts_xmom = quantities['xmomentum'][:,j] |
---|
353 | index_start_x = index_start_stage |
---|
354 | index_end_x = index_start_x + len(urs_e[index_start_urs_e:index_end_urs_e]) |
---|
355 | |
---|
356 | sts_ymom = quantities['ymomentum'][:,j] |
---|
357 | index_start_y = index_start_stage |
---|
358 | index_end_y = index_start_y + len(urs_n[index_start_urs_n:index_end_urs_n]) |
---|
359 | |
---|
360 | # check that urs stage and sts stage are the same |
---|
361 | msg = 'urs stage is not equal to sts stage for for source %i and station %i' %(source_number,j) |
---|
362 | assert num.allclose(urs_stage[index_start_urs_z:index_end_urs_z], |
---|
363 | sts_stage[index_start_stage:index_end_stage], |
---|
364 | rtol=1.0e-6, atol=1.0e-5 ), msg |
---|
365 | |
---|
366 | # check that urs e velocity and sts xmomentum are the same |
---|
367 | msg = 'urs e velocity is not equivalent to sts x momentum for for source %i and station %i' %(source_number,j) |
---|
368 | assert num.allclose(urs_e[index_start_urs_e:index_end_urs_e]*(urs_stage[index_start_urs_e:index_end_urs_e]-elevation[j]), |
---|
369 | sts_xmom[index_start_x:index_end_x], |
---|
370 | rtol=1.0e-5, atol=1.0e-4 ), msg |
---|
371 | |
---|
372 | # check that urs n velocity and sts ymomentum are the same |
---|
373 | #print 'urs n velocity', urs_n[index_start_urs_n:index_end_urs_n]*(urs_stage[index_start_urs_n:index_end_urs_n]-elevation[j]) |
---|
374 | #print 'sts momentum', sts_ymom[index_start_y:index_end_y] |
---|
375 | msg = 'urs n velocity is not equivalent to sts y momentum for source %i and station %i' %(source_number,j) |
---|
376 | assert num.allclose(urs_n[index_start_urs_n:index_end_urs_n]*(urs_stage[index_start_urs_n:index_end_urs_n]-elevation[j]), |
---|
377 | -sts_ymom[index_start_y:index_end_y], |
---|
378 | rtol=1.0e-5, atol=1.0e-4 ), msg |
---|
379 | |
---|
380 | |
---|
381 | fid.close() |
---|
382 | |
---|
383 | os.remove(sts_name_out+'.sts') |
---|
384 | |
---|
385 | def test_urs2sts_combined_sources(self): |
---|
386 | """Test that combined sources compare to actual urs output |
---|
387 | Test that the first recording time is the smallest |
---|
388 | over waveheight, easting and northing velocity |
---|
389 | """ |
---|
390 | |
---|
391 | # combined |
---|
392 | time_start_z = num.array([9.5,10.9,12.4,13.9,17.1]) |
---|
393 | time_start_e = time_start_n = time_start_z |
---|
394 | |
---|
395 | # make sts file for combined sources |
---|
396 | weights = [1., 2., 3.] |
---|
397 | |
---|
398 | path = get_pathname_from_package('anuga.shallow_water') |
---|
399 | |
---|
400 | testdir = os.path.join(path, 'urs_test_data') |
---|
401 | ordering_filename=os.path.join(testdir, 'thinned_bound_order_test.txt') |
---|
402 | |
---|
403 | urs_filenames = [os.path.join(testdir,'1-z.grd'), |
---|
404 | os.path.join(testdir,'2-z.grd'), |
---|
405 | os.path.join(testdir,'3-z.grd')] |
---|
406 | sts_name_out = 'test' |
---|
407 | |
---|
408 | urs2sts(urs_filenames, |
---|
409 | basename_out=sts_name_out, |
---|
410 | ordering_filename=ordering_filename, |
---|
411 | weights=weights, |
---|
412 | mean_stage=0.0, |
---|
413 | verbose=False) |
---|
414 | |
---|
415 | # read in sts file for combined source |
---|
416 | fid = NetCDFFile(sts_name_out+'.sts', netcdf_mode_r) # Open existing file for read |
---|
417 | x = fid.variables['x'][:]+fid.xllcorner # x-coordinates of vertices |
---|
418 | y = fid.variables['y'][:]+fid.yllcorner # y-coordinates of vertices |
---|
419 | elevation = fid.variables['elevation'][:] |
---|
420 | time=fid.variables['time'][:]+fid.starttime |
---|
421 | |
---|
422 | |
---|
423 | # Check that stored permutation is as per default |
---|
424 | permutation = range(len(x)) |
---|
425 | stored_permutation = fid.variables['permutation'][:] |
---|
426 | msg = 'Permutation was not stored correctly. I got ' |
---|
427 | msg += str(stored_permutation) |
---|
428 | assert num.allclose(stored_permutation, permutation), msg |
---|
429 | |
---|
430 | # Get quantity data from sts file |
---|
431 | quantity_names=['stage','xmomentum','ymomentum'] |
---|
432 | quantities = {} |
---|
433 | for i, name in enumerate(quantity_names): |
---|
434 | quantities[name] = fid.variables[name][:] |
---|
435 | |
---|
436 | # For each station, compare urs2sts output to known urs output |
---|
437 | delta_t = 0.1 |
---|
438 | |
---|
439 | # Make sure start time from sts file is the minimum starttime |
---|
440 | # across all stations (for this source) |
---|
441 | starttime = min(time_start_z[:]) |
---|
442 | sts_starttime = fid.starttime[0] |
---|
443 | msg = 'sts starttime was %f. Should have been %f'\ |
---|
444 | %(sts_starttime, starttime) |
---|
445 | assert num.allclose(sts_starttime, starttime), msg |
---|
446 | |
---|
447 | #stations = [1,2,3] |
---|
448 | #for j in stations: |
---|
449 | for j in range(len(x)): |
---|
450 | index_start_urs_z = 0 |
---|
451 | index_end_urs_z = 0 |
---|
452 | index_start_urs_e = 0 |
---|
453 | index_end_urs_e = 0 |
---|
454 | index_start_urs_n = 0 |
---|
455 | index_end_urs_n = 0 |
---|
456 | count = 0 |
---|
457 | |
---|
458 | # read in urs test data for stage, e and n velocity |
---|
459 | urs_file_name_z = 'z_combined_'+str(j)+'.csv' |
---|
460 | dict = load_csv_as_array(os.path.join(testdir, urs_file_name_z)) |
---|
461 | urs_stage = dict['urs_stage'] |
---|
462 | urs_file_name_e = 'e_combined_'+str(j)+'.csv' |
---|
463 | dict = load_csv_as_array(os.path.join(testdir, urs_file_name_e)) |
---|
464 | urs_e = dict['urs_e'] |
---|
465 | urs_file_name_n = 'n_combined_'+str(j)+'.csv' |
---|
466 | dict = load_csv_as_array(os.path.join(testdir, urs_file_name_n)) |
---|
467 | urs_n = dict['urs_n'] |
---|
468 | |
---|
469 | # find start and end time for stage |
---|
470 | for i in range(len(urs_stage)): |
---|
471 | if urs_stage[i] == 0.0: |
---|
472 | index_start_urs_z = i+1 |
---|
473 | if int(urs_stage[i]) == 99 and count <> 1: |
---|
474 | count +=1 |
---|
475 | index_end_urs_z = i |
---|
476 | |
---|
477 | if count == 0: index_end_urs_z = len(urs_stage) |
---|
478 | |
---|
479 | start_times_z = time_start_z[j] |
---|
480 | |
---|
481 | start_times_e = time_start_e[j] |
---|
482 | index_start_urs_e = index_start_urs_z |
---|
483 | |
---|
484 | start_times_n = time_start_n[j] |
---|
485 | index_start_urs_n = index_start_urs_z |
---|
486 | |
---|
487 | # Check that actual start time matches header information for stage |
---|
488 | msg = 'stage start time from urs file is not the same as the ' |
---|
489 | msg += 'header file at station %i' %(j) |
---|
490 | assert num.allclose(index_start_urs_z,start_times_z/delta_t), msg |
---|
491 | |
---|
492 | msg = 'e velocity start time from urs file is not the same as the ' |
---|
493 | msg += 'header file at station %i' %(j) |
---|
494 | assert num.allclose(index_start_urs_e,start_times_e/delta_t), msg |
---|
495 | |
---|
496 | msg = 'n velocity start time from urs file is not the same as the ' |
---|
497 | msg += 'header file at station %i' %(j) |
---|
498 | assert num.allclose(index_start_urs_n,start_times_n/delta_t), msg |
---|
499 | |
---|
500 | # get index for start and end time for sts quantities |
---|
501 | index_start_stage = 0 |
---|
502 | index_end_stage = 0 |
---|
503 | index_start_x = 0 |
---|
504 | index_end_x = 0 |
---|
505 | index_start_y = 0 |
---|
506 | index_end_y = 0 |
---|
507 | count = 0 |
---|
508 | count1 = 0 |
---|
509 | sts_stage = quantities['stage'][:,j] |
---|
510 | for i in range(len(sts_stage)): |
---|
511 | if sts_stage[i] <> 0.0 and count <> 1: |
---|
512 | count += 1 |
---|
513 | index_start_stage = i |
---|
514 | if int(urs_stage[i]) == 99 and count <> 1: |
---|
515 | count +=1 |
---|
516 | index_end_stage = i |
---|
517 | |
---|
518 | index_end_stage = index_start_stage + len(urs_stage[index_start_urs_z:index_end_urs_z]) |
---|
519 | |
---|
520 | index_start_x = index_start_stage |
---|
521 | index_end_x = index_start_x + len(urs_stage[index_start_urs_e:index_end_urs_e]) |
---|
522 | sts_xmom = quantities['ymomentum'][:,j] |
---|
523 | |
---|
524 | index_start_y = index_start_stage |
---|
525 | index_end_y = index_start_y + len(urs_stage[index_start_urs_n:index_end_urs_n]) |
---|
526 | sts_ymom = quantities['ymomentum'][:,j] |
---|
527 | |
---|
528 | # check that urs stage and sts stage are the same |
---|
529 | msg = 'urs stage is not equal to sts stage for station %i' %j |
---|
530 | #print 'urs stage', urs_stage[index_start_urs_z:index_end_urs_z] |
---|
531 | #print 'sts stage', sts_stage[index_start_stage:index_end_stage] |
---|
532 | #print 'diff', max(urs_stage[index_start_urs_z:index_end_urs_z]-sts_stage[index_start_stage:index_end_stage]) |
---|
533 | #print 'index', index_start_stage, index_end_stage, len(sts_stage) |
---|
534 | assert num.allclose(urs_stage[index_start_urs_z:index_end_urs_z], |
---|
535 | sts_stage[index_start_stage:index_end_stage], |
---|
536 | rtol=1.0e-5, atol=1.0e-4 ), msg |
---|
537 | |
---|
538 | # check that urs e velocity and sts xmomentum are the same |
---|
539 | msg = 'urs e velocity is not equivalent to sts xmomentum for station %i' %j |
---|
540 | assert num.allclose(urs_e[index_start_urs_e:index_end_urs_e]*(urs_stage[index_start_urs_e:index_end_urs_e]-elevation[j]), |
---|
541 | sts_xmom[index_start_x:index_end_x], |
---|
542 | rtol=1.0e-5, atol=1.0e-4 ), msg |
---|
543 | |
---|
544 | # check that urs n velocity and sts ymomentum are the same |
---|
545 | msg = 'urs n velocity is not equivalent to sts ymomentum for station %i' %j |
---|
546 | assert num.allclose(urs_n[index_start_urs_n:index_end_urs_n]*(urs_stage[index_start_urs_n:index_end_urs_n]-elevation[j]), |
---|
547 | sts_ymom[index_start_y:index_end_y], |
---|
548 | rtol=1.0e-5, atol=1.0e-4 ), msg |
---|
549 | |
---|
550 | fid.close() |
---|
551 | |
---|
552 | os.remove(sts_name_out+'.sts') |
---|
553 | |
---|
554 | |
---|
555 | |
---|
556 | def test_urs2sts_ordering(self): |
---|
557 | """Test multiple sources with ordering file |
---|
558 | """ |
---|
559 | |
---|
560 | tide = 0.35 |
---|
561 | time_step_count = 6 # I made this a little longer (Ole) |
---|
562 | time_step = 2 |
---|
563 | lat_long_points =[(-21.5,114.5),(-21,114.5),(-21.5,115), (-21.,115.)] |
---|
564 | n=len(lat_long_points) |
---|
565 | first_tstep=num.ones(n,num.int) |
---|
566 | first_tstep[0]+=1 |
---|
567 | first_tstep[2]+=1 |
---|
568 | last_tstep=(time_step_count)*num.ones(n,num.int) |
---|
569 | last_tstep[0]-=1 |
---|
570 | |
---|
571 | gauge_depth=20*num.ones(n,num.float) |
---|
572 | ha=2*num.ones((n,time_step_count),num.float) |
---|
573 | ha[0]=num.arange(0,time_step_count) |
---|
574 | ha[1]=num.arange(time_step_count,2*time_step_count) |
---|
575 | ha[2]=num.arange(2*time_step_count,3*time_step_count) |
---|
576 | ha[3]=num.arange(3*time_step_count,4*time_step_count) |
---|
577 | ua=5*num.ones((n,time_step_count),num.float) |
---|
578 | va=-10*num.ones((n,time_step_count),num.float) |
---|
579 | |
---|
580 | # Create two identical mux files to be combined by urs2sts |
---|
581 | base_nameI, filesI = self.write_mux2(lat_long_points, |
---|
582 | time_step_count, time_step, |
---|
583 | first_tstep, last_tstep, |
---|
584 | depth=gauge_depth, |
---|
585 | ha=ha, |
---|
586 | ua=ua, |
---|
587 | va=va) |
---|
588 | |
---|
589 | base_nameII, filesII = self.write_mux2(lat_long_points, |
---|
590 | time_step_count, time_step, |
---|
591 | first_tstep, last_tstep, |
---|
592 | depth=gauge_depth, |
---|
593 | ha=ha, |
---|
594 | ua=ua, |
---|
595 | va=va) |
---|
596 | |
---|
597 | |
---|
598 | # Create ordering file |
---|
599 | permutation = [3,0,2] |
---|
600 | |
---|
601 | _, ordering_filename = tempfile.mkstemp('') |
---|
602 | order_fid = open(ordering_filename, 'w') |
---|
603 | order_fid.write('index, longitude, latitude\n') |
---|
604 | for index in permutation: |
---|
605 | order_fid.write('%d, %f, %f\n' %(index, |
---|
606 | lat_long_points[index][1], |
---|
607 | lat_long_points[index][0])) |
---|
608 | order_fid.close() |
---|
609 | |
---|
610 | |
---|
611 | |
---|
612 | |
---|
613 | # Call urs2sts with multiple mux files |
---|
614 | urs2sts([base_nameI, base_nameII], |
---|
615 | basename_out=base_nameI, |
---|
616 | ordering_filename=ordering_filename, |
---|
617 | weights=[1.0, 1.0], |
---|
618 | mean_stage=tide, |
---|
619 | verbose=False) |
---|
620 | |
---|
621 | # now I want to check the sts file ... |
---|
622 | sts_file = base_nameI + '.sts' |
---|
623 | |
---|
624 | #Let's interrogate the sts file |
---|
625 | # Note, the sts info is not gridded. It is point data. |
---|
626 | fid = NetCDFFile(sts_file) |
---|
627 | |
---|
628 | # Check that original indices have been stored |
---|
629 | stored_permutation = fid.variables['permutation'][:] |
---|
630 | msg = 'Permutation was not stored correctly. I got ' |
---|
631 | msg += str(stored_permutation) |
---|
632 | assert num.allclose(stored_permutation, permutation), msg |
---|
633 | |
---|
634 | |
---|
635 | # Make x and y absolute |
---|
636 | x = fid.variables['x'][:] |
---|
637 | y = fid.variables['y'][:] |
---|
638 | |
---|
639 | geo_reference = Geo_reference(NetCDFObject=fid) |
---|
640 | points = geo_reference.get_absolute(map(None, x, y)) |
---|
641 | points = ensure_numeric(points) |
---|
642 | |
---|
643 | x = points[:,0] |
---|
644 | y = points[:,1] |
---|
645 | |
---|
646 | #print |
---|
647 | #print x |
---|
648 | #print y |
---|
649 | for i, index in enumerate(permutation): |
---|
650 | # Check that STS points are stored in the correct order |
---|
651 | |
---|
652 | # Work out the UTM coordinates sts point i |
---|
653 | zone, e, n = redfearn(lat_long_points[index][0], |
---|
654 | lat_long_points[index][1]) |
---|
655 | |
---|
656 | #print i, [x[i],y[i]], [e,n] |
---|
657 | assert num.allclose([x[i],y[i]], [e,n]) |
---|
658 | |
---|
659 | |
---|
660 | # Check the time vector |
---|
661 | times = fid.variables['time'][:] |
---|
662 | |
---|
663 | times_actual = [] |
---|
664 | for i in range(time_step_count): |
---|
665 | times_actual.append(time_step * i) |
---|
666 | |
---|
667 | assert num.allclose(ensure_numeric(times), |
---|
668 | ensure_numeric(times_actual)) |
---|
669 | |
---|
670 | |
---|
671 | # Check sts values |
---|
672 | stage = fid.variables['stage'][:] |
---|
673 | xmomentum = fid.variables['xmomentum'][:] |
---|
674 | ymomentum = fid.variables['ymomentum'][:] |
---|
675 | elevation = fid.variables['elevation'][:] |
---|
676 | |
---|
677 | # Set original data used to write mux file to be zero when gauges are |
---|
678 | # not recdoring |
---|
679 | |
---|
680 | ha[0][0]=0.0 |
---|
681 | ha[0][time_step_count-1]=0.0 |
---|
682 | ha[2][0]=0.0 |
---|
683 | ua[0][0]=0.0 |
---|
684 | ua[0][time_step_count-1]=0.0 |
---|
685 | ua[2][0]=0.0 |
---|
686 | va[0][0]=0.0 |
---|
687 | va[0][time_step_count-1]=0.0 |
---|
688 | va[2][0]=0.0; |
---|
689 | |
---|
690 | |
---|
691 | # The stage stored in the .sts file should be the sum of the stage |
---|
692 | # in the two mux2 files because both have weights = 1. In this case |
---|
693 | # the mux2 files are the same so stage == 2.0 * ha |
---|
694 | #print 2.0*num.transpose(ha) - stage |
---|
695 | |
---|
696 | ha_permutation = num.take(ha, permutation, axis=0) |
---|
697 | ua_permutation = num.take(ua, permutation, axis=0) |
---|
698 | va_permutation = num.take(va, permutation, axis=0) |
---|
699 | gauge_depth_permutation = num.take(gauge_depth, permutation, axis=0) |
---|
700 | |
---|
701 | |
---|
702 | assert num.allclose(2.0*num.transpose(ha_permutation)+tide, stage) # Meters |
---|
703 | |
---|
704 | #Check the momentums - ua |
---|
705 | #momentum = velocity*(stage-elevation) |
---|
706 | # elevation = - depth |
---|
707 | #momentum = velocity_ua *(stage+depth) |
---|
708 | |
---|
709 | depth=num.zeros((len(lat_long_points),time_step_count),num.float) |
---|
710 | for i in range(len(lat_long_points)): |
---|
711 | depth[i]=gauge_depth[i]+tide+2.0*ha[i] |
---|
712 | #2.0*ha necessary because using two files with weights=1 are used |
---|
713 | |
---|
714 | depth_permutation = num.take(depth, permutation, axis=0) |
---|
715 | |
---|
716 | |
---|
717 | # The xmomentum stored in the .sts file should be the sum of the ua |
---|
718 | # in the two mux2 files multiplied by the depth. |
---|
719 | assert num.allclose(2.0*num.transpose(ua_permutation*depth_permutation), xmomentum) |
---|
720 | |
---|
721 | #Check the momentums - va |
---|
722 | #momentum = velocity*(stage-elevation) |
---|
723 | # elevation = - depth |
---|
724 | #momentum = velocity_va *(stage+depth) |
---|
725 | |
---|
726 | # The ymomentum stored in the .sts file should be the sum of the va |
---|
727 | # in the two mux2 files multiplied by the depth. |
---|
728 | assert num.allclose(2.0*num.transpose(va_permutation*depth_permutation), ymomentum) |
---|
729 | |
---|
730 | # check the elevation values. |
---|
731 | # -ve since urs measures depth, sww meshers height, |
---|
732 | assert num.allclose(-gauge_depth_permutation, elevation) #Meters |
---|
733 | |
---|
734 | fid.close() |
---|
735 | self.delete_mux(filesI) |
---|
736 | self.delete_mux(filesII) |
---|
737 | os.remove(sts_file) |
---|
738 | |
---|
739 | |
---|
740 | |
---|
741 | |
---|
742 | |
---|
743 | def Xtest_urs2sts_ordering_exception(self): |
---|
744 | """Test that inconsistent lats and lons in ordering file are caught. |
---|
745 | """ |
---|
746 | |
---|
747 | tide=0 |
---|
748 | time_step_count = 3 |
---|
749 | time_step = 2 |
---|
750 | lat_long_points =[(-21.5,114.5),(-21,114.5),(-21.5,115), (-21.,115.)] |
---|
751 | n=len(lat_long_points) |
---|
752 | first_tstep=num.ones(n,num.int) |
---|
753 | first_tstep[0]+=1 |
---|
754 | first_tstep[2]+=1 |
---|
755 | last_tstep=(time_step_count)*num.ones(n,num.int) |
---|
756 | last_tstep[0]-=1 |
---|
757 | |
---|
758 | gauge_depth=20*num.ones(n,num.float) |
---|
759 | ha=2*num.ones((n,time_step_count),num.float) |
---|
760 | ha[0]=num.arange(0,time_step_count) |
---|
761 | ha[1]=num.arange(time_step_count,2*time_step_count) |
---|
762 | ha[2]=num.arange(2*time_step_count,3*time_step_count) |
---|
763 | ha[3]=num.arange(3*time_step_count,4*time_step_count) |
---|
764 | ua=5*num.ones((n,time_step_count),num.float) |
---|
765 | va=-10*num.ones((n,time_step_count),num.float) |
---|
766 | |
---|
767 | # Create two identical mux files to be combined by urs2sts |
---|
768 | base_nameI, filesI = self.write_mux2(lat_long_points, |
---|
769 | time_step_count, time_step, |
---|
770 | first_tstep, last_tstep, |
---|
771 | depth=gauge_depth, |
---|
772 | ha=ha, |
---|
773 | ua=ua, |
---|
774 | va=va) |
---|
775 | |
---|
776 | base_nameII, filesII = self.write_mux2(lat_long_points, |
---|
777 | time_step_count, time_step, |
---|
778 | first_tstep, last_tstep, |
---|
779 | depth=gauge_depth, |
---|
780 | ha=ha, |
---|
781 | ua=ua, |
---|
782 | va=va) |
---|
783 | |
---|
784 | |
---|
785 | # Create ordering file |
---|
786 | permutation = [3,0,2] |
---|
787 | |
---|
788 | # Do it wrongly and check that exception is being raised |
---|
789 | _, ordering_filename = tempfile.mkstemp('') |
---|
790 | order_fid = open(ordering_filename, 'w') |
---|
791 | order_fid.write('index, longitude, latitude\n') |
---|
792 | for index in permutation: |
---|
793 | order_fid.write('%d, %f, %f\n' %(index, |
---|
794 | lat_long_points[index][0], |
---|
795 | lat_long_points[index][1])) |
---|
796 | order_fid.close() |
---|
797 | |
---|
798 | try: |
---|
799 | urs2sts([base_nameI, base_nameII], |
---|
800 | basename_out=base_nameI, |
---|
801 | ordering_filename=ordering_filename, |
---|
802 | weights=[1.0, 1.0], |
---|
803 | mean_stage=tide, |
---|
804 | verbose=False) |
---|
805 | os.remove(ordering_filename) |
---|
806 | except: |
---|
807 | pass |
---|
808 | else: |
---|
809 | msg = 'Should have caught wrong lat longs' |
---|
810 | raise Exception, msg |
---|
811 | |
---|
812 | |
---|
813 | self.delete_mux(filesI) |
---|
814 | self.delete_mux(filesII) |
---|
815 | |
---|
816 | |
---|
817 | |
---|
818 | |
---|
819 | def test_urs2sts_ordering_different_sources(self): |
---|
820 | """Test multiple sources with ordering file, different source files and weights. |
---|
821 | This test also has more variable values than the previous ones |
---|
822 | """ |
---|
823 | |
---|
824 | tide = 1.5 |
---|
825 | time_step_count = 10 |
---|
826 | time_step = 0.2 |
---|
827 | |
---|
828 | times_ref = num.arange(0, time_step_count*time_step, time_step) |
---|
829 | #print 'time vector', times_ref |
---|
830 | |
---|
831 | lat_long_points = [(-21.5,114.5), (-21,114.5), (-21.5,115), (-21.,115.), (-22., 117.)] |
---|
832 | n = len(lat_long_points) |
---|
833 | |
---|
834 | # Create non-trivial weights |
---|
835 | #weights = [0.8, 1.5] # OK |
---|
836 | #weights = [0.8, 10.5] # Fail (up to allclose tolerance) |
---|
837 | #weights = [10.5, 10.5] # OK |
---|
838 | #weights = [0.0, 10.5] # OK |
---|
839 | #weights = [0.8, 0.] # OK |
---|
840 | #weights = [8, 0.1] # OK |
---|
841 | #weights = [0.8, 10.0] # OK |
---|
842 | #weights = [0.8, 10.6] # OK |
---|
843 | weights = [3.8, 7.6] # OK |
---|
844 | #weights = [0.5, 0.5] # OK |
---|
845 | #weights = [2., 2.] # OK |
---|
846 | #weights = [0.0, 0.5] # OK |
---|
847 | #weights = [1.0, 1.0] # OK |
---|
848 | |
---|
849 | |
---|
850 | # Create different timeseries starting and ending at different times |
---|
851 | first_tstep=num.ones(n,num.int) |
---|
852 | first_tstep[0]+=2 # Point 0 starts at 2 |
---|
853 | first_tstep[1]+=4 # Point 1 starts at 4 |
---|
854 | first_tstep[2]+=3 # Point 2 starts at 3 |
---|
855 | |
---|
856 | last_tstep=(time_step_count)*num.ones(n,num.int) |
---|
857 | last_tstep[0]-=1 # Point 0 ends 1 step early |
---|
858 | last_tstep[1]-=2 # Point 1 ends 2 steps early |
---|
859 | last_tstep[4]-=3 # Point 4 ends 3 steps early |
---|
860 | |
---|
861 | #print |
---|
862 | #print 'time_step_count', time_step_count |
---|
863 | #print 'time_step', time_step |
---|
864 | #print 'first_tstep', first_tstep |
---|
865 | #print 'last_tstep', last_tstep |
---|
866 | |
---|
867 | |
---|
868 | # Create varying elevation data (positive values for seafloor) |
---|
869 | gauge_depth=20*num.ones(n,num.float) |
---|
870 | for i in range(n): |
---|
871 | gauge_depth[i] += i**2 |
---|
872 | |
---|
873 | #print 'gauge_depth', gauge_depth |
---|
874 | |
---|
875 | # Create data to be written to first mux file |
---|
876 | ha0=2*num.ones((n,time_step_count),num.float) |
---|
877 | ha0[0]=num.arange(0,time_step_count) |
---|
878 | ha0[1]=num.arange(time_step_count,2*time_step_count) |
---|
879 | ha0[2]=num.arange(2*time_step_count,3*time_step_count) |
---|
880 | ha0[3]=num.arange(3*time_step_count,4*time_step_count) |
---|
881 | ua0=5*num.ones((n,time_step_count),num.float) |
---|
882 | va0=-10*num.ones((n,time_step_count),num.float) |
---|
883 | |
---|
884 | # Ensure data used to write mux file to be zero when gauges are |
---|
885 | # not recording |
---|
886 | for i in range(n): |
---|
887 | # For each point |
---|
888 | |
---|
889 | for j in range(0, first_tstep[i]-1) + range(last_tstep[i], time_step_count): |
---|
890 | # For timesteps before and after recording range |
---|
891 | ha0[i][j] = ua0[i][j] = va0[i][j] = 0.0 |
---|
892 | |
---|
893 | |
---|
894 | |
---|
895 | #print |
---|
896 | #print 'using varying start and end time' |
---|
897 | #print 'ha0', ha0[4] |
---|
898 | #print 'ua0', ua0 |
---|
899 | #print 'va0', va0 |
---|
900 | |
---|
901 | # Write first mux file to be combined by urs2sts |
---|
902 | base_nameI, filesI = self.write_mux2(lat_long_points, |
---|
903 | time_step_count, time_step, |
---|
904 | first_tstep, last_tstep, |
---|
905 | depth=gauge_depth, |
---|
906 | ha=ha0, |
---|
907 | ua=ua0, |
---|
908 | va=va0) |
---|
909 | |
---|
910 | |
---|
911 | |
---|
912 | # Create data to be written to second mux file |
---|
913 | ha1=num.ones((n,time_step_count),num.float) |
---|
914 | ha1[0]=num.sin(times_ref) |
---|
915 | ha1[1]=2*num.sin(times_ref - 3) |
---|
916 | ha1[2]=5*num.sin(4*times_ref) |
---|
917 | ha1[3]=num.sin(times_ref) |
---|
918 | ha1[4]=num.sin(2*times_ref-0.7) |
---|
919 | |
---|
920 | ua1=num.zeros((n,time_step_count),num.float) |
---|
921 | ua1[0]=3*num.cos(times_ref) |
---|
922 | ua1[1]=2*num.sin(times_ref-0.7) |
---|
923 | ua1[2]=num.arange(3*time_step_count,4*time_step_count) |
---|
924 | ua1[4]=2*num.ones(time_step_count) |
---|
925 | |
---|
926 | va1=num.zeros((n,time_step_count),num.float) |
---|
927 | va1[0]=2*num.cos(times_ref-0.87) |
---|
928 | va1[1]=3*num.ones(time_step_count, num.int) #array default# |
---|
929 | va1[3]=2*num.sin(times_ref-0.71) |
---|
930 | |
---|
931 | |
---|
932 | # Ensure data used to write mux file to be zero when gauges are |
---|
933 | # not recording |
---|
934 | for i in range(n): |
---|
935 | # For each point |
---|
936 | |
---|
937 | for j in range(0, first_tstep[i]-1) + range(last_tstep[i], time_step_count): |
---|
938 | # For timesteps before and after recording range |
---|
939 | ha1[i][j] = ua1[i][j] = va1[i][j] = 0.0 |
---|
940 | |
---|
941 | |
---|
942 | #print |
---|
943 | #print 'using varying start and end time' |
---|
944 | #print 'ha1', ha1[4] |
---|
945 | #print 'ua1', ua1 |
---|
946 | #print 'va1', va1 |
---|
947 | |
---|
948 | |
---|
949 | # Write second mux file to be combined by urs2sts |
---|
950 | base_nameII, filesII = self.write_mux2(lat_long_points, |
---|
951 | time_step_count, time_step, |
---|
952 | first_tstep, last_tstep, |
---|
953 | depth=gauge_depth, |
---|
954 | ha=ha1, |
---|
955 | ua=ua1, |
---|
956 | va=va1) |
---|
957 | |
---|
958 | |
---|
959 | # Create ordering file |
---|
960 | permutation = [4,0,2] |
---|
961 | |
---|
962 | _, ordering_filename = tempfile.mkstemp('') |
---|
963 | order_fid = open(ordering_filename, 'w') |
---|
964 | order_fid.write('index, longitude, latitude\n') |
---|
965 | for index in permutation: |
---|
966 | order_fid.write('%d, %f, %f\n' %(index, |
---|
967 | lat_long_points[index][1], |
---|
968 | lat_long_points[index][0])) |
---|
969 | order_fid.close() |
---|
970 | |
---|
971 | |
---|
972 | |
---|
973 | #------------------------------------------------------------ |
---|
974 | # Now read the mux files one by one without weights and test |
---|
975 | |
---|
976 | # Call urs2sts with mux file #0 |
---|
977 | urs2sts([base_nameI], |
---|
978 | basename_out=base_nameI, |
---|
979 | ordering_filename=ordering_filename, |
---|
980 | mean_stage=tide, |
---|
981 | verbose=False) |
---|
982 | |
---|
983 | # Now read the sts file and check that values have been stored correctly. |
---|
984 | sts_file = base_nameI + '.sts' |
---|
985 | fid = NetCDFFile(sts_file) |
---|
986 | |
---|
987 | |
---|
988 | # Check that original indices have been stored |
---|
989 | stored_permutation = fid.variables['permutation'][:] |
---|
990 | msg = 'Permutation was not stored correctly. I got ' |
---|
991 | msg += str(stored_permutation) |
---|
992 | assert num.allclose(stored_permutation, permutation), msg |
---|
993 | |
---|
994 | |
---|
995 | |
---|
996 | |
---|
997 | # Make x and y absolute |
---|
998 | x = fid.variables['x'][:] |
---|
999 | y = fid.variables['y'][:] |
---|
1000 | |
---|
1001 | geo_reference = Geo_reference(NetCDFObject=fid) |
---|
1002 | points = geo_reference.get_absolute(map(None, x, y)) |
---|
1003 | points = ensure_numeric(points) |
---|
1004 | |
---|
1005 | x = points[:,0] |
---|
1006 | y = points[:,1] |
---|
1007 | |
---|
1008 | for i, index in enumerate(permutation): |
---|
1009 | # Check that STS points are stored in the correct order |
---|
1010 | |
---|
1011 | # Work out the UTM coordinates sts point i |
---|
1012 | zone, e, n = redfearn(lat_long_points[index][0], |
---|
1013 | lat_long_points[index][1]) |
---|
1014 | |
---|
1015 | #print i, [x[i],y[i]], [e,n] |
---|
1016 | assert num.allclose([x[i],y[i]], [e,n]) |
---|
1017 | |
---|
1018 | |
---|
1019 | # Check the time vector |
---|
1020 | times = fid.variables['time'][:] |
---|
1021 | assert num.allclose(ensure_numeric(times), |
---|
1022 | ensure_numeric(times_ref)) |
---|
1023 | |
---|
1024 | |
---|
1025 | # Check sts values for mux #0 |
---|
1026 | stage0 = fid.variables['stage'][:].copy() |
---|
1027 | xmomentum0 = fid.variables['xmomentum'][:].copy() |
---|
1028 | ymomentum0 = fid.variables['ymomentum'][:].copy() |
---|
1029 | elevation0 = fid.variables['elevation'][:].copy() |
---|
1030 | |
---|
1031 | |
---|
1032 | #print 'stage', stage0 |
---|
1033 | #print 'xmomentum', xmomentum0 |
---|
1034 | #print 'ymomentum', ymomentum0 |
---|
1035 | #print 'elevation', elevation0 |
---|
1036 | |
---|
1037 | # The quantities stored in the .sts file should be the weighted sum of the |
---|
1038 | # quantities written to the mux2 files subject to the permutation vector. |
---|
1039 | |
---|
1040 | ha_ref = num.take(ha0, permutation, axis=0) |
---|
1041 | ua_ref = num.take(ua0, permutation, axis=0) |
---|
1042 | va_ref = num.take(va0, permutation, axis=0) |
---|
1043 | |
---|
1044 | gauge_depth_ref = num.take(gauge_depth, permutation, axis=0) |
---|
1045 | |
---|
1046 | assert num.allclose(num.transpose(ha_ref)+tide, stage0) # Meters |
---|
1047 | |
---|
1048 | |
---|
1049 | |
---|
1050 | #Check the momentums - ua |
---|
1051 | #momentum = velocity*(stage-elevation) |
---|
1052 | # elevation = - depth |
---|
1053 | #momentum = velocity_ua *(stage+depth) |
---|
1054 | |
---|
1055 | depth_ref = num.zeros((len(permutation), time_step_count), num.float) |
---|
1056 | for i in range(len(permutation)): |
---|
1057 | depth_ref[i]=gauge_depth_ref[i]+tide+ha_ref[i] |
---|
1058 | |
---|
1059 | |
---|
1060 | # The xmomentum stored in the .sts file should be the sum of the ua |
---|
1061 | # in the two mux2 files multiplied by the depth. |
---|
1062 | assert num.allclose(num.transpose(ua_ref*depth_ref), xmomentum0) |
---|
1063 | |
---|
1064 | #Check the momentums - va |
---|
1065 | #momentum = velocity*(stage-elevation) |
---|
1066 | # elevation = - depth |
---|
1067 | #momentum = velocity_va *(stage+depth) |
---|
1068 | |
---|
1069 | # The ymomentum stored in the .sts file should be the sum of the va |
---|
1070 | # in the two mux2 files multiplied by the depth. |
---|
1071 | |
---|
1072 | |
---|
1073 | #print transpose(va_ref*depth_ref) |
---|
1074 | #print ymomentum |
---|
1075 | assert num.allclose(num.transpose(va_ref*depth_ref), ymomentum0) |
---|
1076 | |
---|
1077 | # check the elevation values. |
---|
1078 | # -ve since urs measures depth, sww meshers height, |
---|
1079 | assert num.allclose(-gauge_depth_ref, elevation0) |
---|
1080 | |
---|
1081 | fid.close() |
---|
1082 | os.remove(sts_file) |
---|
1083 | |
---|
1084 | |
---|
1085 | |
---|
1086 | |
---|
1087 | # Call urs2sts with mux file #1 |
---|
1088 | urs2sts([base_nameII], |
---|
1089 | basename_out=base_nameI, |
---|
1090 | ordering_filename=ordering_filename, |
---|
1091 | mean_stage=tide, |
---|
1092 | verbose=False) |
---|
1093 | |
---|
1094 | # Now read the sts file and check that values have been stored correctly. |
---|
1095 | sts_file = base_nameI + '.sts' |
---|
1096 | fid = NetCDFFile(sts_file) |
---|
1097 | |
---|
1098 | |
---|
1099 | # Check that original indices have been stored |
---|
1100 | stored_permutation = fid.variables['permutation'][:] |
---|
1101 | msg = 'Permutation was not stored correctly. I got ' |
---|
1102 | msg += str(stored_permutation) |
---|
1103 | assert num.allclose(stored_permutation, permutation), msg |
---|
1104 | |
---|
1105 | # Make x and y absolute |
---|
1106 | x = fid.variables['x'][:] |
---|
1107 | y = fid.variables['y'][:] |
---|
1108 | |
---|
1109 | geo_reference = Geo_reference(NetCDFObject=fid) |
---|
1110 | points = geo_reference.get_absolute(map(None, x, y)) |
---|
1111 | points = ensure_numeric(points) |
---|
1112 | |
---|
1113 | x = points[:,0] |
---|
1114 | y = points[:,1] |
---|
1115 | |
---|
1116 | for i, index in enumerate(permutation): |
---|
1117 | # Check that STS points are stored in the correct order |
---|
1118 | |
---|
1119 | # Work out the UTM coordinates sts point i |
---|
1120 | zone, e, n = redfearn(lat_long_points[index][0], |
---|
1121 | lat_long_points[index][1]) |
---|
1122 | |
---|
1123 | #print i, [x[i],y[i]], [e,n] |
---|
1124 | assert num.allclose([x[i],y[i]], [e,n]) |
---|
1125 | |
---|
1126 | |
---|
1127 | # Check the time vector |
---|
1128 | times = fid.variables['time'][:] |
---|
1129 | assert num.allclose(ensure_numeric(times), |
---|
1130 | ensure_numeric(times_ref)) |
---|
1131 | |
---|
1132 | |
---|
1133 | # Check sts values for mux #1 |
---|
1134 | stage1 = fid.variables['stage'][:].copy() |
---|
1135 | xmomentum1 = fid.variables['xmomentum'][:].copy() |
---|
1136 | ymomentum1 = fid.variables['ymomentum'][:].copy() |
---|
1137 | elevation1 = fid.variables['elevation'][:].copy() |
---|
1138 | |
---|
1139 | |
---|
1140 | #print 'stage', stage1 |
---|
1141 | #print 'xmomentum', xmomentum1 |
---|
1142 | #print 'ymomentum', ymomentum1 |
---|
1143 | #print 'elevation', elevation1 |
---|
1144 | |
---|
1145 | # The quantities stored in the .sts file should be the weighted sum of the |
---|
1146 | # quantities written to the mux2 files subject to the permutation vector. |
---|
1147 | |
---|
1148 | ha_ref = num.take(ha1, permutation, axis=0) |
---|
1149 | ua_ref = num.take(ua1, permutation, axis=0) |
---|
1150 | va_ref = num.take(va1, permutation, axis=0) |
---|
1151 | |
---|
1152 | gauge_depth_ref = num.take(gauge_depth, permutation, axis=0) |
---|
1153 | |
---|
1154 | |
---|
1155 | #print |
---|
1156 | #print stage1 |
---|
1157 | #print transpose(ha_ref)+tide - stage1 |
---|
1158 | |
---|
1159 | |
---|
1160 | assert num.allclose(num.transpose(ha_ref)+tide, stage1) # Meters |
---|
1161 | #import sys; sys.exit() |
---|
1162 | |
---|
1163 | #Check the momentums - ua |
---|
1164 | #momentum = velocity*(stage-elevation) |
---|
1165 | # elevation = - depth |
---|
1166 | #momentum = velocity_ua *(stage+depth) |
---|
1167 | |
---|
1168 | depth_ref = num.zeros((len(permutation), time_step_count), num.float) |
---|
1169 | for i in range(len(permutation)): |
---|
1170 | depth_ref[i]=gauge_depth_ref[i]+tide+ha_ref[i] |
---|
1171 | |
---|
1172 | |
---|
1173 | # The xmomentum stored in the .sts file should be the sum of the ua |
---|
1174 | # in the two mux2 files multiplied by the depth. |
---|
1175 | assert num.allclose(num.transpose(ua_ref*depth_ref), xmomentum1) |
---|
1176 | |
---|
1177 | #Check the momentums - va |
---|
1178 | #momentum = velocity*(stage-elevation) |
---|
1179 | # elevation = - depth |
---|
1180 | #momentum = velocity_va *(stage+depth) |
---|
1181 | |
---|
1182 | # The ymomentum stored in the .sts file should be the sum of the va |
---|
1183 | # in the two mux2 files multiplied by the depth. |
---|
1184 | |
---|
1185 | |
---|
1186 | #print transpose(va_ref*depth_ref) |
---|
1187 | #print ymomentum |
---|
1188 | assert num.allclose(num.transpose(va_ref*depth_ref), ymomentum1) |
---|
1189 | |
---|
1190 | # check the elevation values. |
---|
1191 | # -ve since urs measures depth, sww meshers height, |
---|
1192 | assert num.allclose(-gauge_depth_ref, elevation1) |
---|
1193 | |
---|
1194 | fid.close() |
---|
1195 | os.remove(sts_file) |
---|
1196 | |
---|
1197 | #---------------------- |
---|
1198 | # Then read the mux files together and test |
---|
1199 | |
---|
1200 | |
---|
1201 | # Call urs2sts with multiple mux files |
---|
1202 | urs2sts([base_nameI, base_nameII], |
---|
1203 | basename_out=base_nameI, |
---|
1204 | ordering_filename=ordering_filename, |
---|
1205 | weights=weights, |
---|
1206 | mean_stage=tide, |
---|
1207 | verbose=False) |
---|
1208 | |
---|
1209 | # Now read the sts file and check that values have been stored correctly. |
---|
1210 | sts_file = base_nameI + '.sts' |
---|
1211 | fid = NetCDFFile(sts_file) |
---|
1212 | |
---|
1213 | # Make x and y absolute |
---|
1214 | x = fid.variables['x'][:] |
---|
1215 | y = fid.variables['y'][:] |
---|
1216 | |
---|
1217 | geo_reference = Geo_reference(NetCDFObject=fid) |
---|
1218 | points = geo_reference.get_absolute(map(None, x, y)) |
---|
1219 | points = ensure_numeric(points) |
---|
1220 | |
---|
1221 | x = points[:,0] |
---|
1222 | y = points[:,1] |
---|
1223 | |
---|
1224 | for i, index in enumerate(permutation): |
---|
1225 | # Check that STS points are stored in the correct order |
---|
1226 | |
---|
1227 | # Work out the UTM coordinates sts point i |
---|
1228 | zone, e, n = redfearn(lat_long_points[index][0], |
---|
1229 | lat_long_points[index][1]) |
---|
1230 | |
---|
1231 | #print i, [x[i],y[i]], [e,n] |
---|
1232 | assert num.allclose([x[i],y[i]], [e,n]) |
---|
1233 | |
---|
1234 | |
---|
1235 | # Check the time vector |
---|
1236 | times = fid.variables['time'][:] |
---|
1237 | assert num.allclose(ensure_numeric(times), |
---|
1238 | ensure_numeric(times_ref)) |
---|
1239 | |
---|
1240 | |
---|
1241 | # Check sts values |
---|
1242 | stage = fid.variables['stage'][:].copy() |
---|
1243 | xmomentum = fid.variables['xmomentum'][:].copy() |
---|
1244 | ymomentum = fid.variables['ymomentum'][:].copy() |
---|
1245 | elevation = fid.variables['elevation'][:].copy() |
---|
1246 | |
---|
1247 | |
---|
1248 | #print 'stage', stage |
---|
1249 | #print 'elevation', elevation |
---|
1250 | |
---|
1251 | # The quantities stored in the .sts file should be the weighted sum of the |
---|
1252 | # quantities written to the mux2 files subject to the permutation vector. |
---|
1253 | |
---|
1254 | ha_ref = (weights[0]*num.take(ha0, permutation, axis=0) |
---|
1255 | + weights[1]*num.take(ha1, permutation, axis=0)) |
---|
1256 | ua_ref = (weights[0]*num.take(ua0, permutation, axis=0) |
---|
1257 | + weights[1]*num.take(ua1, permutation, axis=0)) |
---|
1258 | va_ref = (weights[0]*num.take(va0, permutation, axis=0) |
---|
1259 | + weights[1]*num.take(va1, permutation, axis=0)) |
---|
1260 | |
---|
1261 | gauge_depth_ref = num.take(gauge_depth, permutation, axis=0) |
---|
1262 | |
---|
1263 | |
---|
1264 | #print |
---|
1265 | #print stage |
---|
1266 | #print transpose(ha_ref)+tide - stage |
---|
1267 | |
---|
1268 | assert num.allclose(num.transpose(ha_ref)+tide, stage) # Meters |
---|
1269 | |
---|
1270 | #Check the momentums - ua |
---|
1271 | #momentum = velocity*(stage-elevation) |
---|
1272 | # elevation = - depth |
---|
1273 | #momentum = velocity_ua *(stage+depth) |
---|
1274 | |
---|
1275 | depth_ref = num.zeros((len(permutation), time_step_count), num.float) |
---|
1276 | for i in range(len(permutation)): |
---|
1277 | depth_ref[i]=gauge_depth_ref[i]+tide+ha_ref[i] |
---|
1278 | |
---|
1279 | |
---|
1280 | |
---|
1281 | |
---|
1282 | # The xmomentum stored in the .sts file should be the sum of the ua |
---|
1283 | # in the two mux2 files multiplied by the depth. |
---|
1284 | assert num.allclose(num.transpose(ua_ref*depth_ref), xmomentum) |
---|
1285 | |
---|
1286 | #Check the momentums - va |
---|
1287 | #momentum = velocity*(stage-elevation) |
---|
1288 | # elevation = - depth |
---|
1289 | #momentum = velocity_va *(stage+depth) |
---|
1290 | |
---|
1291 | # The ymomentum stored in the .sts file should be the sum of the va |
---|
1292 | # in the two mux2 files multiplied by the depth. |
---|
1293 | |
---|
1294 | |
---|
1295 | #print transpose(va_ref*depth_ref) |
---|
1296 | #print ymomentum |
---|
1297 | |
---|
1298 | assert num.allclose(num.transpose(va_ref*depth_ref), ymomentum) |
---|
1299 | |
---|
1300 | # check the elevation values. |
---|
1301 | # -ve since urs measures depth, sww meshers height, |
---|
1302 | assert num.allclose(-gauge_depth_ref, elevation) #Meters |
---|
1303 | |
---|
1304 | fid.close() |
---|
1305 | self.delete_mux(filesI) |
---|
1306 | self.delete_mux(filesII) |
---|
1307 | os.remove(sts_file) |
---|
1308 | |
---|
1309 | #--------------- |
---|
1310 | # "Manually" add the timeseries up with weights and test |
---|
1311 | # Tide is discounted from individual results and added back in |
---|
1312 | # |
---|
1313 | |
---|
1314 | stage_man = weights[0]*(stage0-tide) + weights[1]*(stage1-tide) + tide |
---|
1315 | assert num.allclose(stage_man, stage) |
---|
1316 | |
---|
1317 | |
---|
1318 | def test_file_boundary_stsI(self): |
---|
1319 | """test_file_boundary_stsI(self): |
---|
1320 | """ |
---|
1321 | |
---|
1322 | bounding_polygon=[[6.0,97.0],[6.01,97.0],[6.02,97.0],[6.02,97.02],[6.00,97.02]] |
---|
1323 | tide = 0.37 |
---|
1324 | time_step_count = 5 |
---|
1325 | time_step = 2 |
---|
1326 | lat_long_points =bounding_polygon[0:3] |
---|
1327 | n=len(lat_long_points) |
---|
1328 | first_tstep=num.ones(n,num.int) |
---|
1329 | last_tstep=(time_step_count)*num.ones(n,num.int) |
---|
1330 | |
---|
1331 | h = 20 |
---|
1332 | w = 2 |
---|
1333 | u = 10 |
---|
1334 | v = -10 |
---|
1335 | gauge_depth=h*num.ones(n,num.float) |
---|
1336 | ha=w*num.ones((n,time_step_count),num.float) |
---|
1337 | ua=u*num.ones((n,time_step_count),num.float) |
---|
1338 | va=v*num.ones((n,time_step_count),num.float) |
---|
1339 | base_name, files = self.write_mux2(lat_long_points, |
---|
1340 | time_step_count, time_step, |
---|
1341 | first_tstep, last_tstep, |
---|
1342 | depth=gauge_depth, |
---|
1343 | ha=ha, |
---|
1344 | ua=ua, |
---|
1345 | va=va) |
---|
1346 | |
---|
1347 | sts_file=base_name |
---|
1348 | urs2sts(base_name, |
---|
1349 | sts_file, |
---|
1350 | mean_stage=tide, |
---|
1351 | verbose=False) |
---|
1352 | self.delete_mux(files) |
---|
1353 | |
---|
1354 | #print 'start create mesh from regions' |
---|
1355 | for i in range(len(bounding_polygon)): |
---|
1356 | zone,bounding_polygon[i][0],bounding_polygon[i][1]=redfearn(bounding_polygon[i][0],bounding_polygon[i][1]) |
---|
1357 | extent_res=1000000 |
---|
1358 | meshname = 'urs_test_mesh' + '.tsh' |
---|
1359 | interior_regions=None |
---|
1360 | boundary_tags={'ocean': [0,1], 'otherocean': [2,3,4]} |
---|
1361 | create_mesh_from_regions(bounding_polygon, |
---|
1362 | boundary_tags=boundary_tags, |
---|
1363 | maximum_triangle_area=extent_res, |
---|
1364 | filename=meshname, |
---|
1365 | interior_regions=interior_regions, |
---|
1366 | verbose=False) |
---|
1367 | |
---|
1368 | domain_fbound = Domain(meshname) |
---|
1369 | domain_fbound.set_quantity('stage', tide) |
---|
1370 | Bf = File_boundary(sts_file+'.sts', |
---|
1371 | domain_fbound, |
---|
1372 | boundary_polygon=bounding_polygon) |
---|
1373 | Br = Reflective_boundary(domain_fbound) |
---|
1374 | Bd = Dirichlet_boundary([w+tide, u*(w+h+tide), v*(w+h+tide)]) |
---|
1375 | |
---|
1376 | domain_fbound.set_boundary({'ocean': Bf,'otherocean': Br}) |
---|
1377 | |
---|
1378 | # Check boundary object evaluates as it should |
---|
1379 | for i, ((vol_id, edge_id), B) in enumerate(domain_fbound.boundary_objects): |
---|
1380 | if B is Bf: |
---|
1381 | |
---|
1382 | qf = B.evaluate(vol_id, edge_id) # File boundary |
---|
1383 | qd = Bd.evaluate(vol_id, edge_id) # Dirichlet boundary |
---|
1384 | |
---|
1385 | assert num.allclose(qf, qd) |
---|
1386 | |
---|
1387 | |
---|
1388 | # Evolve |
---|
1389 | finaltime=time_step*(time_step_count-1) |
---|
1390 | yieldstep=time_step |
---|
1391 | temp_fbound=num.zeros(int(finaltime/yieldstep)+1,num.float) |
---|
1392 | |
---|
1393 | for i, t in enumerate(domain_fbound.evolve(yieldstep=yieldstep, |
---|
1394 | finaltime=finaltime, |
---|
1395 | skip_initial_step=False)): |
---|
1396 | |
---|
1397 | D = domain_fbound |
---|
1398 | temp_fbound[i]=D.quantities['stage'].centroid_values[2] |
---|
1399 | |
---|
1400 | # Check that file boundary object has populated |
---|
1401 | # boundary array correctly |
---|
1402 | # FIXME (Ole): Do this for the other tests too! |
---|
1403 | for j, val in enumerate(D.get_quantity('stage').boundary_values): |
---|
1404 | |
---|
1405 | (vol_id, edge_id), B = D.boundary_objects[j] |
---|
1406 | if isinstance(B, File_boundary): |
---|
1407 | #print j, val |
---|
1408 | assert num.allclose(val, w + tide) |
---|
1409 | |
---|
1410 | |
---|
1411 | |
---|
1412 | domain_drchlt = Domain(meshname) |
---|
1413 | domain_drchlt.set_quantity('stage', tide) |
---|
1414 | Br = Reflective_boundary(domain_drchlt) |
---|
1415 | |
---|
1416 | domain_drchlt.set_boundary({'ocean': Bd,'otherocean': Br}) |
---|
1417 | temp_drchlt=num.zeros(int(finaltime/yieldstep)+1,num.float) |
---|
1418 | |
---|
1419 | for i, t in enumerate(domain_drchlt.evolve(yieldstep=yieldstep, |
---|
1420 | finaltime=finaltime, |
---|
1421 | skip_initial_step=False)): |
---|
1422 | temp_drchlt[i]=domain_drchlt.quantities['stage'].centroid_values[2] |
---|
1423 | |
---|
1424 | #print domain_fbound.quantities['stage'].vertex_values |
---|
1425 | #print domain_drchlt.quantities['stage'].vertex_values |
---|
1426 | |
---|
1427 | assert num.allclose(temp_fbound,temp_drchlt) |
---|
1428 | |
---|
1429 | assert num.allclose(domain_fbound.quantities['stage'].vertex_values, |
---|
1430 | domain_drchlt.quantities['stage'].vertex_values) |
---|
1431 | |
---|
1432 | assert num.allclose(domain_fbound.quantities['xmomentum'].vertex_values, |
---|
1433 | domain_drchlt.quantities['xmomentum'].vertex_values) |
---|
1434 | |
---|
1435 | assert num.allclose(domain_fbound.quantities['ymomentum'].vertex_values, |
---|
1436 | domain_drchlt.quantities['ymomentum'].vertex_values) |
---|
1437 | |
---|
1438 | |
---|
1439 | os.remove(sts_file+'.sts') |
---|
1440 | os.remove(meshname) |
---|
1441 | |
---|
1442 | |
---|
1443 | def test_file_boundary_stsI_beyond_model_time(self): |
---|
1444 | """test_file_boundary_stsI(self): |
---|
1445 | |
---|
1446 | Test that file_boundary can work when model time |
---|
1447 | exceeds available data. |
---|
1448 | This is optional and requires the user to specify a default |
---|
1449 | boundary object. |
---|
1450 | """ |
---|
1451 | |
---|
1452 | # Don't do warnings in unit test |
---|
1453 | import warnings |
---|
1454 | warnings.simplefilter('ignore') |
---|
1455 | |
---|
1456 | bounding_polygon=[[6.0,97.0],[6.01,97.0],[6.02,97.0], |
---|
1457 | [6.02,97.02],[6.00,97.02]] |
---|
1458 | tide = 0.37 |
---|
1459 | time_step_count = 5 |
---|
1460 | time_step = 2 |
---|
1461 | lat_long_points = bounding_polygon[0:3] |
---|
1462 | n=len(lat_long_points) |
---|
1463 | first_tstep=num.ones(n,num.int) |
---|
1464 | last_tstep=(time_step_count)*num.ones(n,num.int) |
---|
1465 | |
---|
1466 | h = 20 |
---|
1467 | w = 2 |
---|
1468 | u = 10 |
---|
1469 | v = -10 |
---|
1470 | gauge_depth=h*num.ones(n,num.float) |
---|
1471 | ha=w*num.ones((n,time_step_count),num.float) |
---|
1472 | ua=u*num.ones((n,time_step_count),num.float) |
---|
1473 | va=v*num.ones((n,time_step_count),num.float) |
---|
1474 | base_name, files = self.write_mux2(lat_long_points, |
---|
1475 | time_step_count, time_step, |
---|
1476 | first_tstep, last_tstep, |
---|
1477 | depth=gauge_depth, |
---|
1478 | ha=ha, |
---|
1479 | ua=ua, |
---|
1480 | va=va) |
---|
1481 | |
---|
1482 | sts_file=base_name |
---|
1483 | urs2sts(base_name, |
---|
1484 | sts_file, |
---|
1485 | mean_stage=tide, |
---|
1486 | verbose=False) |
---|
1487 | self.delete_mux(files) |
---|
1488 | |
---|
1489 | #print 'start create mesh from regions' |
---|
1490 | for i in range(len(bounding_polygon)): |
---|
1491 | zone,\ |
---|
1492 | bounding_polygon[i][0],\ |
---|
1493 | bounding_polygon[i][1]=redfearn(bounding_polygon[i][0], |
---|
1494 | bounding_polygon[i][1]) |
---|
1495 | |
---|
1496 | extent_res=1000000 |
---|
1497 | meshname = 'urs_test_mesh' + '.tsh' |
---|
1498 | interior_regions=None |
---|
1499 | boundary_tags={'ocean': [0,1], 'otherocean': [2,3,4]} |
---|
1500 | create_mesh_from_regions(bounding_polygon, |
---|
1501 | boundary_tags=boundary_tags, |
---|
1502 | maximum_triangle_area=extent_res, |
---|
1503 | filename=meshname, |
---|
1504 | interior_regions=interior_regions, |
---|
1505 | verbose=False) |
---|
1506 | |
---|
1507 | domain_fbound = Domain(meshname) |
---|
1508 | domain_fbound.set_quantity('stage', tide) |
---|
1509 | |
---|
1510 | Br = Reflective_boundary(domain_fbound) |
---|
1511 | Bd = Dirichlet_boundary([w+tide, u*(w+h+tide), v*(w+h+tide)]) |
---|
1512 | Bf = File_boundary(sts_file+'.sts', |
---|
1513 | domain_fbound, |
---|
1514 | boundary_polygon=bounding_polygon, |
---|
1515 | default_boundary=Bd) # Condition to be used |
---|
1516 | # if model time exceeds |
---|
1517 | # available data |
---|
1518 | |
---|
1519 | domain_fbound.set_boundary({'ocean': Bf,'otherocean': Br}) |
---|
1520 | |
---|
1521 | # Check boundary object evaluates as it should |
---|
1522 | for i, ((vol_id, edge_id), B) in enumerate(domain_fbound.boundary_objects): |
---|
1523 | if B is Bf: |
---|
1524 | |
---|
1525 | qf = B.evaluate(vol_id, edge_id) # File boundary |
---|
1526 | qd = Bd.evaluate(vol_id, edge_id) # Dirichlet boundary |
---|
1527 | |
---|
1528 | assert num.allclose(qf, qd) |
---|
1529 | |
---|
1530 | |
---|
1531 | # Evolve |
---|
1532 | data_finaltime = time_step*(time_step_count-1) |
---|
1533 | finaltime = data_finaltime + 10 # Let model time exceed available data |
---|
1534 | yieldstep = time_step |
---|
1535 | temp_fbound=num.zeros(int(finaltime/yieldstep)+1, num.float) |
---|
1536 | |
---|
1537 | for i, t in enumerate(domain_fbound.evolve(yieldstep=yieldstep, |
---|
1538 | finaltime=finaltime, |
---|
1539 | skip_initial_step=False)): |
---|
1540 | |
---|
1541 | D = domain_fbound |
---|
1542 | temp_fbound[i]=D.quantities['stage'].centroid_values[2] |
---|
1543 | |
---|
1544 | # Check that file boundary object has populated |
---|
1545 | # boundary array correctly |
---|
1546 | # FIXME (Ole): Do this for the other tests too! |
---|
1547 | for j, val in enumerate(D.get_quantity('stage').boundary_values): |
---|
1548 | |
---|
1549 | (vol_id, edge_id), B = D.boundary_objects[j] |
---|
1550 | if isinstance(B, File_boundary): |
---|
1551 | #print j, val |
---|
1552 | assert num.allclose(val, w + tide) |
---|
1553 | |
---|
1554 | |
---|
1555 | domain_drchlt = Domain(meshname) |
---|
1556 | domain_drchlt.set_quantity('stage', tide) |
---|
1557 | Br = Reflective_boundary(domain_drchlt) |
---|
1558 | |
---|
1559 | domain_drchlt.set_boundary({'ocean': Bd,'otherocean': Br}) |
---|
1560 | temp_drchlt=num.zeros(int(finaltime/yieldstep)+1,num.float) |
---|
1561 | |
---|
1562 | for i, t in enumerate(domain_drchlt.evolve(yieldstep=yieldstep, |
---|
1563 | finaltime=finaltime, |
---|
1564 | skip_initial_step=False)): |
---|
1565 | temp_drchlt[i]=domain_drchlt.quantities['stage'].centroid_values[2] |
---|
1566 | |
---|
1567 | #print domain_fbound.quantities['stage'].vertex_values |
---|
1568 | #print domain_drchlt.quantities['stage'].vertex_values |
---|
1569 | |
---|
1570 | assert num.allclose(temp_fbound,temp_drchlt) |
---|
1571 | |
---|
1572 | assert num.allclose(domain_fbound.quantities['stage'].vertex_values, |
---|
1573 | domain_drchlt.quantities['stage'].vertex_values) |
---|
1574 | |
---|
1575 | assert num.allclose(domain_fbound.quantities['xmomentum'].vertex_values, |
---|
1576 | domain_drchlt.quantities['xmomentum'].vertex_values) |
---|
1577 | |
---|
1578 | assert num.allclose(domain_fbound.quantities['ymomentum'].vertex_values, |
---|
1579 | domain_drchlt.quantities['ymomentum'].vertex_values) |
---|
1580 | |
---|
1581 | os.remove(sts_file+'.sts') |
---|
1582 | os.remove(meshname) |
---|
1583 | |
---|
1584 | |
---|
1585 | def test_field_boundary_stsI_beyond_model_time(self): |
---|
1586 | """test_field_boundary(self): |
---|
1587 | |
---|
1588 | Test that field_boundary can work when model time |
---|
1589 | exceeds available data whilst adjusting mean_stage. |
---|
1590 | |
---|
1591 | """ |
---|
1592 | |
---|
1593 | # Don't do warnings in unit test |
---|
1594 | import warnings |
---|
1595 | warnings.simplefilter('ignore') |
---|
1596 | |
---|
1597 | bounding_polygon=[[6.0,97.0],[6.01,97.0],[6.02,97.0], |
---|
1598 | [6.02,97.02],[6.00,97.02]] |
---|
1599 | tide = 0.37 |
---|
1600 | time_step_count = 5 |
---|
1601 | time_step = 2 |
---|
1602 | lat_long_points = bounding_polygon[0:3] |
---|
1603 | n=len(lat_long_points) |
---|
1604 | first_tstep=num.ones(n,num.int) |
---|
1605 | last_tstep=(time_step_count)*num.ones(n,num.int) |
---|
1606 | |
---|
1607 | h = 20 |
---|
1608 | w = 2 |
---|
1609 | u = 10 |
---|
1610 | v = -10 |
---|
1611 | gauge_depth=h*num.ones(n,num.float) |
---|
1612 | ha=w*num.ones((n,time_step_count),num.float) |
---|
1613 | ua=u*num.ones((n,time_step_count),num.float) |
---|
1614 | va=v*num.ones((n,time_step_count),num.float) |
---|
1615 | base_name, files = self.write_mux2(lat_long_points, |
---|
1616 | time_step_count, time_step, |
---|
1617 | first_tstep, last_tstep, |
---|
1618 | depth=gauge_depth, |
---|
1619 | ha=ha, |
---|
1620 | ua=ua, |
---|
1621 | va=va) |
---|
1622 | |
---|
1623 | sts_file=base_name |
---|
1624 | urs2sts(base_name, |
---|
1625 | sts_file, |
---|
1626 | mean_stage=0.0, # Deliberately let Field_boundary do the adjustment |
---|
1627 | verbose=False) |
---|
1628 | self.delete_mux(files) |
---|
1629 | |
---|
1630 | #print 'start create mesh from regions' |
---|
1631 | for i in range(len(bounding_polygon)): |
---|
1632 | zone,\ |
---|
1633 | bounding_polygon[i][0],\ |
---|
1634 | bounding_polygon[i][1]=redfearn(bounding_polygon[i][0], |
---|
1635 | bounding_polygon[i][1]) |
---|
1636 | |
---|
1637 | extent_res=1000000 |
---|
1638 | meshname = 'urs_test_mesh' + '.tsh' |
---|
1639 | interior_regions=None |
---|
1640 | boundary_tags={'ocean': [0,1], 'otherocean': [2,3,4]} |
---|
1641 | create_mesh_from_regions(bounding_polygon, |
---|
1642 | boundary_tags=boundary_tags, |
---|
1643 | maximum_triangle_area=extent_res, |
---|
1644 | filename=meshname, |
---|
1645 | interior_regions=interior_regions, |
---|
1646 | verbose=False) |
---|
1647 | |
---|
1648 | domain_fbound = Domain(meshname) |
---|
1649 | domain_fbound.set_quantity('stage', tide) |
---|
1650 | |
---|
1651 | Br = Reflective_boundary(domain_fbound) |
---|
1652 | Bd = Dirichlet_boundary([w+tide, u*(w+h), v*(w+h)]) |
---|
1653 | Bdefault = Dirichlet_boundary([w, u*(w+h), v*(w+h)]) |
---|
1654 | |
---|
1655 | Bf = Field_boundary(sts_file+'.sts', |
---|
1656 | domain_fbound, |
---|
1657 | mean_stage=tide, # Field boundary to adjust for tide |
---|
1658 | boundary_polygon=bounding_polygon, |
---|
1659 | default_boundary=Bdefault) # Condition to be used |
---|
1660 | # if model time exceeds |
---|
1661 | # available data |
---|
1662 | |
---|
1663 | domain_fbound.set_boundary({'ocean': Bf,'otherocean': Br}) |
---|
1664 | |
---|
1665 | # Check boundary object evaluates as it should |
---|
1666 | for i, ((vol_id, edge_id), B) in enumerate(domain_fbound.boundary_objects): |
---|
1667 | if B is Bf: |
---|
1668 | |
---|
1669 | qf = B.evaluate(vol_id, edge_id) # Field boundary |
---|
1670 | qd = Bd.evaluate(vol_id, edge_id) # Dirichlet boundary |
---|
1671 | |
---|
1672 | msg = 'Got %s, should have been %s' %(qf, qd) |
---|
1673 | assert num.allclose(qf, qd), msg |
---|
1674 | |
---|
1675 | # Evolve |
---|
1676 | data_finaltime = time_step*(time_step_count-1) |
---|
1677 | finaltime = data_finaltime + 10 # Let model time exceed available data |
---|
1678 | yieldstep = time_step |
---|
1679 | temp_fbound=num.zeros(int(finaltime/yieldstep)+1, num.float) |
---|
1680 | |
---|
1681 | for i, t in enumerate(domain_fbound.evolve(yieldstep=yieldstep, |
---|
1682 | finaltime=finaltime, |
---|
1683 | skip_initial_step=False)): |
---|
1684 | |
---|
1685 | D = domain_fbound |
---|
1686 | temp_fbound[i]=D.quantities['stage'].centroid_values[2] |
---|
1687 | |
---|
1688 | # Check that file boundary object has populated |
---|
1689 | # boundary array correctly |
---|
1690 | # FIXME (Ole): Do this for the other tests too! |
---|
1691 | for j, val in enumerate(D.get_quantity('stage').boundary_values): |
---|
1692 | |
---|
1693 | (vol_id, edge_id), B = D.boundary_objects[j] |
---|
1694 | if isinstance(B, Field_boundary): |
---|
1695 | msg = 'Got %f should have been %f' %(val, w+tide) |
---|
1696 | assert num.allclose(val, w + tide), msg |
---|
1697 | |
---|
1698 | |
---|
1699 | def test_file_boundary_stsII(self): |
---|
1700 | """test_file_boundary_stsII(self): |
---|
1701 | |
---|
1702 | mux2 file has points not included in boundary |
---|
1703 | mux2 gauges are not stored with the same order as they are |
---|
1704 | found in bounding_polygon. This does not matter as long as bounding |
---|
1705 | polygon passed to file_function contains the mux2 points needed (in |
---|
1706 | the correct order). |
---|
1707 | """ |
---|
1708 | |
---|
1709 | bounding_polygon=[[6.01,97.0],[6.02,97.0],[6.02,97.02],[6.00,97.02],[6.0,97.0]] |
---|
1710 | tide = -2.20 |
---|
1711 | time_step_count = 20 |
---|
1712 | time_step = 2 |
---|
1713 | lat_long_points=bounding_polygon[0:2] |
---|
1714 | lat_long_points.insert(0,bounding_polygon[len(bounding_polygon)-1]) |
---|
1715 | lat_long_points.insert(0,[6.0,97.01]) |
---|
1716 | n=len(lat_long_points) |
---|
1717 | first_tstep=num.ones(n,num.int) |
---|
1718 | last_tstep=(time_step_count)*num.ones(n,num.int) |
---|
1719 | gauge_depth=20*num.ones(n,num.float) |
---|
1720 | ha=2*num.ones((n,time_step_count),num.float) |
---|
1721 | ua=10*num.ones((n,time_step_count),num.float) |
---|
1722 | va=-10*num.ones((n,time_step_count),num.float) |
---|
1723 | base_name, files = self.write_mux2(lat_long_points, |
---|
1724 | time_step_count, |
---|
1725 | time_step, |
---|
1726 | first_tstep, |
---|
1727 | last_tstep, |
---|
1728 | depth=gauge_depth, |
---|
1729 | ha=ha, |
---|
1730 | ua=ua, |
---|
1731 | va=va) |
---|
1732 | |
---|
1733 | sts_file=base_name |
---|
1734 | urs2sts(base_name,sts_file,mean_stage=tide,verbose=False) |
---|
1735 | self.delete_mux(files) |
---|
1736 | |
---|
1737 | #print 'start create mesh from regions' |
---|
1738 | for i in range(len(bounding_polygon)): |
---|
1739 | zone,\ |
---|
1740 | bounding_polygon[i][0],\ |
---|
1741 | bounding_polygon[i][1]=redfearn(bounding_polygon[i][0], |
---|
1742 | bounding_polygon[i][1]) |
---|
1743 | |
---|
1744 | extent_res=1000000 |
---|
1745 | meshname = 'urs_test_mesh' + '.tsh' |
---|
1746 | interior_regions=None |
---|
1747 | boundary_tags={'ocean': [0,1], 'otherocean': [2,3,4]} |
---|
1748 | # have to change boundary tags from last example because now bounding |
---|
1749 | # polygon starts in different place. |
---|
1750 | create_mesh_from_regions(bounding_polygon,boundary_tags=boundary_tags, |
---|
1751 | maximum_triangle_area=extent_res,filename=meshname, |
---|
1752 | interior_regions=interior_regions,verbose=False) |
---|
1753 | |
---|
1754 | domain_fbound = Domain(meshname) |
---|
1755 | domain_fbound.set_quantity('stage', tide) |
---|
1756 | Bf = File_boundary(sts_file+'.sts', |
---|
1757 | domain_fbound, |
---|
1758 | boundary_polygon=bounding_polygon) |
---|
1759 | Br = Reflective_boundary(domain_fbound) |
---|
1760 | |
---|
1761 | domain_fbound.set_boundary({'ocean': Bf,'otherocean': Br}) |
---|
1762 | finaltime=time_step*(time_step_count-1) |
---|
1763 | yieldstep=time_step |
---|
1764 | temp_fbound=num.zeros(int(finaltime/yieldstep)+1,num.float) |
---|
1765 | |
---|
1766 | for i, t in enumerate(domain_fbound.evolve(yieldstep=yieldstep, |
---|
1767 | finaltime=finaltime, |
---|
1768 | skip_initial_step = False)): |
---|
1769 | temp_fbound[i]=domain_fbound.quantities['stage'].centroid_values[2] |
---|
1770 | |
---|
1771 | domain_drchlt = Domain(meshname) |
---|
1772 | domain_drchlt.set_quantity('stage', tide) |
---|
1773 | Br = Reflective_boundary(domain_drchlt) |
---|
1774 | Bd = Dirichlet_boundary([2.0+tide,220+10*tide,-220-10*tide]) |
---|
1775 | domain_drchlt.set_boundary({'ocean': Bd,'otherocean': Br}) |
---|
1776 | temp_drchlt=num.zeros(int(finaltime/yieldstep)+1,num.float) |
---|
1777 | |
---|
1778 | for i, t in enumerate(domain_drchlt.evolve(yieldstep=yieldstep, |
---|
1779 | finaltime=finaltime, |
---|
1780 | skip_initial_step = False)): |
---|
1781 | temp_drchlt[i]=domain_drchlt.quantities['stage'].centroid_values[2] |
---|
1782 | |
---|
1783 | |
---|
1784 | assert num.allclose(temp_fbound,temp_drchlt) |
---|
1785 | |
---|
1786 | #print domain_fbound.quantities['stage'].vertex_values |
---|
1787 | #print domain_drchlt.quantities['stage'].vertex_values |
---|
1788 | |
---|
1789 | |
---|
1790 | assert num.allclose(domain_fbound.quantities['stage'].vertex_values, |
---|
1791 | domain_drchlt.quantities['stage'].vertex_values) |
---|
1792 | |
---|
1793 | assert num.allclose(domain_fbound.quantities['xmomentum'].vertex_values, |
---|
1794 | domain_drchlt.quantities['xmomentum'].vertex_values) |
---|
1795 | |
---|
1796 | assert num.allclose(domain_fbound.quantities['ymomentum'].vertex_values, |
---|
1797 | domain_drchlt.quantities['ymomentum'].vertex_values) |
---|
1798 | |
---|
1799 | |
---|
1800 | |
---|
1801 | os.remove(sts_file+'.sts') |
---|
1802 | os.remove(meshname) |
---|
1803 | |
---|
1804 | |
---|
1805 | |
---|
1806 | def test_file_boundary_stsIII_ordering(self): |
---|
1807 | """test_file_boundary_stsIII_ordering(self): |
---|
1808 | Read correct points from ordering file and apply sts to boundary |
---|
1809 | """ |
---|
1810 | |
---|
1811 | lat_long_points=[[6.01,97.0],[6.02,97.0],[6.05,96.9],[6.0,97.0]] |
---|
1812 | bounding_polygon=[[6.0,97.0],[6.01,97.0],[6.02,97.0], |
---|
1813 | [6.02,97.02],[6.00,97.02]] |
---|
1814 | tide = 3.0 |
---|
1815 | time_step_count = 50 |
---|
1816 | time_step = 2 |
---|
1817 | n=len(lat_long_points) |
---|
1818 | first_tstep=num.ones(n,num.int) |
---|
1819 | last_tstep=(time_step_count)*num.ones(n,num.int) |
---|
1820 | gauge_depth=20*num.ones(n,num.float) |
---|
1821 | ha=2*num.ones((n,time_step_count),num.float) |
---|
1822 | ua=10*num.ones((n,time_step_count),num.float) |
---|
1823 | va=-10*num.ones((n,time_step_count),num.float) |
---|
1824 | base_name, files = self.write_mux2(lat_long_points, |
---|
1825 | time_step_count, |
---|
1826 | time_step, |
---|
1827 | first_tstep, |
---|
1828 | last_tstep, |
---|
1829 | depth=gauge_depth, |
---|
1830 | ha=ha, |
---|
1831 | ua=ua, |
---|
1832 | va=va) |
---|
1833 | # base name will not exist, but 3 other files are created |
---|
1834 | |
---|
1835 | # Write order file |
---|
1836 | file_handle, order_base_name = tempfile.mkstemp("") |
---|
1837 | os.close(file_handle) |
---|
1838 | os.remove(order_base_name) |
---|
1839 | d="," |
---|
1840 | order_file=order_base_name+'order.txt' |
---|
1841 | fid=open(order_file,'w') |
---|
1842 | |
---|
1843 | # Write Header |
---|
1844 | header='index, longitude, latitude\n' |
---|
1845 | fid.write(header) |
---|
1846 | indices=[3,0,1] |
---|
1847 | for i in indices: |
---|
1848 | line=str(i)+d+str(lat_long_points[i][1])+d+\ |
---|
1849 | str(lat_long_points[i][0])+"\n" |
---|
1850 | fid.write(line) |
---|
1851 | fid.close() |
---|
1852 | |
---|
1853 | sts_file=base_name |
---|
1854 | urs2sts(base_name, |
---|
1855 | basename_out=sts_file, |
---|
1856 | ordering_filename=order_file, |
---|
1857 | mean_stage=tide, |
---|
1858 | verbose=False) |
---|
1859 | self.delete_mux(files) |
---|
1860 | |
---|
1861 | assert(os.access(sts_file+'.sts', os.F_OK)) |
---|
1862 | |
---|
1863 | boundary_polygon = create_sts_boundary(base_name) |
---|
1864 | |
---|
1865 | os.remove(order_file) |
---|
1866 | |
---|
1867 | # Append the remaining part of the boundary polygon to be defined by |
---|
1868 | # the user |
---|
1869 | bounding_polygon_utm=[] |
---|
1870 | for point in bounding_polygon: |
---|
1871 | zone,easting,northing=redfearn(point[0],point[1]) |
---|
1872 | bounding_polygon_utm.append([easting,northing]) |
---|
1873 | |
---|
1874 | boundary_polygon.append(bounding_polygon_utm[3]) |
---|
1875 | boundary_polygon.append(bounding_polygon_utm[4]) |
---|
1876 | |
---|
1877 | plot=False |
---|
1878 | if plot: |
---|
1879 | from pylab import plot,show,axis |
---|
1880 | boundary_polygon=ensure_numeric(boundary_polygon) |
---|
1881 | bounding_polygon_utm=ensure_numeric(bounding_polygon_utm) |
---|
1882 | #lat_long_points=ensure_numeric(lat_long_points) |
---|
1883 | #plot(lat_long_points[:,0],lat_long_points[:,1],'o') |
---|
1884 | plot(boundary_polygon[:,0], boundary_polygon[:,1],'d') |
---|
1885 | plot(bounding_polygon_utm[:,0],bounding_polygon_utm[:,1],'o') |
---|
1886 | show() |
---|
1887 | |
---|
1888 | assert num.allclose(bounding_polygon_utm,boundary_polygon) |
---|
1889 | |
---|
1890 | |
---|
1891 | extent_res=1000000 |
---|
1892 | meshname = 'urs_test_mesh' + '.tsh' |
---|
1893 | interior_regions=None |
---|
1894 | boundary_tags={'ocean': [0,1], 'otherocean': [2,3,4]} |
---|
1895 | |
---|
1896 | # have to change boundary tags from last example because now bounding |
---|
1897 | # polygon starts in different place. |
---|
1898 | create_mesh_from_regions(boundary_polygon,boundary_tags=boundary_tags, |
---|
1899 | maximum_triangle_area=extent_res,filename=meshname, |
---|
1900 | interior_regions=interior_regions,verbose=False) |
---|
1901 | |
---|
1902 | domain_fbound = Domain(meshname) |
---|
1903 | domain_fbound.set_quantity('stage', tide) |
---|
1904 | Bf = File_boundary(sts_file+'.sts', |
---|
1905 | domain_fbound, |
---|
1906 | boundary_polygon=boundary_polygon) |
---|
1907 | Br = Reflective_boundary(domain_fbound) |
---|
1908 | |
---|
1909 | domain_fbound.set_boundary({'ocean': Bf,'otherocean': Br}) |
---|
1910 | finaltime=time_step*(time_step_count-1) |
---|
1911 | yieldstep=time_step |
---|
1912 | temp_fbound=num.zeros(int(finaltime/yieldstep)+1,num.float) |
---|
1913 | |
---|
1914 | for i, t in enumerate(domain_fbound.evolve(yieldstep=yieldstep, |
---|
1915 | finaltime=finaltime, |
---|
1916 | skip_initial_step = False)): |
---|
1917 | temp_fbound[i]=domain_fbound.quantities['stage'].centroid_values[2] |
---|
1918 | |
---|
1919 | |
---|
1920 | domain_drchlt = Domain(meshname) |
---|
1921 | domain_drchlt.set_quantity('stage', tide) |
---|
1922 | Br = Reflective_boundary(domain_drchlt) |
---|
1923 | Bd = Dirichlet_boundary([2.0+tide,220+10*tide,-220-10*tide]) |
---|
1924 | domain_drchlt.set_boundary({'ocean': Bd,'otherocean': Br}) |
---|
1925 | temp_drchlt=num.zeros(int(finaltime/yieldstep)+1,num.float) |
---|
1926 | |
---|
1927 | for i, t in enumerate(domain_drchlt.evolve(yieldstep=yieldstep, |
---|
1928 | finaltime=finaltime, |
---|
1929 | skip_initial_step = False)): |
---|
1930 | temp_drchlt[i]=domain_drchlt.quantities['stage'].centroid_values[2] |
---|
1931 | |
---|
1932 | |
---|
1933 | #print domain_fbound.quantities['stage'].vertex_values |
---|
1934 | #print domain_drchlt.quantities['stage'].vertex_values |
---|
1935 | |
---|
1936 | assert num.allclose(temp_fbound,temp_drchlt) |
---|
1937 | |
---|
1938 | |
---|
1939 | assert num.allclose(domain_fbound.quantities['stage'].vertex_values, |
---|
1940 | domain_drchlt.quantities['stage'].vertex_values) |
---|
1941 | |
---|
1942 | assert num.allclose(domain_fbound.quantities['xmomentum'].vertex_values, |
---|
1943 | domain_drchlt.quantities['xmomentum'].vertex_values) |
---|
1944 | |
---|
1945 | assert num.allclose(domain_fbound.quantities['ymomentum'].vertex_values, |
---|
1946 | domain_drchlt.quantities['ymomentum'].vertex_values) |
---|
1947 | |
---|
1948 | # Use known Dirichlet condition (if sufficient timesteps have been taken) |
---|
1949 | |
---|
1950 | #FIXME: What do these assertions test? Also do they assume tide =0 |
---|
1951 | #print domain_fbound.quantities['stage'].vertex_values |
---|
1952 | #assert allclose(domain_drchlt.quantities['stage'].vertex_values[6], 2) |
---|
1953 | #assert allclose(domain_fbound.quantities['stage'].vertex_values[6], 2) |
---|
1954 | |
---|
1955 | if not sys.platform == 'win32': |
---|
1956 | os.remove(sts_file+'.sts') |
---|
1957 | |
---|
1958 | os.remove(meshname) |
---|
1959 | |
---|
1960 | |
---|
1961 | |
---|
1962 | #------------------------------------------------------------- |
---|
1963 | |
---|
1964 | if __name__ == "__main__": |
---|
1965 | suite = unittest.makeSuite(Test_Urs2Sts,'test') |
---|
1966 | runner = unittest.TextTestRunner() |
---|
1967 | runner.run(suite) |
---|