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

Last change on this file since 7573 was 7573, checked in by steve, 14 years ago

Committing a version of shallow_water_balanced which passes it unit test
using a version of edge limiting which doesn't limit boundary edges. THis
is useful to allow linear functions to be reconstructed.

Had to play with the transmissive BC and use centroid values which is
set via domain set_centroid_transmissive_bc

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.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
872if __name__ == "__main__":
873    suite = unittest.makeSuite(Test_swb_boundary_condition, 'test')
874    runner = unittest.TextTestRunner(verbosity=1)
875    runner.run(suite)
Note: See TracBrowser for help on using the repository browser.