source: anuga_core/source/anuga/shallow_water_balanced/test_swb_boundary_condition.py @ 7559

Last change on this file since 7559 was 7559, checked in by steve, 10 years ago

Broke up test_swb_domain so that we can concentrate on new functionality

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