source: trunk/anuga_core/source/anuga/shallow_water/test_sww_interrogate.py @ 7870

Last change on this file since 7870 was 7778, checked in by hudson, 14 years ago

Cleaned up unit tests to use API.

File size: 33.3 KB
Line 
1import unittest
2import copy
3import os
4import numpy as num
5
6from anuga.coordinate_transforms.geo_reference import Geo_reference
7from anuga.geometry.polygon import is_inside_polygon
8from anuga.abstract_2d_finite_volumes.util import file_function
9from anuga.config import netcdf_mode_r, netcdf_mode_w, netcdf_mode_a
10from anuga.config import g
11
12from boundaries import Reflective_boundary, \
13            Field_boundary, Transmissive_momentum_set_stage_boundary, \
14            Transmissive_stage_zero_momentum_boundary
15from anuga.abstract_2d_finite_volumes.generic_boundary_conditions\
16     import Transmissive_boundary, Dirichlet_boundary, \
17            Time_boundary, File_boundary, AWI_boundary
18
19from anuga.file.sww import get_mesh_and_quantities_from_file
20           
21from shallow_water_domain import Domain
22
23from anuga.abstract_2d_finite_volumes.mesh_factory \
24        import rectangular_cross
25from sww_interrogate import get_maximum_inundation_elevation, \
26            get_maximum_inundation_location, get_maximum_inundation_data, \
27            get_flow_through_cross_section, get_energy_through_cross_section
28           
29           
30               
31
32class Test_sww_Interrogate(unittest.TestCase):
33    def test_get_maximum_inundation(self):
34        """Test that sww information can be converted correctly to maximum
35        runup elevation and location (without and with georeferencing)
36
37        This test creates a slope and a runup which is maximal (~11m) at around 10s
38        and levels out to the boundary condition (1m) at about 30s.
39        """
40
41        import time, os
42        from Scientific.IO.NetCDF import NetCDFFile
43
44        #Setup
45
46        from mesh_factory import rectangular
47
48        # Create basic mesh (100m x 100m)
49        points, vertices, boundary = rectangular(20, 5, 100, 50)
50
51        # Create shallow water domain
52        domain = Domain(points, vertices, boundary)
53        domain.default_order = 2
54        domain.set_minimum_storable_height(0.01)
55
56        domain.set_name('runuptest')
57        swwfile = domain.get_name() + '.sww'
58
59        domain.set_datadir('.')
60        domain.format = 'sww'
61        domain.smooth = True
62
63        # FIXME (Ole): Backwards compatibility
64        # Look at sww file and see what happens when
65        # domain.tight_slope_limiters = 1
66        domain.tight_slope_limiters = 0
67        domain.use_centroid_velocities = 0 # Backwards compatibility (7/5/8)       
68       
69        Br = Reflective_boundary(domain)
70        Bd = Dirichlet_boundary([1.0,0,0])
71
72
73        #---------- First run without geo referencing
74       
75        domain.set_quantity('elevation', lambda x,y: -0.2*x + 14) # Slope
76        domain.set_quantity('stage', -6)
77        domain.set_boundary( {'left': Br, 'right': Bd, 'top': Br, 'bottom': Br})
78
79        for t in domain.evolve(yieldstep=1, finaltime = 50):
80            pass
81
82
83        # Check maximal runup
84        runup = get_maximum_inundation_elevation(swwfile)
85        location = get_maximum_inundation_location(swwfile)
86        #print 'Runup, location', runup, location
87        assert num.allclose(runup, 11) or num.allclose(runup, 12) # old limiters
88        assert num.allclose(location[0], 15) or num.allclose(location[0], 10)
89
90        # Check final runup
91        runup = get_maximum_inundation_elevation(swwfile, time_interval=[45,50])
92        location = get_maximum_inundation_location(swwfile, time_interval=[45,50])
93        # print 'Runup, location:',runup, location       
94        assert num.allclose(runup, 1)
95        assert num.allclose(location[0], 65)
96
97        # Check runup restricted to a polygon
98        p = [[50,1], [99,1], [99,49], [50,49]]
99        runup = get_maximum_inundation_elevation(swwfile, polygon=p)
100        location = get_maximum_inundation_location(swwfile, polygon=p)
101        #print runup, location       
102        assert num.allclose(runup, 4)
103        assert num.allclose(location[0], 50)               
104
105        # Check that mimimum_storable_height works
106        fid = NetCDFFile(swwfile, netcdf_mode_r) # Open existing file
107       
108        stage = fid.variables['stage'][:]
109        z = fid.variables['elevation'][:]
110        xmomentum = fid.variables['xmomentum'][:]
111        ymomentum = fid.variables['ymomentum'][:]       
112
113       
114       
115        for i in range(stage.shape[0]):
116            h = stage[i]-z # depth vector at time step i
117           
118            # Check every node location
119            for j in range(stage.shape[1]):
120                # Depth being either exactly zero implies
121                # momentum being zero.
122                # Or else depth must be greater than or equal to
123                # the minimal storable height
124                if h[j] == 0.0:
125                    assert xmomentum[i,j] == 0.0
126                    assert ymomentum[i,j] == 0.0               
127                else:
128                    assert h[j] >= domain.minimum_storable_height
129       
130        fid.close()
131
132        # Cleanup
133        os.remove(swwfile)
134       
135
136
137        #------------- Now the same with georeferencing
138
139        domain.time=0.0
140        E = 308500
141        N = 6189000
142        #E = N = 0
143        domain.geo_reference = Geo_reference(56, E, N)
144
145        domain.set_quantity('elevation', lambda x,y: -0.2*x + 14) # Slope
146        domain.set_quantity('stage', -6)
147        domain.set_boundary( {'left': Br, 'right': Bd, 'top': Br, 'bottom': Br})
148
149        for t in domain.evolve(yieldstep=1, finaltime = 50):
150            pass
151
152        # Check maximal runup
153        runup = get_maximum_inundation_elevation(swwfile)
154        location = get_maximum_inundation_location(swwfile)
155        assert num.allclose(runup, 11) or num.allclose(runup, 12) # old limiters
156        assert num.allclose(location[0], 15+E) or num.allclose(location[0], 10+E)
157
158        # Check final runup
159        runup = get_maximum_inundation_elevation(swwfile, time_interval=[45,50])
160        location = get_maximum_inundation_location(swwfile, time_interval=[45,50])
161        assert num.allclose(runup, 1)
162        assert num.allclose(location[0], 65+E)
163
164        # Check runup restricted to a polygon
165        p = num.array([[50,1], [99,1], [99,49], [50,49]], num.int) + num.array([E, N], num.int)      #array default#
166
167        runup = get_maximum_inundation_elevation(swwfile, polygon=p)
168        location = get_maximum_inundation_location(swwfile, polygon=p)
169        assert num.allclose(runup, 4)
170        assert num.allclose(location[0], 50+E)               
171
172
173        # Cleanup
174        os.remove(swwfile)
175
176
177
178    def test_get_flow_through_cross_section(self):
179        """test_get_flow_through_cross_section(self):
180
181        Test that the total flow through a cross section can be
182        correctly obtained from an sww file.
183       
184        This test creates a flat bed with a known flow through it and tests
185        that the function correctly returns the expected flow.
186
187        The specifics are
188        u = 2 m/s
189        h = 1 m
190        w = 3 m (width of channel)
191
192        q = u*h*w = 6 m^3/s
193
194        #---------- First run without geo referencing       
195       
196        """
197
198        import time, os
199        from Scientific.IO.NetCDF import NetCDFFile
200
201        # Setup
202        from mesh_factory import rectangular
203
204        # Create basic mesh (20m x 3m)
205        width = 3
206        length = 20
207        t_end = 3
208        points, vertices, boundary = rectangular(length, width,
209                                                 length, width)
210
211        # Create shallow water domain
212        domain = Domain(points, vertices, boundary)
213        domain.default_order = 2
214        domain.set_minimum_storable_height(0.01)
215
216        domain.set_name('flowtest')
217        swwfile = domain.get_name() + '.sww'
218
219        domain.set_datadir('.')
220        domain.format = 'sww'
221        domain.smooth = True
222
223        h = 1.0
224        u = 2.0
225        uh = u*h
226
227        Br = Reflective_boundary(domain)     # Side walls
228        Bd = Dirichlet_boundary([h, uh, 0])  # 2 m/s across the 3 m inlet:
229
230
231       
232        domain.set_quantity('elevation', 0.0)
233        domain.set_quantity('stage', h)
234        domain.set_quantity('xmomentum', uh)
235        domain.set_boundary( {'left': Bd, 'right': Bd, 'top': Br, 'bottom': Br})
236
237        for t in domain.evolve(yieldstep=1, finaltime = t_end):
238            pass
239
240        # Check that momentum is as it should be in the interior
241
242        I = [[0, width/2.],
243             [length/2., width/2.],
244             [length, width/2.]]
245       
246        f = file_function(swwfile,
247                          quantities=['stage', 'xmomentum', 'ymomentum'],
248                          interpolation_points=I,
249                          verbose=False)
250        for t in range(t_end+1):
251            for i in range(3):
252                assert num.allclose(f(t, i), [1, 2, 0], atol=1.0e-6)
253           
254
255        # Check flows through the middle
256        for i in range(5):
257            x = length/2. + i*0.23674563 # Arbitrary
258            cross_section = [[x, 0], [x, width]]
259            time, Q = get_flow_through_cross_section(swwfile,
260                                                     cross_section,
261                                                     verbose=False)
262
263            assert num.allclose(Q, uh*width)
264
265
266       
267        # Try the same with partial lines
268        x = length/2.
269        for i in range(5):
270            start_point = [length/2., i*width/5.]
271            #print start_point
272                           
273            cross_section = [start_point, [length/2., width]]
274            time, Q = get_flow_through_cross_section(swwfile,
275                                                     cross_section,
276                                                     verbose=False)
277
278            #print i, Q, (width-start_point[1])
279            assert num.allclose(Q, uh*(width-start_point[1]))
280
281
282        # Verify no flow when line is parallel to flow
283        cross_section = [[length/2.-10, width/2.], [length/2.+10, width/2.]]
284        time, Q = get_flow_through_cross_section(swwfile,
285                                                 cross_section,
286                                                 verbose=False)
287
288        #print i, Q
289        assert num.allclose(Q, 0, atol=1.0e-5)       
290
291
292        # Try with lines on an angle (all flow still runs through here)
293        cross_section = [[length/2., 0], [length/2.+width, width]]
294        time, Q = get_flow_through_cross_section(swwfile,
295                                                 cross_section,
296                                                 verbose=False)
297
298        assert num.allclose(Q, uh*width)       
299       
300
301
302                                     
303    def test_get_flow_through_cross_section_with_geo(self):
304        """test_get_flow_through_cross_section(self):
305
306        Test that the total flow through a cross section can be
307        correctly obtained from an sww file.
308       
309        This test creates a flat bed with a known flow through it and tests
310        that the function correctly returns the expected flow.
311
312        The specifics are
313        u = 2 m/s
314        h = 2 m
315        w = 3 m (width of channel)
316
317        q = u*h*w = 12 m^3/s
318
319
320        This run tries it with georeferencing and with elevation = -1
321       
322        """
323
324        import time, os
325        from Scientific.IO.NetCDF import NetCDFFile
326
327        # Setup
328        from mesh_factory import rectangular
329
330        # Create basic mesh (20m x 3m)
331        width = 3
332        length = 20
333        t_end = 1
334        points, vertices, boundary = rectangular(length, width,
335                                                 length, width)
336
337        # Create shallow water domain
338        domain = Domain(points, vertices, boundary,
339                        geo_reference = Geo_reference(56,308500,6189000))
340
341        domain.default_order = 2
342        domain.set_minimum_storable_height(0.01)
343
344        domain.set_name('flowtest')
345        swwfile = domain.get_name() + '.sww'
346
347        domain.set_datadir('.')
348        domain.format = 'sww'
349        domain.smooth = True
350
351        e = -1.0
352        w = 1.0
353        h = w-e
354        u = 2.0
355        uh = u*h
356
357        Br = Reflective_boundary(domain)     # Side walls
358        Bd = Dirichlet_boundary([w, uh, 0])  # 2 m/s across the 3 m inlet:
359
360
361
362       
363        domain.set_quantity('elevation', e)
364        domain.set_quantity('stage', w)
365        domain.set_quantity('xmomentum', uh)
366        domain.set_boundary( {'left': Bd, 'right': Bd, 'top': Br, 'bottom': Br})
367
368        for t in domain.evolve(yieldstep=1, finaltime = t_end):
369            pass
370
371        # Check that momentum is as it should be in the interior
372
373        I = [[0, width/2.],
374             [length/2., width/2.],
375             [length, width/2.]]
376       
377        I = domain.geo_reference.get_absolute(I)
378        f = file_function(swwfile,
379                          quantities=['stage', 'xmomentum', 'ymomentum'],
380                          interpolation_points=I,
381                          verbose=False)
382
383        for t in range(t_end+1):
384            for i in range(3):
385                #print i, t, f(t, i)           
386                assert num.allclose(f(t, i), [w, uh, 0], atol=1.0e-6)
387           
388
389        # Check flows through the middle
390        for i in range(5):
391            x = length/2. + i*0.23674563 # Arbitrary
392            cross_section = [[x, 0], [x, width]]
393
394            cross_section = domain.geo_reference.get_absolute(cross_section)           
395            time, Q = get_flow_through_cross_section(swwfile,
396                                                     cross_section,
397                                                     verbose=False)
398
399            assert num.allclose(Q, uh*width)
400
401
402           
403    def test_get_energy_through_cross_section(self):
404        """test_get_energy_through_cross_section(self):
405
406        Test that the specific and total energy through a cross section can be
407        correctly obtained from an sww file.
408       
409        This test creates a flat bed with a known flow through it and tests
410        that the function correctly returns the expected energies.
411
412        The specifics are
413        u = 2 m/s
414        h = 1 m
415        w = 3 m (width of channel)
416
417        q = u*h*w = 6 m^3/s
418        Es = h + 0.5*v*v/g  # Specific energy head [m]
419        Et = w + 0.5*v*v/g  # Total energy head [m]       
420
421
422        This test uses georeferencing
423       
424        """
425
426        import time, os
427        from Scientific.IO.NetCDF import NetCDFFile
428
429        # Setup
430        from mesh_factory import rectangular
431
432        # Create basic mesh (20m x 3m)
433        width = 3
434        length = 20
435        t_end = 1
436        points, vertices, boundary = rectangular(length, width,
437                                                 length, width)
438
439        # Create shallow water domain
440        domain = Domain(points, vertices, boundary,
441                        geo_reference = Geo_reference(56,308500,6189000))
442
443        domain.default_order = 2
444        domain.set_minimum_storable_height(0.01)
445
446        domain.set_name('flowtest')
447        swwfile = domain.get_name() + '.sww'
448
449        domain.set_datadir('.')
450        domain.format = 'sww'
451        domain.smooth = True
452
453        e = -1.0
454        w = 1.0
455        h = w-e
456        u = 2.0
457        uh = u*h
458
459        Br = Reflective_boundary(domain)     # Side walls
460        Bd = Dirichlet_boundary([w, uh, 0])  # 2 m/s across the 3 m inlet:
461
462       
463        domain.set_quantity('elevation', e)
464        domain.set_quantity('stage', w)
465        domain.set_quantity('xmomentum', uh)
466        domain.set_boundary( {'left': Bd, 'right': Bd, 'top': Br, 'bottom': Br})
467
468        for t in domain.evolve(yieldstep=1, finaltime = t_end):
469            pass
470
471        # Check that momentum is as it should be in the interior
472
473        I = [[0, width/2.],
474             [length/2., width/2.],
475             [length, width/2.]]
476       
477        I = domain.geo_reference.get_absolute(I)
478        f = file_function(swwfile,
479                          quantities=['stage', 'xmomentum', 'ymomentum'],
480                          interpolation_points=I,
481                          verbose=False)
482
483        for t in range(t_end+1):
484            for i in range(3):
485                #print i, t, f(t, i)
486                assert num.allclose(f(t, i), [w, uh, 0], atol=1.0e-6)
487           
488
489        # Check energies through the middle
490        for i in range(5):
491            x = length/2. + i*0.23674563 # Arbitrary
492            cross_section = [[x, 0], [x, width]]
493
494            cross_section = domain.geo_reference.get_absolute(cross_section)           
495           
496            time, Es = get_energy_through_cross_section(swwfile,
497                                                       cross_section,
498                                                       kind='specific',
499                                                       verbose=False)
500            assert num.allclose(Es, h + 0.5*u*u/g)
501           
502            time, Et = get_energy_through_cross_section(swwfile,
503                                                        cross_section,
504                                                        kind='total',
505                                                        verbose=False)
506            assert num.allclose(Et, w + 0.5*u*u/g)           
507
508
509
510    def test_get_maximum_inundation_from_sww(self):
511        """test_get_maximum_inundation_from_sww(self)
512
513        Test of get_maximum_inundation_elevation()
514        and get_maximum_inundation_location() from data_manager.py
515
516        This is based on test_get_maximum_inundation_3(self) but works with the
517        stored results instead of with the internal data structure.
518
519        This test uses the underlying get_maximum_inundation_data for tests
520        """
521
522
523        initial_runup_height = -0.4
524        final_runup_height = -0.3
525
526        #--------------------------------------------------------------
527        # Setup computational domain
528        #--------------------------------------------------------------
529        N = 10
530        points, vertices, boundary = rectangular_cross(N, N)
531        domain = Domain(points, vertices, boundary)
532        domain.set_name('runup_test')
533        domain.set_maximum_allowed_speed(1.0)
534
535        # FIXME: This works better with old limiters so far
536        domain.tight_slope_limiters = 0
537
538        #--------------------------------------------------------------
539        # Setup initial conditions
540        #--------------------------------------------------------------
541        def topography(x, y):
542            return -x/2                             # linear bed slope
543
544        # Use function for elevation
545        domain.set_quantity('elevation', topography)
546        domain.set_quantity('friction', 0.)                # Zero friction
547        # Constant negative initial stage
548        domain.set_quantity('stage', initial_runup_height)
549
550        #--------------------------------------------------------------
551        # Setup boundary conditions
552        #--------------------------------------------------------------
553        Br = Reflective_boundary(domain)                       # Reflective wall
554        Bd = Dirichlet_boundary([final_runup_height, 0, 0])    # Constant inflow
555
556        # All reflective to begin with (still water)
557        domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
558
559        #--------------------------------------------------------------
560        # Test initial inundation height
561        #--------------------------------------------------------------
562        indices = domain.get_wet_elements()
563        z = domain.get_quantity('elevation').\
564                get_values(location='centroids', indices=indices)
565        assert num.alltrue(z < initial_runup_height)
566
567        q_ref = domain.get_maximum_inundation_elevation()
568        # First order accuracy
569        assert num.allclose(q_ref, initial_runup_height, rtol=1.0/N)
570
571        #--------------------------------------------------------------
572        # Let triangles adjust
573        #--------------------------------------------------------------
574        for t in domain.evolve(yieldstep = 0.1, finaltime = 1.0):
575            pass
576
577        #--------------------------------------------------------------
578        # Test inundation height again
579        #--------------------------------------------------------------
580        q_ref = domain.get_maximum_inundation_elevation()
581        q = get_maximum_inundation_elevation('runup_test.sww')
582        msg = 'We got %f, should have been %f' % (q, q_ref)
583        assert num.allclose(q, q_ref, rtol=1.0/N), msg
584
585        q = get_maximum_inundation_elevation('runup_test.sww')
586        msg = 'We got %f, should have been %f' % (q, initial_runup_height)
587        assert num.allclose(q, initial_runup_height, rtol = 1.0/N), msg
588
589        # Test error condition if time interval is out
590        try:
591            q = get_maximum_inundation_elevation('runup_test.sww',
592                                                 time_interval=[2.0, 3.0])
593        except ValueError:
594            pass
595        else:
596            msg = 'should have caught wrong time interval'
597            raise Exception, msg
598
599        # Check correct time interval
600        q, loc = get_maximum_inundation_data('runup_test.sww',
601                                             time_interval=[0.0, 3.0])
602        msg = 'We got %f, should have been %f' % (q, initial_runup_height)
603        assert num.allclose(q, initial_runup_height, rtol = 1.0/N), msg
604        assert num.allclose(-loc[0]/2, q)    # From topography formula
605
606        #--------------------------------------------------------------
607        # Update boundary to allow inflow
608        #--------------------------------------------------------------
609        domain.set_boundary({'right': Bd})
610
611        #--------------------------------------------------------------
612        # Evolve system through time
613        #--------------------------------------------------------------
614        q_max = None
615        for t in domain.evolve(yieldstep = 0.1, finaltime = 3.0,
616                               skip_initial_step = True):
617            q = domain.get_maximum_inundation_elevation()
618            if q > q_max:
619                q_max = q
620
621        #--------------------------------------------------------------
622        # Test inundation height again
623        #--------------------------------------------------------------
624        indices = domain.get_wet_elements()
625        z = domain.get_quantity('elevation').\
626                get_values(location='centroids', indices=indices)
627
628        assert num.alltrue(z < final_runup_height)
629
630        q = domain.get_maximum_inundation_elevation()
631        # First order accuracy
632        assert num.allclose(q, final_runup_height, rtol=1.0/N)
633
634        q, loc = get_maximum_inundation_data('runup_test.sww',
635                                             time_interval=[3.0, 3.0])
636        msg = 'We got %f, should have been %f' % (q, final_runup_height)
637        assert num.allclose(q, final_runup_height, rtol=1.0/N), msg
638        assert num.allclose(-loc[0]/2, q)    # From topography formula
639
640        q = get_maximum_inundation_elevation('runup_test.sww')
641        loc = get_maximum_inundation_location('runup_test.sww')
642        msg = 'We got %f, should have been %f' % (q, q_max)
643        assert num.allclose(q, q_max, rtol=1.0/N), msg
644        assert num.allclose(-loc[0]/2, q)    # From topography formula
645
646        q = get_maximum_inundation_elevation('runup_test.sww',
647                                             time_interval=[0, 3])
648        msg = 'We got %f, should have been %f' % (q, q_max)
649        assert num.allclose(q, q_max, rtol=1.0/N), msg
650
651        # Check polygon mode
652        # Runup region
653        polygon = [[0.3, 0.0], [0.9, 0.0], [0.9, 1.0], [0.3, 1.0]]
654        q = get_maximum_inundation_elevation('runup_test.sww',
655                                             polygon = polygon,
656                                             time_interval=[0, 3])
657        msg = 'We got %f, should have been %f' % (q, q_max)
658        assert num.allclose(q, q_max, rtol=1.0/N), msg
659
660        # Offshore region
661        polygon = [[0.9, 0.0], [1.0, 0.0], [1.0, 1.0], [0.9, 1.0]]
662        q, loc = get_maximum_inundation_data('runup_test.sww',
663                                             polygon = polygon,
664                                             time_interval=[0, 3])
665        msg = 'We got %f, should have been %f' % (q, -0.475)
666        assert num.allclose(q, -0.475, rtol=1.0/N), msg
667        assert is_inside_polygon(loc, polygon)
668        assert num.allclose(-loc[0]/2, q)    # From topography formula
669
670        # Dry region
671        polygon = [[0.0, 0.0], [0.4, 0.0], [0.4, 1.0], [0.0, 1.0]]
672        q, loc = get_maximum_inundation_data('runup_test.sww',
673                                             polygon = polygon,
674                                             time_interval=[0, 3])
675        msg = 'We got %s, should have been None' % (q)
676        assert q is None, msg
677        msg = 'We got %s, should have been None' % (loc)
678        assert loc is None, msg
679
680        # Check what happens if no time point is within interval
681        try:
682            q = get_maximum_inundation_elevation('runup_test.sww',
683                                                 time_interval=[2.75, 2.75])
684        except AssertionError:
685            pass
686        else:
687            msg = 'Time interval should have raised an exception'
688            raise Exception, msg
689
690        # Cleanup
691        try:
692            os.remove(domain.get_name() + '.sww')
693        except:
694            pass
695            #FIXME(Ole): Windows won't allow removal of this
696
697
698    def test_get_maximum_inundation_from_sww(self):
699        """test_get_maximum_inundation_from_sww(self)
700
701        Test of get_maximum_inundation_elevation()
702        and get_maximum_inundation_location().
703   
704        This is based on test_get_maximum_inundation_3(self) but works with the
705        stored results instead of with the internal data structure.
706
707        This test uses the underlying get_maximum_inundation_data for tests
708        """
709
710        initial_runup_height = -0.4
711        final_runup_height = -0.3
712
713        #--------------------------------------------------------------
714        # Setup computational domain
715        #--------------------------------------------------------------
716        N = 10
717        points, vertices, boundary = rectangular_cross(N, N)
718        domain = Domain(points, vertices, boundary)
719        domain.set_name('runup_test')
720        domain.set_maximum_allowed_speed(1.0)
721
722        # FIXME: This works better with old limiters so far
723        domain.tight_slope_limiters = 0
724
725        #--------------------------------------------------------------
726        # Setup initial conditions
727        #--------------------------------------------------------------
728        def topography(x, y):
729            return -x/2                             # linear bed slope
730
731        # Use function for elevation
732        domain.set_quantity('elevation', topography)
733        domain.set_quantity('friction', 0.)                # Zero friction
734        # Constant negative initial stage
735        domain.set_quantity('stage', initial_runup_height)
736
737        #--------------------------------------------------------------
738        # Setup boundary conditions
739        #--------------------------------------------------------------
740        Br = Reflective_boundary(domain)                       # Reflective wall
741        Bd = Dirichlet_boundary([final_runup_height, 0, 0])    # Constant inflow
742
743        # All reflective to begin with (still water)
744        domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
745
746        #--------------------------------------------------------------
747        # Test initial inundation height
748        #--------------------------------------------------------------
749        indices = domain.get_wet_elements()
750        z = domain.get_quantity('elevation').\
751                get_values(location='centroids', indices=indices)
752        assert num.alltrue(z < initial_runup_height)
753
754        q_ref = domain.get_maximum_inundation_elevation()
755        # First order accuracy
756        assert num.allclose(q_ref, initial_runup_height, rtol=1.0/N)
757
758        #--------------------------------------------------------------
759        # Let triangles adjust
760        #--------------------------------------------------------------
761        for t in domain.evolve(yieldstep = 0.1, finaltime = 1.0):
762            pass
763
764        #--------------------------------------------------------------
765        # Test inundation height again
766        #--------------------------------------------------------------
767        q_ref = domain.get_maximum_inundation_elevation()
768        q = get_maximum_inundation_elevation('runup_test.sww')
769        msg = 'We got %f, should have been %f' % (q, q_ref)
770        assert num.allclose(q, q_ref, rtol=1.0/N), msg
771
772        q = get_maximum_inundation_elevation('runup_test.sww')
773        msg = 'We got %f, should have been %f' % (q, initial_runup_height)
774        assert num.allclose(q, initial_runup_height, rtol = 1.0/N), msg
775
776        # Test error condition if time interval is out
777        try:
778            q = get_maximum_inundation_elevation('runup_test.sww',
779                                                 time_interval=[2.0, 3.0])
780        except ValueError:
781            pass
782        else:
783            msg = 'should have caught wrong time interval'
784            raise Exception, msg
785
786        # Check correct time interval
787        q, loc = get_maximum_inundation_data('runup_test.sww',
788                                             time_interval=[0.0, 3.0])
789        msg = 'We got %f, should have been %f' % (q, initial_runup_height)
790        assert num.allclose(q, initial_runup_height, rtol = 1.0/N), msg
791        assert num.allclose(-loc[0]/2, q)    # From topography formula
792
793        #--------------------------------------------------------------
794        # Update boundary to allow inflow
795        #--------------------------------------------------------------
796        domain.set_boundary({'right': Bd})
797
798        #--------------------------------------------------------------
799        # Evolve system through time
800        #--------------------------------------------------------------
801        q_max = None
802        for t in domain.evolve(yieldstep = 0.1, finaltime = 3.0,
803                               skip_initial_step = True):
804            q = domain.get_maximum_inundation_elevation()
805            if q > q_max:
806                q_max = q
807
808        #--------------------------------------------------------------
809        # Test inundation height again
810        #--------------------------------------------------------------
811        indices = domain.get_wet_elements()
812        z = domain.get_quantity('elevation').\
813                get_values(location='centroids', indices=indices)
814
815        assert num.alltrue(z < final_runup_height)
816
817        q = domain.get_maximum_inundation_elevation()
818        # First order accuracy
819        assert num.allclose(q, final_runup_height, rtol=1.0/N)
820
821        q, loc = get_maximum_inundation_data('runup_test.sww',
822                                             time_interval=[3.0, 3.0])
823        msg = 'We got %f, should have been %f' % (q, final_runup_height)
824        assert num.allclose(q, final_runup_height, rtol=1.0/N), msg
825        assert num.allclose(-loc[0]/2, q)    # From topography formula
826
827        q = get_maximum_inundation_elevation('runup_test.sww')
828        loc = get_maximum_inundation_location('runup_test.sww')
829        msg = 'We got %f, should have been %f' % (q, q_max)
830        assert num.allclose(q, q_max, rtol=1.0/N), msg
831        assert num.allclose(-loc[0]/2, q)    # From topography formula
832
833        q = get_maximum_inundation_elevation('runup_test.sww',
834                                             time_interval=[0, 3])
835        msg = 'We got %f, should have been %f' % (q, q_max)
836        assert num.allclose(q, q_max, rtol=1.0/N), msg
837
838        # Check polygon mode
839        # Runup region
840        polygon = [[0.3, 0.0], [0.9, 0.0], [0.9, 1.0], [0.3, 1.0]]
841        q = get_maximum_inundation_elevation('runup_test.sww',
842                                             polygon = polygon,
843                                             time_interval=[0, 3])
844        msg = 'We got %f, should have been %f' % (q, q_max)
845        assert num.allclose(q, q_max, rtol=1.0/N), msg
846
847        # Offshore region
848        polygon = [[0.9, 0.0], [1.0, 0.0], [1.0, 1.0], [0.9, 1.0]]
849        q, loc = get_maximum_inundation_data('runup_test.sww',
850                                             polygon = polygon,
851                                             time_interval=[0, 3])
852        msg = 'We got %f, should have been %f' % (q, -0.475)
853        assert num.allclose(q, -0.475, rtol=1.0/N), msg
854        assert is_inside_polygon(loc, polygon)
855        assert num.allclose(-loc[0]/2, q)    # From topography formula
856
857        # Dry region
858        polygon = [[0.0, 0.0], [0.4, 0.0], [0.4, 1.0], [0.0, 1.0]]
859        q, loc = get_maximum_inundation_data('runup_test.sww',
860                                             polygon = polygon,
861                                             time_interval=[0, 3])
862        msg = 'We got %s, should have been None' % (q)
863        assert q is None, msg
864        msg = 'We got %s, should have been None' % (loc)
865        assert loc is None, msg
866
867        # Check what happens if no time point is within interval
868        try:
869            q = get_maximum_inundation_elevation('runup_test.sww',
870                                                 time_interval=[2.75, 2.75])
871        except AssertionError:
872            pass
873        else:
874            msg = 'Time interval should have raised an exception'
875            raise Exception, msg
876
877        # Cleanup
878        try:
879            os.remove(domain.get_name() + '.sww')
880        except:
881            pass
882            #FIXME(Ole): Windows won't allow removal of this
883
884 
885 
886 
887if __name__ == "__main__":
888    suite = unittest.makeSuite(Test_sww_Interrogate, 'test')
889    runner = unittest.TextTestRunner() #verbosity=2)
890    runner.run(suite)
891               
Note: See TracBrowser for help on using the repository browser.