source: trunk/anuga_core/source/anuga/shallow_water_balanced/test_swb_boundary_condition.py @ 7869

Last change on this file since 7869 was 7869, checked in by hudson, 12 years ago

Shallow water balanced tests fail, but no longer have errors.

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