1 | #!/usr/bin/env python |
---|
2 | |
---|
3 | import unittest, os |
---|
4 | import os.path |
---|
5 | from math import pi, sqrt |
---|
6 | import tempfile |
---|
7 | |
---|
8 | from anuga.config import g, epsilon |
---|
9 | from anuga.config import netcdf_mode_r, netcdf_mode_w, netcdf_mode_a |
---|
10 | from anuga.utilities.numerical_tools import mean |
---|
11 | from anuga.utilities.polygon import is_inside_polygon |
---|
12 | from anuga.coordinate_transforms.geo_reference import Geo_reference |
---|
13 | from anuga.abstract_2d_finite_volumes.quantity import Quantity |
---|
14 | from anuga.geospatial_data.geospatial_data import Geospatial_data |
---|
15 | from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular_cross |
---|
16 | |
---|
17 | from anuga.utilities.system_tools import get_pathname_from_package |
---|
18 | from swb_domain import * |
---|
19 | |
---|
20 | import numpy as num |
---|
21 | |
---|
22 | # Get gateway to C implementation of flux function for direct testing |
---|
23 | from shallow_water_ext import flux_function_central as flux_function |
---|
24 | |
---|
25 | |
---|
26 | |
---|
27 | |
---|
28 | class Test_swb_boundary_condition(unittest.TestCase): |
---|
29 | def setUp(self): |
---|
30 | pass |
---|
31 | |
---|
32 | def tearDown(self): |
---|
33 | pass |
---|
34 | def test_boundary_conditions(self): |
---|
35 | a = [0.0, 0.0] |
---|
36 | b = [0.0, 2.0] |
---|
37 | c = [2.0, 0.0] |
---|
38 | d = [0.0, 4.0] |
---|
39 | e = [2.0, 2.0] |
---|
40 | f = [4.0, 0.0] |
---|
41 | |
---|
42 | points = [a, b, c, d, e, f] |
---|
43 | # bac, bce, ecf, dbe |
---|
44 | vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]] |
---|
45 | boundary = {(0, 0): 'Third', |
---|
46 | (0, 2): 'First', |
---|
47 | (2, 0): 'Second', |
---|
48 | (2, 1): 'Second', |
---|
49 | (3, 1): 'Second', |
---|
50 | (3, 2): 'Third'} |
---|
51 | |
---|
52 | domain = Domain(points, vertices, boundary) |
---|
53 | domain.check_integrity() |
---|
54 | |
---|
55 | domain.set_quantity('stage', [[1,2,3], [5,5,5], [0,0,9], [-6,3,3]]) |
---|
56 | |
---|
57 | domain.set_quantity('xmomentum', [[1,1,1], [2,2,2], [3,3,3], [4,4,4]]) |
---|
58 | |
---|
59 | domain.set_quantity('ymomentum', [[10,10,10], [20,20,20], |
---|
60 | [30,30,30], [40,40,40]]) |
---|
61 | |
---|
62 | D = Dirichlet_boundary([5, 2, 1]) |
---|
63 | T = Transmissive_boundary(domain) |
---|
64 | R = Reflective_boundary(domain) |
---|
65 | domain.set_boundary({'First': D, 'Second': T, 'Third': R}) |
---|
66 | |
---|
67 | domain.set_centroid_transmissive_bc(False) |
---|
68 | domain.update_boundary() |
---|
69 | |
---|
70 | # Stage |
---|
71 | assert domain.quantities['stage'].boundary_values[0] == 2.5 |
---|
72 | # Reflective (2.5) |
---|
73 | assert (domain.quantities['stage'].boundary_values[0] == |
---|
74 | domain.get_conserved_quantities(0, edge=0)[0]) |
---|
75 | # Dirichlet |
---|
76 | assert domain.quantities['stage'].boundary_values[1] == 5. |
---|
77 | |
---|
78 | |
---|
79 | # Transmissive (4.5) |
---|
80 | assert (domain.quantities['stage'].boundary_values[2] == |
---|
81 | domain.get_conserved_quantities(2, edge=0)[0]) |
---|
82 | # Transmissive (4.5) |
---|
83 | assert (domain.quantities['stage'].boundary_values[3] == |
---|
84 | domain.get_conserved_quantities(2, edge=1)[0]) |
---|
85 | # Transmissive (-1.5) |
---|
86 | assert (domain.quantities['stage'].boundary_values[4] == |
---|
87 | domain.get_conserved_quantities(3, edge=1)[0]) |
---|
88 | # Reflective (-1.5) |
---|
89 | assert (domain.quantities['stage'].boundary_values[5] == |
---|
90 | domain.get_conserved_quantities(3, edge=2)[0]) |
---|
91 | |
---|
92 | # Xmomentum |
---|
93 | # Reflective |
---|
94 | assert domain.quantities['xmomentum'].boundary_values[0] == 1.0 |
---|
95 | # Dirichlet |
---|
96 | assert domain.quantities['xmomentum'].boundary_values[1] == 2. |
---|
97 | # Transmissive |
---|
98 | assert (domain.quantities['xmomentum'].boundary_values[2] == |
---|
99 | domain.get_conserved_quantities(2, edge=0)[1]) |
---|
100 | # Transmissive |
---|
101 | assert (domain.quantities['xmomentum'].boundary_values[3] == |
---|
102 | domain.get_conserved_quantities(2, edge=1)[1]) |
---|
103 | # Transmissive |
---|
104 | assert (domain.quantities['xmomentum'].boundary_values[4] == |
---|
105 | domain.get_conserved_quantities(3, edge=1)[1]) |
---|
106 | # Reflective |
---|
107 | assert domain.quantities['xmomentum'].boundary_values[5] == -4.0 |
---|
108 | |
---|
109 | # Ymomentum |
---|
110 | # Reflective |
---|
111 | assert domain.quantities['ymomentum'].boundary_values[0] == -10.0 |
---|
112 | # Dirichlet |
---|
113 | assert domain.quantities['ymomentum'].boundary_values[1] == 1. |
---|
114 | # Transmissive |
---|
115 | assert domain.quantities['ymomentum'].boundary_values[2] == 30. |
---|
116 | # Transmissive |
---|
117 | assert domain.quantities['ymomentum'].boundary_values[3] == 30. |
---|
118 | # Transmissive |
---|
119 | assert domain.quantities['ymomentum'].boundary_values[4] == 40. |
---|
120 | # Reflective |
---|
121 | assert domain.quantities['ymomentum'].boundary_values[5] == 40. |
---|
122 | |
---|
123 | def test_boundary_conditionsII(self): |
---|
124 | a = [0.0, 0.0] |
---|
125 | b = [0.0, 2.0] |
---|
126 | c = [2.0, 0.0] |
---|
127 | d = [0.0, 4.0] |
---|
128 | e = [2.0, 2.0] |
---|
129 | f = [4.0, 0.0] |
---|
130 | |
---|
131 | points = [a, b, c, d, e, f] |
---|
132 | # bac, bce, ecf, dbe |
---|
133 | vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]] |
---|
134 | boundary = {(0, 0): 'Third', |
---|
135 | (0, 2): 'First', |
---|
136 | (2, 0): 'Second', |
---|
137 | (2, 1): 'Second', |
---|
138 | (3, 1): 'Second', |
---|
139 | (3, 2): 'Third', |
---|
140 | (0, 1): 'Internal'} |
---|
141 | |
---|
142 | domain = Domain(points, vertices, boundary) |
---|
143 | domain.check_integrity() |
---|
144 | |
---|
145 | domain.set_quantity('stage', [[1,2,3], [5,5,5], [0,0,9], [-6,3,3]]) |
---|
146 | |
---|
147 | domain.set_quantity('xmomentum', [[1,1,1], [2,2,2], [3,3,3], [4,4,4]]) |
---|
148 | |
---|
149 | domain.set_quantity('ymomentum', [[10,10,10], [20,20,20], |
---|
150 | [30,30,30], [40,40,40]]) |
---|
151 | |
---|
152 | D = Dirichlet_boundary([5, 2, 1]) |
---|
153 | T = Transmissive_boundary(domain) |
---|
154 | R = Reflective_boundary(domain) |
---|
155 | domain.set_boundary({'First': D, 'Second': T, |
---|
156 | 'Third': R, 'Internal': None}) |
---|
157 | |
---|
158 | domain.update_boundary() |
---|
159 | domain.check_integrity() |
---|
160 | |
---|
161 | def test_boundary_conditionsIII(self): |
---|
162 | """test_boundary_conditionsIII |
---|
163 | |
---|
164 | Test Transmissive_stage_zero_momentum_boundary |
---|
165 | """ |
---|
166 | |
---|
167 | a = [0.0, 0.0] |
---|
168 | b = [0.0, 2.0] |
---|
169 | c = [2.0, 0.0] |
---|
170 | d = [0.0, 4.0] |
---|
171 | e = [2.0, 2.0] |
---|
172 | f = [4.0, 0.0] |
---|
173 | |
---|
174 | points = [a, b, c, d, e, f] |
---|
175 | # bac, bce, ecf, dbe |
---|
176 | vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]] |
---|
177 | boundary = {(0, 0): 'Third', |
---|
178 | (0, 2): 'First', |
---|
179 | (2, 0): 'Second', |
---|
180 | (2, 1): 'Second', |
---|
181 | (3, 1): 'Second', |
---|
182 | (3, 2): 'Third'} |
---|
183 | |
---|
184 | domain = Domain(points, vertices, boundary) |
---|
185 | domain.check_integrity() |
---|
186 | |
---|
187 | domain.set_quantity('stage', [[1,2,3], [5,5,5], [0,0,9], [-6,3,3]]) |
---|
188 | |
---|
189 | domain.set_quantity('xmomentum', [[1,1,1], [2,2,2], [3,3,3], [4,4,4]]) |
---|
190 | |
---|
191 | domain.set_quantity('ymomentum', [[10,10,10], [20,20,20], |
---|
192 | [30,30,30], [40,40,40]]) |
---|
193 | |
---|
194 | D = Dirichlet_boundary([5, 2, 1]) |
---|
195 | T = Transmissive_stage_zero_momentum_boundary(domain) |
---|
196 | R = Reflective_boundary(domain) |
---|
197 | domain.set_boundary({'First': D, 'Second': T, 'Third': R}) |
---|
198 | |
---|
199 | domain.update_boundary() |
---|
200 | |
---|
201 | # Stage |
---|
202 | # Reflective (2.5) |
---|
203 | assert domain.quantities['stage'].boundary_values[0] == 2.5 |
---|
204 | assert (domain.quantities['stage'].boundary_values[0] == |
---|
205 | domain.get_conserved_quantities(0, edge=0)[0]) |
---|
206 | # Dirichlet |
---|
207 | assert domain.quantities['stage'].boundary_values[1] == 5. |
---|
208 | # Transmissive (4.5) |
---|
209 | assert (domain.quantities['stage'].boundary_values[2] == |
---|
210 | domain.get_conserved_quantities(2, edge=0)[0]) |
---|
211 | # Transmissive (4.5) |
---|
212 | assert (domain.quantities['stage'].boundary_values[3] == |
---|
213 | domain.get_conserved_quantities(2, edge=1)[0]) |
---|
214 | # Transmissive (-1.5) |
---|
215 | assert (domain.quantities['stage'].boundary_values[4] == |
---|
216 | domain.get_conserved_quantities(3, edge=1)[0]) |
---|
217 | # Reflective (-1.5) |
---|
218 | assert (domain.quantities['stage'].boundary_values[5] == |
---|
219 | domain.get_conserved_quantities(3, edge=2)[0]) |
---|
220 | |
---|
221 | # Xmomentum |
---|
222 | # Reflective |
---|
223 | assert domain.quantities['xmomentum'].boundary_values[0] == 1.0 |
---|
224 | # Dirichlet |
---|
225 | assert domain.quantities['xmomentum'].boundary_values[1] == 2. |
---|
226 | assert domain.quantities['xmomentum'].boundary_values[2] == 0.0 |
---|
227 | assert domain.quantities['xmomentum'].boundary_values[3] == 0.0 |
---|
228 | assert domain.quantities['xmomentum'].boundary_values[4] == 0.0 |
---|
229 | # Reflective |
---|
230 | assert domain.quantities['xmomentum'].boundary_values[5] == -4.0 |
---|
231 | |
---|
232 | # Ymomentum |
---|
233 | # Reflective |
---|
234 | assert domain.quantities['ymomentum'].boundary_values[0] == -10.0 |
---|
235 | # Dirichlet |
---|
236 | assert domain.quantities['ymomentum'].boundary_values[1] == 1. |
---|
237 | assert domain.quantities['ymomentum'].boundary_values[2] == 0.0 |
---|
238 | assert domain.quantities['ymomentum'].boundary_values[3] == 0.0 |
---|
239 | assert domain.quantities['ymomentum'].boundary_values[4] == 0.0 |
---|
240 | # Reflective |
---|
241 | assert domain.quantities['ymomentum'].boundary_values[5] == 40. |
---|
242 | |
---|
243 | def test_boundary_condition_time(self): |
---|
244 | """test_boundary_condition_time |
---|
245 | |
---|
246 | This tests that boundary conditions are evaluated |
---|
247 | using the right time from domain. |
---|
248 | """ |
---|
249 | |
---|
250 | # Setup computational domain |
---|
251 | from anuga.abstract_2d_finite_volumes.mesh_factory \ |
---|
252 | import rectangular_cross |
---|
253 | |
---|
254 | #-------------------------------------------------------------- |
---|
255 | # Setup computational domain |
---|
256 | #-------------------------------------------------------------- |
---|
257 | N = 5 |
---|
258 | points, vertices, boundary = rectangular_cross(N, N) |
---|
259 | domain = Domain(points, vertices, boundary) |
---|
260 | |
---|
261 | #-------------------------------------------------------------- |
---|
262 | # Setup initial conditions |
---|
263 | #-------------------------------------------------------------- |
---|
264 | domain.set_quantity('elevation', 0.0) |
---|
265 | domain.set_quantity('friction', 0.0) |
---|
266 | domain.set_quantity('stage', 0.0) |
---|
267 | |
---|
268 | #-------------------------------------------------------------- |
---|
269 | # Setup boundary conditions |
---|
270 | #-------------------------------------------------------------- |
---|
271 | # Time dependent boundary |
---|
272 | Bt = Time_boundary(domain=domain, f=lambda t: [t, 0.0, 0.0]) |
---|
273 | |
---|
274 | Br = Reflective_boundary(domain) # Reflective wall |
---|
275 | |
---|
276 | domain.set_boundary({'left': Bt, 'right': Br, 'top': Br, 'bottom': Br}) |
---|
277 | |
---|
278 | for t in domain.evolve(yieldstep = 10, finaltime = 20.0): |
---|
279 | q = Bt.evaluate() |
---|
280 | |
---|
281 | # FIXME (Ole): This test would not have passed in |
---|
282 | # changeset:5846. |
---|
283 | msg = 'Time boundary not evaluated correctly' |
---|
284 | assert num.allclose(t, q[0]), msg |
---|
285 | |
---|
286 | |
---|
287 | def test_spatio_temporal_boundary_1(self): |
---|
288 | """Test that boundary values can be read from file and interpolated |
---|
289 | in both time and space. |
---|
290 | |
---|
291 | Verify that the same steady state solution is arrived at and that |
---|
292 | time interpolation works. |
---|
293 | |
---|
294 | The full solution history is not exactly the same as |
---|
295 | file boundary must read and interpolate from *smoothed* version |
---|
296 | as stored in sww. |
---|
297 | """ |
---|
298 | |
---|
299 | # Create sww file of simple propagation from left to right |
---|
300 | # through rectangular domain |
---|
301 | |
---|
302 | import time |
---|
303 | from mesh_factory import rectangular |
---|
304 | |
---|
305 | # Create basic mesh |
---|
306 | points, vertices, boundary = rectangular(3, 3) |
---|
307 | |
---|
308 | # Create shallow water domain |
---|
309 | domain1 = Domain(points, vertices, boundary) |
---|
310 | |
---|
311 | domain1.reduction = mean |
---|
312 | domain1.smooth = False # Exact result |
---|
313 | |
---|
314 | domain1.default_order = 2 |
---|
315 | domain1.store = True |
---|
316 | domain1.set_datadir('.') |
---|
317 | domain1.set_name('spatio_temporal_boundary_source' + str(time.time())) |
---|
318 | |
---|
319 | # FIXME: This is extremely important! |
---|
320 | # How can we test if they weren't stored? |
---|
321 | |
---|
322 | |
---|
323 | # Bed-slope and friction at vertices (and interpolated elsewhere) |
---|
324 | domain1.set_quantity('elevation', 0) |
---|
325 | domain1.set_quantity('friction', 0) |
---|
326 | |
---|
327 | # Boundary conditions |
---|
328 | Br = Reflective_boundary(domain1) |
---|
329 | Bd = Dirichlet_boundary([0.3,0,0]) |
---|
330 | domain1.set_boundary({'left': Bd, 'top': Bd, 'right': Br, 'bottom': Br}) |
---|
331 | |
---|
332 | # Initial condition |
---|
333 | domain1.set_quantity('stage', 0) |
---|
334 | domain1.check_integrity() |
---|
335 | |
---|
336 | finaltime = 5 |
---|
337 | |
---|
338 | # Evolution (full domain - large steps) |
---|
339 | for t in domain1.evolve(yieldstep=0.671, finaltime=finaltime): |
---|
340 | pass |
---|
341 | |
---|
342 | cv1 = domain1.quantities['stage'].centroid_values |
---|
343 | |
---|
344 | # Create a triangle shaped domain (reusing coordinates from domain 1), |
---|
345 | # formed from the lower and right hand boundaries and |
---|
346 | # the sw-ne diagonal |
---|
347 | # from domain 1. Call it domain2 |
---|
348 | |
---|
349 | points = [[0,0], [1.0/3,0], [1.0/3,1.0/3], |
---|
350 | [2.0/3,0], [2.0/3,1.0/3], [2.0/3,2.0/3], |
---|
351 | [1,0], [1,1.0/3], [1,2.0/3], [1,1]] |
---|
352 | |
---|
353 | vertices = [[1,2,0], [3,4,1], [2,1,4], [4,5,2], |
---|
354 | [6,7,3], [4,3,7], [7,8,4], [5,4,8], [8,9,5]] |
---|
355 | |
---|
356 | boundary = {(0,1): 'bottom', (1,1): 'bottom', (4,1): 'bottom', |
---|
357 | (4,2): 'right', (6,2): 'right', (8,2): 'right', |
---|
358 | (0,0): 'diagonal', (3,0): 'diagonal', (8,0): 'diagonal'} |
---|
359 | |
---|
360 | domain2 = Domain(points, vertices, boundary) |
---|
361 | |
---|
362 | domain2.reduction = domain1.reduction |
---|
363 | domain2.smooth = False |
---|
364 | domain2.default_order = 2 |
---|
365 | |
---|
366 | # Bed-slope and friction at vertices (and interpolated elsewhere) |
---|
367 | domain2.set_quantity('elevation', 0) |
---|
368 | domain2.set_quantity('friction', 0) |
---|
369 | domain2.set_quantity('stage', 0) |
---|
370 | |
---|
371 | # Boundary conditions |
---|
372 | Br = Reflective_boundary(domain2) |
---|
373 | Bf = Field_boundary(domain1.get_name() + '.sww', domain2) |
---|
374 | domain2.set_boundary({'right':Br, 'bottom':Br, 'diagonal':Bf}) |
---|
375 | domain2.check_integrity() |
---|
376 | |
---|
377 | # Evolution (small steps) |
---|
378 | for t in domain2.evolve(yieldstep=0.0711, finaltime=finaltime): |
---|
379 | pass |
---|
380 | |
---|
381 | # Use output from domain1 as spatio-temporal boundary for domain2 |
---|
382 | # and verify that results at right hand side are close. |
---|
383 | cv2 = domain2.quantities['stage'].centroid_values |
---|
384 | |
---|
385 | |
---|
386 | assert num.allclose(num.take(cv1, (0,8,16), axis=0), |
---|
387 | num.take(cv2, (0,3,8), axis=0),atol=1.0e-2) # Diag |
---|
388 | assert num.allclose(num.take(cv1, (0,6,12), axis=0), |
---|
389 | num.take(cv2, (0,1,4), axis=0),atol=1.0e-2) # Bottom |
---|
390 | assert num.allclose(num.take(cv1, (12,14,16), axis=0), |
---|
391 | num.take(cv2, (4,6,8), axis=0),atol=1.0e-2) # RHS |
---|
392 | |
---|
393 | # Cleanup |
---|
394 | os.remove(domain1.get_name() + '.sww') |
---|
395 | |
---|
396 | def test_spatio_temporal_boundary_2(self): |
---|
397 | """Test that boundary values can be read from file and interpolated |
---|
398 | in both time and space. |
---|
399 | This is a more basic test, verifying that boundary object |
---|
400 | produces the expected results |
---|
401 | """ |
---|
402 | |
---|
403 | import time |
---|
404 | from mesh_factory import rectangular |
---|
405 | |
---|
406 | # Create sww file of simple propagation from left to right |
---|
407 | # through rectangular domain |
---|
408 | |
---|
409 | # Create basic mesh |
---|
410 | points, vertices, boundary = rectangular(3, 3) |
---|
411 | |
---|
412 | #Create shallow water domain |
---|
413 | domain1 = Domain(points, vertices, boundary) |
---|
414 | |
---|
415 | domain1.reduction = mean |
---|
416 | domain1.smooth = True # To mimic MOST output |
---|
417 | |
---|
418 | domain1.default_order = 2 |
---|
419 | domain1.store = True |
---|
420 | domain1.set_datadir('.') |
---|
421 | domain1.set_name('spatio_temporal_boundary_source' + str(time.time())) |
---|
422 | |
---|
423 | # FIXME: This is extremely important! |
---|
424 | # How can we test if they weren't stored? |
---|
425 | domain1.quantities_to_be_stored = {'elevation': 1, |
---|
426 | 'stage': 2, |
---|
427 | 'xmomentum': 2, |
---|
428 | 'ymomentum': 2} |
---|
429 | |
---|
430 | # Bed-slope and friction at vertices (and interpolated elsewhere) |
---|
431 | domain1.set_quantity('elevation', 0) |
---|
432 | domain1.set_quantity('friction', 0) |
---|
433 | |
---|
434 | # Boundary conditions |
---|
435 | Br = Reflective_boundary(domain1) |
---|
436 | Bd = Dirichlet_boundary([0.3,0,0]) |
---|
437 | domain1.set_boundary({'left': Bd, 'top': Bd, 'right': Br, 'bottom': Br}) |
---|
438 | |
---|
439 | # Initial condition |
---|
440 | domain1.set_quantity('stage', 0) |
---|
441 | domain1.check_integrity() |
---|
442 | |
---|
443 | finaltime = 5 |
---|
444 | |
---|
445 | # Evolution (full domain - large steps) |
---|
446 | for t in domain1.evolve(yieldstep=1, finaltime=finaltime): |
---|
447 | pass |
---|
448 | |
---|
449 | # Create an triangle shaped domain (coinciding with some |
---|
450 | # coordinates from domain 1), |
---|
451 | # formed from the lower and right hand boundaries and |
---|
452 | # the sw-ne diagonal |
---|
453 | # from domain 1. Call it domain2 |
---|
454 | points = [[0,0], |
---|
455 | [1.0/3,0], [1.0/3,1.0/3], |
---|
456 | [2.0/3,0], [2.0/3,1.0/3], [2.0/3,2.0/3], |
---|
457 | [1,0], [1,1.0/3], [1,2.0/3], [1,1]] |
---|
458 | |
---|
459 | vertices = [[1,2,0], [3,4,1], [2,1,4], [4,5,2], |
---|
460 | [6,7,3], [4,3,7], [7,8,4], [5,4,8], [8,9,5]] |
---|
461 | |
---|
462 | boundary = {(0,1): 'bottom', (1,1): 'bottom', (4,1): 'bottom', |
---|
463 | (4,2): 'right', (6,2): 'right', (8,2): 'right', |
---|
464 | (0,0): 'diagonal', (3,0): 'diagonal', (8,0): 'diagonal'} |
---|
465 | |
---|
466 | domain2 = Domain(points, vertices, boundary) |
---|
467 | |
---|
468 | domain2.reduction = domain1.reduction |
---|
469 | domain2.smooth = False |
---|
470 | domain2.default_order = 2 |
---|
471 | |
---|
472 | # Bed-slope and friction at vertices (and interpolated elsewhere) |
---|
473 | domain2.set_quantity('elevation', 0) |
---|
474 | domain2.set_quantity('friction', 0) |
---|
475 | domain2.set_quantity('stage', 0) |
---|
476 | |
---|
477 | # Read results for specific timesteps t=1 and t=2 |
---|
478 | from Scientific.IO.NetCDF import NetCDFFile |
---|
479 | fid = NetCDFFile(domain1.get_name() + '.sww') |
---|
480 | |
---|
481 | x = fid.variables['x'][:] |
---|
482 | y = fid.variables['y'][:] |
---|
483 | s1 = fid.variables['stage'][1,:] |
---|
484 | s2 = fid.variables['stage'][2,:] |
---|
485 | fid.close() |
---|
486 | |
---|
487 | shp = (len(x), 1) |
---|
488 | points = num.concatenate((num.reshape(x, shp), num.reshape(y, shp)), |
---|
489 | axis=1) |
---|
490 | |
---|
491 | # The diagonal points of domain 1 are 0, 5, 10, 15 |
---|
492 | msg = ('value was\n%s\nshould be\n' |
---|
493 | '[[0,0], [1.0/3, 1.0/3],\n' |
---|
494 | '[2.0/3, 2.0/3], [1,1]]' |
---|
495 | % str(num.take(points, [0,5,10,15], axis=0))) |
---|
496 | assert num.allclose(num.take(points, [0,5,10,15], axis=0), |
---|
497 | [[0,0], [1.0/3, 1.0/3], [2.0/3, 2.0/3], [1,1]]), msg |
---|
498 | |
---|
499 | # Boundary conditions |
---|
500 | Br = Reflective_boundary(domain2) |
---|
501 | Bf = Field_boundary(domain1.get_name() + '.sww', |
---|
502 | domain2, verbose=False) |
---|
503 | domain2.set_boundary({'right':Br, 'bottom':Br, 'diagonal':Bf}) |
---|
504 | domain2.check_integrity() |
---|
505 | |
---|
506 | # Test that interpolation points are the mid points of the all boundary |
---|
507 | # segments |
---|
508 | boundary_midpoints = [[1.0/6, 0], [1.0/2, 0], [5.0/6,0], |
---|
509 | [1.0, 1.0/6], [1.0, 1.0/2], [1.0, 5.0/6], |
---|
510 | [1.0/6, 1.0/6], [0.5, 0.5], [5.0/6, 5.0/6]] |
---|
511 | |
---|
512 | boundary_midpoints.sort() |
---|
513 | R = Bf.F.interpolation_points.tolist() |
---|
514 | R.sort() |
---|
515 | assert num.allclose(boundary_midpoints, R) |
---|
516 | |
---|
517 | # Check spatially interpolated output at time == 1 |
---|
518 | domain2.time = 1 |
---|
519 | |
---|
520 | # First diagonal midpoint |
---|
521 | R0 = Bf.evaluate(0, 0) |
---|
522 | assert num.allclose(R0[0], (s1[0] + s1[5])/2) |
---|
523 | |
---|
524 | # Second diagonal midpoint |
---|
525 | R0 = Bf.evaluate(3, 0) |
---|
526 | assert num.allclose(R0[0], (s1[5] + s1[10])/2) |
---|
527 | |
---|
528 | # First diagonal midpoint |
---|
529 | R0 = Bf.evaluate(8, 0) |
---|
530 | assert num.allclose(R0[0], (s1[10] + s1[15])/2) |
---|
531 | |
---|
532 | # Check spatially interpolated output at time == 2 |
---|
533 | domain2.time = 2 |
---|
534 | |
---|
535 | # First diagonal midpoint |
---|
536 | R0 = Bf.evaluate(0, 0) |
---|
537 | assert num.allclose(R0[0], (s2[0] + s2[5])/2) |
---|
538 | |
---|
539 | # Second diagonal midpoint |
---|
540 | R0 = Bf.evaluate(3, 0) |
---|
541 | assert num.allclose(R0[0], (s2[5] + s2[10])/2) |
---|
542 | |
---|
543 | # First diagonal midpoint |
---|
544 | R0 = Bf.evaluate(8, 0) |
---|
545 | assert num.allclose(R0[0], (s2[10] + s2[15])/2) |
---|
546 | |
---|
547 | # Now check temporal interpolation |
---|
548 | domain2.time = 1 + 2.0/3 |
---|
549 | |
---|
550 | # First diagonal midpoint |
---|
551 | R0 = Bf.evaluate(0,0) |
---|
552 | assert num.allclose(R0[0], |
---|
553 | ((s1[0] + s1[5])/2 + 2.0*(s2[0] + s2[5])/2)/3) |
---|
554 | |
---|
555 | # Second diagonal midpoint |
---|
556 | R0 = Bf.evaluate(3, 0) |
---|
557 | assert num.allclose(R0[0], |
---|
558 | ((s1[5] + s1[10])/2 + 2.0*(s2[5] + s2[10])/2)/3) |
---|
559 | |
---|
560 | # First diagonal midpoint |
---|
561 | R0 = Bf.evaluate(8, 0) |
---|
562 | assert num.allclose(R0[0], |
---|
563 | ((s1[10] + s1[15])/2 + 2.0*(s2[10] + s2[15])/2)/3) |
---|
564 | |
---|
565 | # Cleanup |
---|
566 | os.remove(domain1.get_name() + '.sww') |
---|
567 | |
---|
568 | def test_spatio_temporal_boundary_3(self): |
---|
569 | """Test that boundary values can be read from file and interpolated |
---|
570 | in both time and space. |
---|
571 | This is a more basic test, verifying that boundary object |
---|
572 | produces the expected results |
---|
573 | |
---|
574 | This tests adjusting using mean_stage |
---|
575 | """ |
---|
576 | |
---|
577 | #Create sww file of simple propagation from left to right |
---|
578 | #through rectangular domain |
---|
579 | |
---|
580 | import time |
---|
581 | from mesh_factory import rectangular |
---|
582 | |
---|
583 | mean_stage = 5.2 # Adjust stage by this amount in boundary |
---|
584 | |
---|
585 | # Create basic mesh |
---|
586 | points, vertices, boundary = rectangular(3, 3) |
---|
587 | |
---|
588 | # Create shallow water domain |
---|
589 | domain1 = Domain(points, vertices, boundary) |
---|
590 | |
---|
591 | domain1.reduction = mean |
---|
592 | domain1.smooth = True # To mimic MOST output |
---|
593 | |
---|
594 | domain1.default_order = 2 |
---|
595 | domain1.store = True |
---|
596 | domain1.set_datadir('.') |
---|
597 | domain1.set_name('spatio_temporal_boundary_source' + str(time.time())) |
---|
598 | |
---|
599 | # FIXME: This is extremely important! |
---|
600 | # How can we test if they weren't stored? |
---|
601 | |
---|
602 | # Bed-slope and friction at vertices (and interpolated elsewhere) |
---|
603 | domain1.set_quantity('elevation', 0) |
---|
604 | domain1.set_quantity('friction', 0) |
---|
605 | |
---|
606 | # Boundary conditions |
---|
607 | Br = Reflective_boundary(domain1) |
---|
608 | Bd = Dirichlet_boundary([0.3, 0, 0]) |
---|
609 | domain1.set_boundary({'left': Bd, 'top': Bd, 'right': Br, 'bottom': Br}) |
---|
610 | |
---|
611 | # Initial condition |
---|
612 | domain1.set_quantity('stage', 0) |
---|
613 | domain1.check_integrity() |
---|
614 | |
---|
615 | finaltime = 5 |
---|
616 | |
---|
617 | # Evolution (full domain - large steps) |
---|
618 | for t in domain1.evolve(yieldstep=1, finaltime=finaltime): |
---|
619 | pass |
---|
620 | |
---|
621 | # Create an triangle shaped domain (coinciding with some |
---|
622 | # coordinates from domain 1), |
---|
623 | # formed from the lower and right hand boundaries and |
---|
624 | # the sw-ne diagonal |
---|
625 | # from domain 1. Call it domain2 |
---|
626 | points = [[0,0], |
---|
627 | [1.0/3,0], [1.0/3,1.0/3], |
---|
628 | [2.0/3,0], [2.0/3,1.0/3], [2.0/3,2.0/3], |
---|
629 | [1,0], [1,1.0/3], [1,2.0/3], [1,1]] |
---|
630 | |
---|
631 | vertices = [[1,2,0], |
---|
632 | [3,4,1], [2,1,4], [4,5,2], |
---|
633 | [6,7,3], [4,3,7], [7,8,4], [5,4,8], [8,9,5]] |
---|
634 | |
---|
635 | boundary = {(0,1): 'bottom', (1,1): 'bottom', (4,1): 'bottom', |
---|
636 | (4,2): 'right', (6,2): 'right', (8,2): 'right', |
---|
637 | (0,0): 'diagonal', (3,0): 'diagonal', (8,0): 'diagonal'} |
---|
638 | |
---|
639 | domain2 = Domain(points, vertices, boundary) |
---|
640 | |
---|
641 | domain2.reduction = domain1.reduction |
---|
642 | domain2.smooth = False |
---|
643 | domain2.default_order = 2 |
---|
644 | |
---|
645 | # Bed-slope and friction at vertices (and interpolated elsewhere) |
---|
646 | domain2.set_quantity('elevation', 0) |
---|
647 | domain2.set_quantity('friction', 0) |
---|
648 | domain2.set_quantity('stage', 0) |
---|
649 | |
---|
650 | # Read results for specific timesteps t=1 and t=2 |
---|
651 | from Scientific.IO.NetCDF import NetCDFFile |
---|
652 | fid = NetCDFFile(domain1.get_name() + '.sww') |
---|
653 | |
---|
654 | x = fid.variables['x'][:] |
---|
655 | y = fid.variables['y'][:] |
---|
656 | s1 = fid.variables['stage'][1,:] |
---|
657 | s2 = fid.variables['stage'][2,:] |
---|
658 | fid.close() |
---|
659 | |
---|
660 | shp = (len(x), 1) |
---|
661 | points = num.concatenate((num.reshape(x, shp), num.reshape(y, shp)), |
---|
662 | axis=1) |
---|
663 | #The diagonal points of domain 1 are 0, 5, 10, 15 |
---|
664 | |
---|
665 | msg = ('values was\n%s\nshould be\n' |
---|
666 | '[[0,0], [1.0/3, 1.0/3],\n' |
---|
667 | '[2.0/3, 2.0/3], [1,1]]' |
---|
668 | % str(num.take(points, [0,5,10,15], axis=0))) |
---|
669 | assert num.allclose(num.take(points, [0,5,10,15], axis=0), |
---|
670 | [[0,0], [1.0/3, 1.0/3], [2.0/3, 2.0/3], [1,1]]), msg |
---|
671 | |
---|
672 | # Boundary conditions |
---|
673 | Br = Reflective_boundary(domain2) |
---|
674 | Bf = Field_boundary(domain1.get_name() + '.sww', |
---|
675 | domain2, mean_stage=mean_stage, verbose=False) |
---|
676 | |
---|
677 | domain2.set_boundary({'right':Br, 'bottom':Br, 'diagonal':Bf}) |
---|
678 | domain2.check_integrity() |
---|
679 | |
---|
680 | # Test that interpolation points are the mid points of the all boundary |
---|
681 | # segments |
---|
682 | boundary_midpoints = [[1.0/6, 0], [1.0/2, 0], [5.0/6,0], |
---|
683 | [1.0, 1.0/6], [1.0, 1.0/2], [1.0, 5.0/6], |
---|
684 | [1.0/6, 1.0/6], [0.5, 0.5], [5.0/6, 5.0/6]] |
---|
685 | |
---|
686 | boundary_midpoints.sort() |
---|
687 | R = Bf.F.interpolation_points.tolist() |
---|
688 | R.sort() |
---|
689 | assert num.allclose(boundary_midpoints, R) |
---|
690 | |
---|
691 | # Check spatially interpolated output at time == 1 |
---|
692 | domain2.time = 1 |
---|
693 | |
---|
694 | # First diagonal midpoint |
---|
695 | R0 = Bf.evaluate(0, 0) |
---|
696 | assert num.allclose(R0[0], (s1[0] + s1[5])/2 + mean_stage) |
---|
697 | |
---|
698 | # Second diagonal midpoint |
---|
699 | R0 = Bf.evaluate(3, 0) |
---|
700 | assert num.allclose(R0[0], (s1[5] + s1[10])/2 + mean_stage) |
---|
701 | |
---|
702 | # First diagonal midpoint |
---|
703 | R0 = Bf.evaluate(8, 0) |
---|
704 | assert num.allclose(R0[0], (s1[10] + s1[15])/2 + mean_stage) |
---|
705 | |
---|
706 | # Check spatially interpolated output at time == 2 |
---|
707 | domain2.time = 2 |
---|
708 | |
---|
709 | # First diagonal midpoint |
---|
710 | R0 = Bf.evaluate(0, 0) |
---|
711 | assert num.allclose(R0[0], (s2[0] + s2[5])/2 + mean_stage) |
---|
712 | |
---|
713 | # Second diagonal midpoint |
---|
714 | R0 = Bf.evaluate(3, 0) |
---|
715 | assert num.allclose(R0[0], (s2[5] + s2[10])/2 + mean_stage) |
---|
716 | |
---|
717 | # First diagonal midpoint |
---|
718 | R0 = Bf.evaluate(8, 0) |
---|
719 | assert num.allclose(R0[0], (s2[10] + s2[15])/2 + mean_stage) |
---|
720 | |
---|
721 | #Now check temporal interpolation |
---|
722 | domain2.time = 1 + 2.0/3 |
---|
723 | |
---|
724 | # First diagonal midpoint |
---|
725 | R0 = Bf.evaluate(0, 0) |
---|
726 | assert num.allclose(R0[0], |
---|
727 | ((s1[0] + s1[5])/2 + 2.0*(s2[0] + s2[5])/2)/3 + |
---|
728 | mean_stage) |
---|
729 | |
---|
730 | # Second diagonal midpoint |
---|
731 | R0 = Bf.evaluate(3, 0) |
---|
732 | assert num.allclose(R0[0], |
---|
733 | ((s1[5] + s1[10])/2 + 2.0*(s2[5] + s2[10])/2)/3 + |
---|
734 | mean_stage) |
---|
735 | |
---|
736 | # First diagonal midpoint |
---|
737 | R0 = Bf.evaluate(8, 0) |
---|
738 | assert num.allclose(R0[0], |
---|
739 | ((s1[10] + s1[15])/2 + 2.0*(s2[10] + s2[15])/2)/3 + |
---|
740 | mean_stage) |
---|
741 | |
---|
742 | # Cleanup |
---|
743 | os.remove(domain1.get_name() + '.sww') |
---|
744 | os.remove(domain2.get_name() + '.sww') |
---|
745 | |
---|
746 | |
---|
747 | |
---|
748 | def test_spatio_temporal_boundary_outside(self): |
---|
749 | """Test that field_boundary catches if a point is outside the sww |
---|
750 | that defines it |
---|
751 | """ |
---|
752 | |
---|
753 | import time |
---|
754 | from mesh_factory import rectangular |
---|
755 | |
---|
756 | # Create sww file of simple propagation from left to right |
---|
757 | # through rectangular domain |
---|
758 | |
---|
759 | # Create basic mesh |
---|
760 | points, vertices, boundary = rectangular(3, 3) |
---|
761 | |
---|
762 | # Create shallow water domain |
---|
763 | domain1 = Domain(points, vertices, boundary) |
---|
764 | |
---|
765 | domain1.reduction = mean |
---|
766 | domain1.smooth = True # To mimic MOST output |
---|
767 | |
---|
768 | domain1.default_order = 2 |
---|
769 | domain1.store = True |
---|
770 | domain1.set_datadir('.') |
---|
771 | domain1.set_name('spatio_temporal_boundary_source' + str(time.time())) |
---|
772 | |
---|
773 | # FIXME: This is extremely important! |
---|
774 | # How can we test if they weren't stored? |
---|
775 | domain1.set_quantities_to_be_stored({'elevation': 1, |
---|
776 | 'stage': 2, |
---|
777 | 'xmomentum': 2, |
---|
778 | 'ymomentum': 2}) |
---|
779 | |
---|
780 | |
---|
781 | # Bed-slope and friction at vertices (and interpolated elsewhere) |
---|
782 | domain1.set_quantity('elevation', 0) |
---|
783 | domain1.set_quantity('friction', 0) |
---|
784 | |
---|
785 | # Boundary conditions |
---|
786 | Br = Reflective_boundary(domain1) |
---|
787 | Bd = Dirichlet_boundary([0.3, 0, 0]) |
---|
788 | domain1.set_boundary({'left': Bd, 'top': Bd, 'right': Br, 'bottom': Br}) |
---|
789 | |
---|
790 | # Initial condition |
---|
791 | domain1.set_quantity('stage', 0) |
---|
792 | domain1.check_integrity() |
---|
793 | |
---|
794 | finaltime = 5 |
---|
795 | |
---|
796 | # Evolution (full domain - large steps) |
---|
797 | for t in domain1.evolve(yieldstep=1, finaltime=finaltime): |
---|
798 | pass |
---|
799 | |
---|
800 | # Create an triangle shaped domain (coinciding with some |
---|
801 | # coordinates from domain 1, but one edge outside!), |
---|
802 | # formed from the lower and right hand boundaries and |
---|
803 | # the sw-ne diagonal as in the previous test but scaled |
---|
804 | # in the x direction by a factor of 2 |
---|
805 | points = [[0,0], |
---|
806 | [2.0/3,0], [2.0/3,1.0/3], |
---|
807 | [4.0/3,0], [4.0/3,1.0/3], [4.0/3,2.0/3], |
---|
808 | [2,0], [2,1.0/3], [2,2.0/3], [2,1]] |
---|
809 | |
---|
810 | vertices = [[1,2,0], |
---|
811 | [3,4,1], [2,1,4], [4,5,2], |
---|
812 | [6,7,3], [4,3,7], [7,8,4], [5,4,8], [8,9,5]] |
---|
813 | |
---|
814 | boundary = {(0,1): 'bottom', (1,1): 'bottom', (4,1): 'bottom', |
---|
815 | (4,2): 'right', (6,2): 'right', (8,2): 'right', |
---|
816 | (0,0): 'diagonal', (3,0): 'diagonal', (8,0): 'diagonal'} |
---|
817 | |
---|
818 | domain2 = Domain(points, vertices, boundary) |
---|
819 | |
---|
820 | domain2.reduction = domain1.reduction |
---|
821 | domain2.smooth = False |
---|
822 | domain2.default_order = 2 |
---|
823 | |
---|
824 | # Bed-slope and friction at vertices (and interpolated elsewhere) |
---|
825 | domain2.set_quantity('elevation', 0) |
---|
826 | domain2.set_quantity('friction', 0) |
---|
827 | domain2.set_quantity('stage', 0) |
---|
828 | |
---|
829 | # Read results for specific timesteps t=1 and t=2 |
---|
830 | from Scientific.IO.NetCDF import NetCDFFile |
---|
831 | fid = NetCDFFile(domain1.get_name() + '.sww') |
---|
832 | |
---|
833 | x = fid.variables['x'][:] |
---|
834 | y = fid.variables['y'][:] |
---|
835 | s1 = fid.variables['stage'][1,:] |
---|
836 | s2 = fid.variables['stage'][2,:] |
---|
837 | fid.close() |
---|
838 | |
---|
839 | shp = (len(x), 1) |
---|
840 | points = num.concatenate((num.reshape(x, shp), num.reshape(y, shp)), |
---|
841 | axis=1) |
---|
842 | #The diagonal points of domain 1 are 0, 5, 10, 15 |
---|
843 | |
---|
844 | assert num.allclose(num.take(points, [0,5,10,15], axis=0), |
---|
845 | [[0,0], [1.0/3,1.0/3], [2.0/3,2.0/3], [1,1]]) |
---|
846 | |
---|
847 | # Boundary conditions |
---|
848 | Br = Reflective_boundary(domain2) |
---|
849 | Bf = Field_boundary(domain1.get_name() + '.sww', |
---|
850 | domain2, mean_stage=1, verbose=False) |
---|
851 | |
---|
852 | domain2.set_boundary({'right':Br, 'bottom':Br, 'diagonal':Bf}) |
---|
853 | domain2.check_integrity() |
---|
854 | |
---|
855 | try: |
---|
856 | for t in domain2.evolve(yieldstep=1, finaltime=finaltime): |
---|
857 | pass |
---|
858 | except: |
---|
859 | pass |
---|
860 | else: |
---|
861 | msg = 'This should have caught NAN at boundary' |
---|
862 | raise Exception, msg |
---|
863 | |
---|
864 | #Cleanup |
---|
865 | os.remove(domain1.get_name() + '.sww') |
---|
866 | os.remove(domain2.get_name() + '.sww') |
---|
867 | |
---|
868 | |
---|
869 | |
---|
870 | ################################################################################# |
---|
871 | |
---|
872 | if __name__ == "__main__": |
---|
873 | suite = unittest.makeSuite(Test_swb_boundary_condition, 'test') |
---|
874 | runner = unittest.TextTestRunner(verbosity=1) |
---|
875 | runner.run(suite) |
---|