source: trunk/anuga_core/source/anuga/shallow_water_balanced/test_swb_basic.py @ 7876

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

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

File size: 92.0 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# For test_fitting_using_shallow_water_domain example
29def linear_function(point):
30    point = num.array(point)
31    return point[:,0]+point[:,1]
32
33class Weir:
34    """Set a bathymetry for weir with a hole and a downstream gutter
35    x,y are assumed to be in the unit square
36    """
37
38    def __init__(self, stage):
39        self.inflow_stage = stage
40
41    def __call__(self, x, y):
42        N = len(x)
43        assert N == len(y)
44
45        z = num.zeros(N, num.float)
46        for i in range(N):
47            z[i] = -x[i]/2  #General slope
48
49            #Flattish bit to the left
50            if x[i] < 0.3:
51                z[i] = -x[i]/10
52
53            #Weir
54            if x[i] >= 0.3 and x[i] < 0.4:
55                z[i] = -x[i]+0.9
56
57            #Dip
58            x0 = 0.6
59            depth = -1.0
60            plateaux = -0.6
61            if y[i] < 0.7:
62                if x[i] > x0 and x[i] < 0.9:
63                    z[i] = depth
64                #RHS plateaux
65                if x[i] >= 0.9:
66                    z[i] = plateaux
67            elif y[i] >= 0.7 and y[i] < 1.5:
68                #Restrict and deepen
69                if x[i] >= x0 and x[i] < 0.8:
70                    z[i] = depth - (y[i]/3 - 0.3)
71                elif x[i] >= 0.8:
72                    #RHS plateaux
73                    z[i] = plateaux
74            elif y[i] >= 1.5:
75                if x[i] >= x0 and x[i] < 0.8 + (y[i]-1.5)/1.2:
76                    #Widen up and stay at constant depth
77                    z[i] = depth-1.5/5
78                elif x[i] >= 0.8 + (y[i]-1.5)/1.2:
79                    #RHS plateaux
80                    z[i] = plateaux
81
82            #Hole in weir (slightly higher than inflow condition)
83            if x[i] >= 0.3 and x[i] < 0.4 and y[i] > 0.2 and y[i] < 0.4:
84                z[i] = -x[i]+self.inflow_stage + 0.02
85
86            #Channel behind weir
87            x0 = 0.5
88            if x[i] >= 0.4 and x[i] < x0 and y[i] > 0.2 and y[i] < 0.4:
89                z[i] = -x[i]+self.inflow_stage + 0.02
90
91            if x[i] >= x0 and x[i] < 0.6 and y[i] > 0.2 and y[i] < 0.4:
92                #Flatten it out towards the end
93                z[i] = -x0+self.inflow_stage + 0.02 + (x0-x[i])/5
94
95            # Hole to the east
96            x0 = 1.1
97            y0 = 0.35
98            if num.sqrt((2*(x[i]-x0))**2 + (2*(y[i]-y0))**2) < 0.2:
99                z[i] = num.sqrt(((x[i]-x0))**2 + ((y[i]-y0))**2)-1.0
100
101            #Tiny channel draining hole
102            if x[i] >= 1.14 and x[i] < 1.2 and y[i] >= 0.4 and y[i] < 0.6:
103                z[i] = -0.9 #North south
104
105            if x[i] >= 0.9 and x[i] < 1.18 and y[i] >= 0.58 and y[i] < 0.65:
106                z[i] = -1.0 + (x[i]-0.9)/3 #East west
107
108            # Stuff not in use
109
110            # Upward slope at inlet to the north west
111            # if x[i] < 0.0: # and y[i] > 0.5:
112            #    #z[i] = -y[i]+0.5  #-x[i]/2
113            #    z[i] = x[i]/4 - y[i]**2 + 0.5
114
115            # Hole to the west
116            # x0 = -0.4; y0 = 0.35 # center
117            # if sqrt((2*(x[i]-x0))**2 + (2*(y[i]-y0))**2) < 0.2:
118            #    z[i] = sqrt(((x[i]-x0))**2 + ((y[i]-y0))**2)-0.2
119
120        return z/2
121
122
123
124#########################################################
125
126class Test_swb_basic(unittest.TestCase):
127    def setUp(self):
128        pass
129
130    def tearDown(self):
131        pass
132
133
134    def test_rotate(self):
135        from anuga.shallow_water.shallow_water_ext import rotate
136        normal = num.array([0.0, -1.0])
137
138        q = num.array([1.0, 2.0, 3.0])
139
140        r = rotate(q, normal, direction = 1)
141        assert r[0] == 1
142        assert r[1] == -3
143        assert r[2] == 2
144
145        w = rotate(r, normal, direction = -1)
146        assert num.allclose(w, q)
147
148        # Check error check
149        try:
150            rotate(r, num.array([1, 1, 1]))
151        except:
152            pass
153        else:
154            raise Exception, 'Should have raised an exception'
155
156    # Individual flux tests
157    def test_flux_zero_case(self):
158        ql = num.zeros(3, num.float)
159        qr = num.zeros(3, num.float)
160        normal = num.zeros(2, num.float)
161        edgeflux = num.zeros(3, num.float)
162        zl = zr = 0.
163        H0 = 1.0e-3 # As suggested in the manual
164       
165        max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux, epsilon, g, H0)
166
167        assert num.allclose(edgeflux, [0,0,0])
168        assert max_speed == 0.
169
170    def test_flux_constants(self):
171        w = 2.0
172
173        normal = num.array([1.,0])
174        ql = num.array([w, 0, 0])
175        qr = num.array([w, 0, 0])
176        edgeflux = num.zeros(3, num.float)
177        zl = zr = 0.
178        h = w - (zl+zr)/2
179        H0 = 0.0
180
181        max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux, epsilon, g, H0)       
182        assert num.allclose(edgeflux, [0., 0.5*g*h**2, 0.])
183        assert max_speed == num.sqrt(g*h)
184
185    #def test_flux_slope(self):
186    #    #FIXME: TODO
187    #    w = 2.0
188    #
189    #    normal = array([1.,0])
190    #    ql = array([w, 0, 0])
191    #    qr = array([w, 0, 0])
192    #    zl = zr = 0.
193    #    h = w - (zl+zr)/2
194    #
195    #    flux, max_speed = flux_function(normal, ql, qr, zl, zr)
196    #
197    #    assert allclose(flux, [0., 0.5*g*h**2, 0.])
198    #    assert max_speed == sqrt(g*h)
199
200    def test_flux1(self):
201        # Use data from previous version of abstract_2d_finite_volumes
202        normal = num.array([1., 0])
203        ql = num.array([-0.2, 2, 3])
204        qr = num.array([-0.2, 2, 3])
205        zl = zr = -0.5
206        edgeflux = num.zeros(3, num.float)
207
208        H0 = 0.0
209
210        max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux,
211                                  epsilon, g, H0)
212
213        assert num.allclose(edgeflux, [2., 13.77433333, 20.])
214        assert num.allclose(max_speed, 8.38130948661)
215
216    def test_flux2(self):
217        # Use data from previous version of abstract_2d_finite_volumes
218        normal = num.array([0., -1.])
219        ql = num.array([-0.075, 2, 3])
220        qr = num.array([-0.075, 2, 3])
221        zl = zr = -0.375
222
223        edgeflux = num.zeros(3, num.float)
224        H0 = 0.0
225        max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux,
226                                  epsilon, g, H0)
227
228        assert num.allclose(edgeflux, [-3., -20.0, -30.441])
229        assert num.allclose(max_speed, 11.7146428199)
230
231    def test_flux3(self):
232        # Use data from previous version of abstract_2d_finite_volumes
233        normal = num.array([-sqrt(2)/2, sqrt(2)/2])
234        ql = num.array([-0.075, 2, 3])
235        qr = num.array([-0.075, 2, 3])
236        zl = zr = -0.375
237
238        edgeflux = num.zeros(3, num.float)
239        H0 = 0.0
240        max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux,
241                                  epsilon, g, H0)
242
243        assert num.allclose(edgeflux, [sqrt(2)/2, 4.40221112, 7.3829019])
244        assert num.allclose(max_speed, 4.0716654239)
245
246    def test_flux4(self):
247        # Use data from previous version of abstract_2d_finite_volumes
248        normal = num.array([-sqrt(2)/2, sqrt(2)/2])
249        ql = num.array([-0.34319278, 0.10254161, 0.07273855])
250        qr = num.array([-0.30683287, 0.1071986, 0.05930515])
251        zl = zr = -0.375
252
253        edgeflux = num.zeros(3, num.float)
254        H0 = 0.0
255        max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux,
256                                  epsilon, g, H0)
257
258        assert num.allclose(edgeflux, [-0.04072676, -0.07096636, -0.01604364])
259        assert num.allclose(max_speed, 1.31414103233)
260
261    def test_flux_computation(self):
262        """test flux calculation (actual C implementation)
263
264        This one tests the constant case where only the pressure term
265        contributes to each edge and cancels out once the total flux has
266        been summed up.
267        """
268
269        a = [0.0, 0.0]
270        b = [0.0, 2.0]
271        c = [2.0, 0.0]
272        d = [0.0, 4.0]
273        e = [2.0, 2.0]
274        f = [4.0, 0.0]
275
276        points = [a, b, c, d, e, f]
277        #              bac,     bce,     ecf,     dbe
278        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
279
280        domain = Domain(points, vertices)
281        domain.check_integrity()
282
283        # The constant case
284        domain.set_quantity('elevation', -1)
285        domain.set_quantity('stage', 1)
286
287        domain.compute_fluxes()
288        # Central triangle
289        assert num.allclose(domain.get_quantity('stage').explicit_update[1], 0)
290
291        # The more general case
292        def surface(x, y):
293            return -x/2
294
295        domain.set_quantity('elevation', -10)
296        domain.set_quantity('stage', surface)
297        domain.set_quantity('xmomentum', 1)
298
299        domain.compute_fluxes()
300
301        #print domain.get_quantity('stage').explicit_update
302        # FIXME (Ole): TODO the general case
303        #assert allclose(domain.get_quantity('stage').explicit_update[1], ...??)
304
305    def test_sw_domain_simple(self):
306        a = [0.0, 0.0]
307        b = [0.0, 2.0]
308        c = [2.0, 0.0]
309        d = [0.0, 4.0]
310        e = [2.0, 2.0]
311        f = [4.0, 0.0]
312
313        points = [a, b, c, d, e, f]
314        #             bac,     bce,     ecf,     dbe
315        vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
316
317        #from anuga.abstract_2d_finite_volumes.domain import Domain as Generic_domain
318        #msg = 'The class %s is not a subclass of the generic domain class %s'\
319        #      %(DomainClass, Domain)
320        #assert issubclass(DomainClass, Domain), msg
321
322        domain = Domain(points, vertices)
323        domain.check_integrity()
324
325        for name in ['stage', 'xmomentum', 'ymomentum',
326                     'elevation', 'friction']:
327            assert domain.quantities.has_key(name)
328
329        assert num.alltrue(domain.get_conserved_quantities(0, edge=1) == 0.)
330
331    def xtest_catching_negative_heights(self):
332        #OBSOLETE
333
334        a = [0.0, 0.0]
335        b = [0.0, 2.0]
336        c = [2.0, 0.0]
337        d = [0.0, 4.0]
338        e = [2.0, 2.0]
339        f = [4.0, 0.0]
340
341        points = [a, b, c, d, e, f]
342        #             bac,     bce,     ecf,     dbe
343        vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
344
345        domain = Domain(points, vertices)
346        val0 = 2. + 2.0/3
347        val1 = 4. + 4.0/3
348        val2 = 8. + 2.0/3
349        val3 = 2. + 8.0/3
350
351        zl = zr = 4    # Too large
352        domain.set_quantity('elevation', zl*num.ones((4, 3), num.int)) #array default#
353        domain.set_quantity('stage', [[val0, val0-1, val0-2],
354                                      [val1, val1+1, val1],
355                                      [val2, val2-2, val2],
356                                      [val3-0.5, val3, val3]])
357
358        #Should fail
359        try:
360            domain.check_integrity()
361        except:
362            pass
363
364    def test_get_wet_elements(self):
365        a = [0.0, 0.0]
366        b = [0.0, 2.0]
367        c = [2.0, 0.0]
368        d = [0.0, 4.0]
369        e = [2.0, 2.0]
370        f = [4.0, 0.0]
371
372        points = [a, b, c, d, e, f]
373        #             bac,     bce,     ecf,     dbe
374        vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
375
376        domain = Domain(points, vertices)
377
378        val0 = 2. + 2.0/3
379        val1 = 4. + 4.0/3
380        val2 = 8. + 2.0/3
381        val3 = 2. + 8.0/3
382
383        zl = zr = 5
384        domain.set_quantity('elevation', zl*num.ones((4, 3), num.int)) #array default#
385        domain.set_quantity('stage', [[val0, val0-1, val0-2],
386                                      [val1, val1+1, val1],
387                                      [val2, val2-2, val2],
388                                      [val3-0.5, val3, val3]])
389
390        domain.check_integrity()
391
392        indices = domain.get_wet_elements()
393        assert num.allclose(indices, [1, 2])
394
395        indices = domain.get_wet_elements(indices=[0, 1, 3])
396        assert num.allclose(indices, [1])
397
398    def test_get_maximum_inundation_1(self):
399        a = [0.0, 0.0]
400        b = [0.0, 2.0]
401        c = [2.0, 0.0]
402        d = [0.0, 4.0]
403        e = [2.0, 2.0]
404        f = [4.0, 0.0]
405
406        points = [a, b, c, d, e, f]
407        #             bac,     bce,     ecf,     dbe
408        vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
409
410        domain = Domain(points, vertices)
411
412        domain.set_quantity('elevation', lambda x, y: x+2*y)    # 2 4 4 6
413        domain.set_quantity('stage', 3)
414
415        domain.check_integrity()
416
417        indices = domain.get_wet_elements()
418        assert num.allclose(indices, [0])
419
420        q = domain.get_maximum_inundation_elevation()
421        assert num.allclose(q, domain.get_quantity('elevation').\
422                                   get_values(location='centroids')[0])
423
424        x, y = domain.get_maximum_inundation_location()
425        assert num.allclose([x, y], domain.get_centroid_coordinates()[0])
426
427    def test_get_maximum_inundation_2(self):
428        """test_get_maximum_inundation_2(self)
429
430        Test multiple wet cells with same elevation
431        """
432
433        a = [0.0, 0.0]
434        b = [0.0, 2.0]
435        c = [2.0, 0.0]
436        d = [0.0, 4.0]
437        e = [2.0, 2.0]
438        f = [4.0, 0.0]
439
440        points = [a, b, c, d, e, f]
441        #             bac,     bce,     ecf,     dbe
442        vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
443
444        domain = Domain(points, vertices)
445
446        domain.set_quantity('elevation', lambda x, y: x+2*y)    # 2 4 4 6
447        domain.set_quantity('stage', 4.1)
448
449        domain.check_integrity()
450
451        indices = domain.get_wet_elements()
452        assert num.allclose(indices, [0, 1, 2])
453
454        q = domain.get_maximum_inundation_elevation()
455        assert num.allclose(q, 4)
456
457        x, y = domain.get_maximum_inundation_location()
458        assert num.allclose([x, y], domain.get_centroid_coordinates()[1])
459
460    def test_get_maximum_inundation_3(self):
461        """test_get_maximum_inundation_3(self)
462
463        Test of real runup example:
464        """
465
466        from anuga.abstract_2d_finite_volumes.mesh_factory \
467                import rectangular_cross
468
469        initial_runup_height = -0.4
470        final_runup_height = -0.3
471
472        #--------------------------------------------------------------
473        # Setup computational domain
474        #--------------------------------------------------------------
475        N = 5
476        points, vertices, boundary = rectangular_cross(N, N)
477        domain = Domain(points, vertices, boundary)
478        domain.set_maximum_allowed_speed(1.0)
479
480        #--------------------------------------------------------------
481        # Setup initial conditions
482        #--------------------------------------------------------------
483        def topography(x, y):
484            return -x/2                             # linear bed slope
485
486        # Use function for elevation
487        domain.set_quantity('elevation', topography)
488        domain.set_quantity('friction', 0.)                # Zero friction
489        # Constant negative initial stage
490        domain.set_quantity('stage', initial_runup_height)
491
492        #--------------------------------------------------------------
493        # Setup boundary conditions
494        #--------------------------------------------------------------
495        Br = anuga.Reflective_boundary(domain)                 # Reflective wall
496        Bd = anuga.Dirichlet_boundary([final_runup_height, 0, 0])# Steady inflow
497
498        # All reflective to begin with (still water)
499        domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
500
501        #--------------------------------------------------------------
502        # Test initial inundation height
503        #--------------------------------------------------------------
504
505        indices = domain.get_wet_elements()
506        z = domain.get_quantity('elevation').get_values(location='centroids',
507                                                        indices=indices)
508        assert num.alltrue(z < initial_runup_height)
509
510        q = domain.get_maximum_inundation_elevation()
511        # First order accuracy
512        assert num.allclose(q, initial_runup_height, rtol=1.0/N)
513
514        x, y = domain.get_maximum_inundation_location()
515
516        qref = domain.get_quantity('elevation').\
517                     get_values(interpolation_points=[[x, y]])
518        assert num.allclose(q, qref)
519
520        wet_elements = domain.get_wet_elements()
521        wet_elevations = domain.get_quantity('elevation').\
522                                    get_values(location='centroids',
523                                               indices=wet_elements)
524        assert num.alltrue(wet_elevations < initial_runup_height)
525        assert num.allclose(wet_elevations, z)
526
527        #--------------------------------------------------------------
528        # Let triangles adjust
529        #--------------------------------------------------------------
530        for t in domain.evolve(yieldstep = 0.1, finaltime = 1.0):
531            pass
532
533        #--------------------------------------------------------------
534        # Test inundation height again
535        #--------------------------------------------------------------
536        indices = domain.get_wet_elements()
537        z = domain.get_quantity('elevation').get_values(location='centroids',
538                                                        indices=indices)
539
540        assert num.alltrue(z < initial_runup_height)
541
542        q = domain.get_maximum_inundation_elevation()
543        # First order accuracy
544        assert num.allclose(q, initial_runup_height, rtol=1.0/N)
545
546        x, y = domain.get_maximum_inundation_location()
547        qref = domain.get_quantity('elevation').\
548                        get_values(interpolation_points=[[x, y]])
549        assert num.allclose(q, qref)
550
551        #--------------------------------------------------------------
552        # Update boundary to allow inflow
553        #--------------------------------------------------------------
554        domain.set_boundary({'right': Bd})
555
556        #--------------------------------------------------------------
557        # Evolve system through time
558        #--------------------------------------------------------------
559        for t in domain.evolve(yieldstep = 0.1, finaltime = 3.0):
560            pass
561
562        #--------------------------------------------------------------
563        # Test inundation height again
564        #--------------------------------------------------------------
565        indices = domain.get_wet_elements()
566        z = domain.get_quantity('elevation').\
567                    get_values(location='centroids', indices=indices)
568
569        assert num.alltrue(z < final_runup_height)
570
571        q = domain.get_maximum_inundation_elevation()
572        # First order accuracy
573        assert num.allclose(q, final_runup_height, rtol=1.0/N)
574
575        x, y = domain.get_maximum_inundation_location()
576        qref = domain.get_quantity('elevation').\
577                        get_values(interpolation_points=[[x, y]])
578        assert num.allclose(q, qref)
579
580        wet_elements = domain.get_wet_elements()
581        wet_elevations = domain.get_quantity('elevation').\
582                             get_values(location='centroids',
583                                        indices=wet_elements)
584        assert num.alltrue(wet_elevations < final_runup_height)
585        assert num.allclose(wet_elevations, z)
586
587    def test_get_maximum_inundation_from_sww(self):
588        """test_get_maximum_inundation_from_sww(self)
589
590        Test of get_maximum_inundation_elevation()
591        and get_maximum_inundation_location() from data_manager.py
592
593        This is based on test_get_maximum_inundation_3(self) but works with the
594        stored results instead of with the internal data structure.
595
596        This test uses the underlying get_maximum_inundation_data for tests
597        """
598
599        from anuga.abstract_2d_finite_volumes.mesh_factory \
600                import rectangular_cross
601        from anuga.shallow_water.sww_interrogate import \
602            get_maximum_inundation_elevation, get_maximum_inundation_location, \
603            get_maximum_inundation_data
604
605        initial_runup_height = -0.4
606        final_runup_height = -0.3
607
608        #--------------------------------------------------------------
609        # Setup computational domain
610        #--------------------------------------------------------------
611        N = 10
612        points, vertices, boundary = rectangular_cross(N, N)
613        domain = Domain(points, vertices, boundary)
614        domain.set_name('runup_test')
615        #domain.set_maximum_allowed_speed(1.0)
616
617        # FIXME: This works better with old limiters so far
618        #domain.tight_slope_limiters = 0
619
620        #--------------------------------------------------------------
621        # Setup initial conditions
622        #--------------------------------------------------------------
623        def topography(x, y):
624            return -x/2                             # linear bed slope
625
626        # Use function for elevation
627        domain.set_quantity('elevation', topography)
628        domain.set_quantity('friction', 0.)                # Zero friction
629        # Constant negative initial stage
630        domain.set_quantity('stage', initial_runup_height)
631
632        #--------------------------------------------------------------
633        # Setup boundary conditions
634        #--------------------------------------------------------------
635        Br = anuga.Reflective_boundary(domain)                 # Reflective wall
636        Bd = anuga.Dirichlet_boundary([final_runup_height, 0, 0])# Steady inflow
637
638        # All reflective to begin with (still water)
639        domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
640
641        #--------------------------------------------------------------
642        # Test initial inundation height
643        #--------------------------------------------------------------
644        indices = domain.get_wet_elements()
645        z = domain.get_quantity('elevation').\
646                get_values(location='centroids', indices=indices)
647        assert num.alltrue(z < initial_runup_height)
648
649        q_ref = domain.get_maximum_inundation_elevation()
650        # First order accuracy
651        assert num.allclose(q_ref, initial_runup_height, rtol=1.0/N)
652
653        #--------------------------------------------------------------
654        # Let triangles adjust
655        #--------------------------------------------------------------
656        for t in domain.evolve(yieldstep = 0.1, finaltime = 1.0):
657            pass
658
659        #--------------------------------------------------------------
660        # Test inundation height again
661        #--------------------------------------------------------------
662        q_ref = domain.get_maximum_inundation_elevation()
663        q = get_maximum_inundation_elevation('runup_test.sww')
664        msg = 'We got %f, should have been %f' % (q, q_ref)
665        assert num.allclose(q, q_ref, rtol=1.0/N), msg
666
667        q = get_maximum_inundation_elevation('runup_test.sww')
668        msg = 'We got %f, should have been %f' % (q, initial_runup_height)
669        assert num.allclose(q, initial_runup_height, rtol = 1.0/N), msg
670
671        # Test error condition if time interval is out
672        try:
673            q = get_maximum_inundation_elevation('runup_test.sww',
674                                                 time_interval=[2.0, 3.0])
675        except ValueError:
676            pass
677        else:
678            msg = 'should have caught wrong time interval'
679            raise Exception, msg
680
681        # Check correct time interval
682        q, loc = get_maximum_inundation_data('runup_test.sww',
683                                             time_interval=[0.0, 3.0])
684        msg = 'We got %f, should have been %f' % (q, initial_runup_height)
685        assert num.allclose(q, initial_runup_height, rtol = 1.0/N), msg
686        assert num.allclose(-loc[0]/2, q)    # From topography formula
687
688        #--------------------------------------------------------------
689        # Update boundary to allow inflow
690        #--------------------------------------------------------------
691        domain.set_boundary({'right': Bd})
692
693        #--------------------------------------------------------------
694        # Evolve system through time
695        #--------------------------------------------------------------
696        q_max = None
697        for t in domain.evolve(yieldstep = 0.1, finaltime = 3.0,
698                               skip_initial_step = True):
699            q = domain.get_maximum_inundation_elevation()
700            if q > q_max:
701                q_max = q
702
703        #--------------------------------------------------------------
704        # Test inundation height again
705        #--------------------------------------------------------------
706        indices = domain.get_wet_elements()
707        z = domain.get_quantity('elevation').\
708                get_values(location='centroids', indices=indices)
709
710        assert num.alltrue(z < final_runup_height)
711
712        q = domain.get_maximum_inundation_elevation()
713        # First order accuracy
714        assert num.allclose(q, final_runup_height, rtol=1.0/N)
715
716        q, loc = get_maximum_inundation_data('runup_test.sww',
717                                             time_interval=[3.0, 3.0])
718        msg = 'We got %f, should have been %f' % (q, final_runup_height)
719        assert num.allclose(q, final_runup_height, rtol=1.0/N), msg
720        assert num.allclose(-loc[0]/2, q)    # From topography formula
721
722        q = get_maximum_inundation_elevation('runup_test.sww')
723        loc = get_maximum_inundation_location('runup_test.sww')
724        msg = 'We got %f, should have been %f' % (q, q_max)
725        assert num.allclose(q, q_max, rtol=1.0/N), msg
726        assert num.allclose(-loc[0]/2, q)    # From topography formula
727
728        q = get_maximum_inundation_elevation('runup_test.sww',
729                                             time_interval=[0, 3])
730        msg = 'We got %f, should have been %f' % (q, q_max)
731        assert num.allclose(q, q_max, rtol=1.0/N), msg
732
733        # Check polygon mode
734        # Runup region
735        polygon = [[0.3, 0.0], [0.9, 0.0], [0.9, 1.0], [0.3, 1.0]]
736        q = get_maximum_inundation_elevation('runup_test.sww',
737                                             polygon = polygon,
738                                             time_interval=[0, 3])
739        msg = 'We got %f, should have been %f' % (q, q_max)
740        assert num.allclose(q, q_max, rtol=1.0/N), msg
741
742        # Offshore region
743        polygon = [[0.9, 0.0], [1.0, 0.0], [1.0, 1.0], [0.9, 1.0]]
744        q, loc = get_maximum_inundation_data('runup_test.sww',
745                                             polygon = polygon,
746                                             time_interval=[0, 3])
747        msg = 'We got %f, should have been %f' % (q, -0.475)
748        assert num.allclose(q, -0.475, rtol=1.0/N), msg
749        assert is_inside_polygon(loc, polygon)
750        assert num.allclose(-loc[0]/2, q)    # From topography formula
751
752        # Dry region
753        polygon = [[0.0, 0.0], [0.2, 0.0], [0.2, 1.0], [0.0, 1.0]]
754        q, loc = get_maximum_inundation_data('runup_test.sww',
755                                             polygon = polygon,
756                                             time_interval=[0, 3])
757        msg = 'We got %s, should have been None' % (q)
758        assert q is None, msg
759        msg = 'We got %s, should have been None' % (loc)
760        assert loc is None, msg
761
762        # Check what happens if no time point is within interval
763        try:
764            q = get_maximum_inundation_elevation('runup_test.sww',
765                                                 time_interval=[2.75, 2.75])
766        except AssertionError:
767            pass
768        else:
769            msg = 'Time interval should have raised an exception'
770            raise Exception, msg
771
772        # Cleanup
773        try:
774            os.remove(domain.get_name() + '.sww')
775        except:
776            pass
777            #FIXME(Ole): Windows won't allow removal of this
778
779
780
781
782
783    def test_another_runup_example(self):
784        """test_another_runup_example
785
786        Test runup example where actual timeseries at interpolated
787        points are tested.
788        """
789
790        #-----------------------------------------------------------------
791        # Setup computational domain
792        #-----------------------------------------------------------------
793        points, vertices, boundary = anuga.rectangular_cross(10, 10)# Basic mesh
794        domain = anuga.Domain(points, vertices, boundary) # Create domain
795        domain.set_default_order(2)
796        domain.set_quantities_to_be_stored(None)
797        domain.H0 = 1.0e-3
798
799        #-----------------------------------------------------------------
800        # Setup initial conditions
801        #-----------------------------------------------------------------
802        def topography(x, y):
803            return -x/2                              # linear bed slope
804
805        domain.set_quantity('elevation', topography)
806        domain.set_quantity('friction', 0.0)
807        domain.set_quantity('stage', expression='elevation')
808
809        #----------------------------------------------------------------
810        # Setup boundary conditions
811        #----------------------------------------------------------------
812        Br = anuga.Reflective_boundary(domain)           # Solid reflective wall
813        Bd = anuga.Dirichlet_boundary([-0.2, 0., 0.])    # Constant
814        domain.set_boundary({'left': Br, 'right': Bd, 'top': Br, 'bottom': Br})
815
816        #----------------------------------------------------------------
817        # Evolve system through time
818        #----------------------------------------------------------------
819        interpolation_points = [[0.4,0.5], [0.6,0.5], [0.8,0.5], [0.9,0.5]]
820        gauge_values = []
821        for _ in interpolation_points:
822            gauge_values.append([])
823
824        time = []
825        for t in domain.evolve(yieldstep=0.1, finaltime=5.0):
826            # Record time series at known points
827            time.append(domain.get_time())
828
829            stage = domain.get_quantity('stage')
830            w = stage.get_values(interpolation_points=interpolation_points)
831
832            for i, _ in enumerate(interpolation_points):
833                gauge_values[i].append(w[i])
834
835        #Reference (nautilus 26/6/2008)
836       
837        G0 = [-0.20000000000000001, -0.20000000000000001, -0.20000000000000001, -0.1958465301767274, -0.19059602372752493, -0.18448466250700923, -0.16979321333876071, -0.15976372740651074, -0.1575649333345176, -0.15710373731900584, -0.1530922283220747, -0.18836084336565725, -0.19921529311644628, -0.19923451799698919, -0.19923795414410964, -0.19923178806924047, -0.19925157557666154, -0.19930747801697429, -0.1993266428576112, -0.19932004930281799, -0.19929691326931867, -0.19926285267313795, -0.19916645449780995, -0.1988942593296438, -0.19900620256621993, -0.19914327423060865, -0.19918708440899577, -0.19921557252449132, -0.1992404368022069, -0.19925070370697717, -0.19925975477892274, -0.1992671090445659, -0.19927254203777162, -0.19927631910959256, -0.19927843552031504, -0.19927880339239365, -0.19927763204453783, -0.19927545249577633, -0.19927289590622824, -0.19927076261495152, -0.19926974334736983, -0.19927002562760332, -0.19927138236272213, -0.1992734501064522, -0.19927573251318192, -0.19927778936001547, -0.1992793776883893, -0.19928040577720926, -0.19928092586206753, -0.19928110982948721, -0.19928118887248453]
838
839        G1 = [-0.29999999999999993, -0.29999999999999993, -0.29139135018319512, -0.27257394456094503, -0.24141437432643265, -0.22089173942479151, -0.20796171092975532, -0.19874580192293825, -0.19014580508752857, -0.18421165368665365, -0.18020808282748838, -0.17518824759550247, -0.16436633464497749, -0.18714479115225544, -0.2045242886738807, -0.21011244240826329, -0.21151316017424124, -0.21048112933621732, -0.20772920477355789, -0.20489184334204144, -0.20286043930678221, -0.20094305756540246, -0.19948172752345467, -0.19886725178309209, -0.1986680808256765, -0.19860718133373548, -0.19862076543539733, -0.19888949069732539, -0.19932190310819023, -0.19982482967777454, -0.20036045468470615, -0.20064263130562704, -0.2007255389410077, -0.20068815669152493, -0.20057471332984647, -0.20042203940851802, -0.20026779013499779, -0.20015995671464712, -0.2000684005446525, -0.20001486753189174, -0.20000743467898013, -0.20003739771775905, -0.20008784600912621, -0.20013758305093884, -0.20017277554845025, -0.20018629233766885, -0.20018106288462198, -0.20016648079299326, -0.20015155958426531, -0.20014259747382254, -0.20014096648362509]
840       
841       
842        G2 = [-0.40000000000000002, -0.38885199453206343, -0.33425057028323962, -0.30154253721772117, -0.27624597383474103, -0.26037811196890087, -0.24415404585285019, -0.22941383121091052, -0.21613496492144549, -0.20418199946908885, -0.19506212965221825, -0.18851924999737435, -0.18271143344718843, -0.16910750701722474, -0.17963775224176018, -0.19442870510406052, -0.20164216917300118, -0.20467219452758528, -0.20608246104917802, -0.20640259931640445, -0.2054749739152594, -0.20380549124050265, -0.20227296931678532, -0.20095834856297176, -0.20000430919304379, -0.19946673053844086, -0.1990733499211611, -0.19882136174363013, -0.19877442300323914, -0.19905182154377868, -0.19943266521643804, -0.19988524183849191, -0.20018905307631765, -0.20031895675727809, -0.20033991149804931, -0.20031574232920274, -0.20027004750680638, -0.20020472427796293, -0.20013382447039607, -0.2000635018536408, -0.20001515436367642, -0.19998427691514989, -0.19997263083178107, -0.19998545383896535, -0.20000134502238734, -0.2000127243362736, -0.20001564474711939, -0.20001267360809977, -0.20002707579781318, -0.20004087961702843, -0.20004212947389177]
843       
844        G3 = [-0.45000000000000001, -0.38058172993544187, -0.33756059941741273, -0.31015371357441396, -0.29214769368562965, -0.27545447937118606, -0.25871585649808154, -0.24254276680581988, -0.22758633129006092, -0.21417276895743134, -0.20237184768790789, -0.19369491041576814, -0.18721625333717057, -0.1794243868465818, -0.17052113574042196, -0.18534300640363346, -0.19601184621026671, -0.20185028431829469, -0.20476187496918136, -0.20602933256960082, -0.20598569228739247, -0.20492643756666848, -0.20339695828762758, -0.20196440373022595, -0.20070304052919338, -0.19986227854052355, -0.19933020476408528, -0.19898034831018035, -0.19878317651286193, -0.19886879323961787, -0.19915860801206181, -0.19953675278099042, -0.19992828019602107, -0.20015957043092364, -0.20025268671087426, -0.20028559516444974, -0.20027084877341045, -0.20022991487243985, -0.20016234295579871, -0.20009131445092507, -0.20003149397006523, -0.19998473356473795, -0.19996011913447218, -0.19995647168667186, -0.19996526209120422, -0.19996600297827097, -0.19997268800221216, -0.19998682275066659, -0.20000372259781876, -0.20001628681983963, -0.2000173314086407]
845       
846        assert num.allclose(gauge_values[0], G0)
847        assert num.allclose(gauge_values[1], G1)
848        assert num.allclose(gauge_values[2], G2)
849        assert num.allclose(gauge_values[3], G3)
850
851    #####################################################
852
853
854    def test_initial_condition(self):
855        """test_initial_condition
856
857        Test that initial condition is output at time == 0 and that
858        computed values change as system evolves
859        """
860
861        from anuga.config import g
862        import copy
863
864        a = [0.0, 0.0]
865        b = [0.0, 2.0]
866        c = [2.0, 0.0]
867        d = [0.0, 4.0]
868        e = [2.0, 2.0]
869        f = [4.0, 0.0]
870
871        points = [a, b, c, d, e, f]
872        #             bac,     bce,     ecf,     dbe
873        vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
874
875        domain = anuga.Domain(points, vertices)
876
877        #Set up for a gradient of (3,0) at mid triangle (bce)
878        def slope(x, y):
879            return 3*x
880
881        h = 0.1
882        def stage(x, y):
883            return slope(x, y) + h
884
885        domain.set_quantity('elevation', slope)
886        domain.set_quantity('stage', stage)
887
888        # Allow slope limiters to work
889        # (FIXME (Ole): Shouldn't this be automatic in ANUGA?)
890        domain.distribute_to_vertices_and_edges()
891
892        initial_stage = copy.copy(domain.quantities['stage'].vertex_values)
893
894        domain.set_boundary({'exterior': anuga.Reflective_boundary(domain)})
895
896        domain.optimise_dry_cells = True
897
898        #Evolution
899        for t in domain.evolve(yieldstep=0.5, finaltime=2.0):
900            stage = domain.quantities['stage'].vertex_values
901
902            if t == 0.0:
903                assert num.allclose(stage, initial_stage)
904            else:
905                assert not num.allclose(stage, initial_stage)
906
907        os.remove(domain.get_name() + '.sww')
908
909    #####################################################
910
911    def test_second_order_flat_bed_onestep(self):
912
913        #Create basic mesh
914        points, vertices, boundary = anuga.rectangular(6, 6)
915
916        #Create shallow water domain
917        domain = anuga.Domain(points, vertices, boundary)
918        domain.set_default_order(2)
919
920        # Boundary conditions
921        Br = anuga.Reflective_boundary(domain)
922        Bd = anuga.Dirichlet_boundary([0.1, 0., 0.])
923        domain.set_boundary({'left': Bd, 'right': Br, 'top': Br, 'bottom': Br})
924
925        domain.check_integrity()
926
927        # Evolution
928        for t in domain.evolve(yieldstep=0.05, finaltime=0.05):
929            pass
930
931
932        # Data from earlier version of abstract_2d_finite_volumes
933        assert num.allclose(domain.recorded_min_timestep, 0.0396825396825) or \
934               num.allclose(domain.recorded_min_timestep, 0.0235282801879)
935               
936        assert num.allclose(domain.recorded_max_timestep, 0.0396825396825) or \
937               num.allclose(domain.recorded_max_timestep, 0.0235282801879)
938
939
940               
941        assert num.allclose(domain.quantities['stage'].centroid_values[:12],
942                            [0.00171396, 0.02656103, 0.00241523, 0.02656103,
943                             0.00241523, 0.02656103, 0.00241523, 0.02656103,
944                             0.00241523, 0.02656103, 0.00241523, 0.0272623], atol=1.0e-3) or \
945               num.allclose(domain.quantities['stage'].centroid_values[:12],
946                            [ 0.00053119,  0.02900893,  0.00077912,  0.02900893,
947                              0.00077912,  0.02900893,  0.00077912,  0.02900893,
948                              0.00077912,  0.02900893,  0.00077912,  0.02873746], atol=1.0e-3)
949
950        domain.distribute_to_vertices_and_edges()
951
952 
953
954        assert num.allclose(domain.quantities['stage'].vertex_values[:12,0],
955                            [ -1.96794125e-03,   2.65610347e-02,   0.00000000e+00,   2.65610347e-02,
956                              -8.67361738e-19,   2.65610347e-02,   4.33680869e-19,   2.65610347e-02,
957                              -2.16840434e-18,   2.65610347e-02,  -9.44042339e-05,   2.72623006e-02],
958                            atol =1.0e-3) or \
959                num.allclose(domain.quantities['stage'].vertex_values[:12,0],
960                            [ -5.51381419e-04,   5.74866732e-02,   1.00006808e-15,   5.72387383e-02,
961                              9.99851243e-16,   5.72387383e-02,   1.00050176e-15,   5.72387383e-02,
962                              9.99417563e-16,   5.72387383e-02,   1.09882029e-05,   5.66957956e-02],
963                             atol=1.0e-3)
964
965
966        assert num.allclose(domain.quantities['stage'].vertex_values[:12,1],
967                            [  5.14188587e-03,   2.65610347e-02,   0.00000000e+00,   2.65610347e-02,
968                               8.67361738e-19,   2.65610347e-02,  -4.33680869e-19,   2.65610347e-02,
969                               1.30104261e-18,   2.65610347e-02,   9.44042339e-05,   2.72623006e-02],
970                            atol =1.0e-3) or \
971               num.allclose(domain.quantities['stage'].vertex_values[:12,1],
972                           [  1.59356551e-03,   5.72387383e-02,   1.00006808e-15,   5.72387383e-02,
973                              1.00006808e-15,   5.72387383e-02,   9.99634403e-16,   5.72387383e-02,
974                              1.00050176e-15,   5.72387383e-02,  -1.09882029e-05,   1.47582915e-02],
975                            atol =1.0e-3)
976       
977        assert num.allclose(domain.quantities['stage'].vertex_values[:12,2],
978                            [ 0.00196794,  0.02656103,  0.00724568,  0.02656103,
979                              0.00724568,  0.02656103,  0.00724568,  0.02656103,
980                              0.00724568,  0.02656103,  0.00724568,  0.0272623 ], atol =1.0e-3) or \
981               num.allclose(domain.quantities['stage'].vertex_values[:12,2],
982                            [ 0.00055138, -0.02769862,  0.00233737, -0.02745068,
983                              0.00233737, -0.02745068,  0.00233737, -0.02745068,
984                              0.00233737, -0.02745068,  0.00233737,  0.01475829], atol =1.0e-3)
985
986
987        assert num.allclose(domain.quantities['xmomentum'].centroid_values[:12],
988                            [0.00113961, 0.01302432, 0.00148672,
989                             0.01302432, 0.00148672, 0.01302432,
990                             0.00148672, 0.01302432, 0.00148672 ,
991                             0.01302432, 0.00148672, 0.01337143], atol=1.0e-3) or \
992               num.allclose(domain.quantities['xmomentum'].centroid_values[:12],
993                        [ 0.00019529,  0.01425863,  0.00025665,
994                          0.01425863,  0.00025665,  0.01425863,
995                          0.00025665,  0.01425863,  0.00025665,
996                          0.01425863,  0.00025665,  0.014423  ], atol=1.0e-3)
997       
998        assert num.allclose(domain.quantities['ymomentum'].centroid_values[:12],
999                            [-2.91240050e-004 , 1.22721531e-004,
1000                             -1.22721531e-004,  1.22721531e-004 ,
1001                             -1.22721531e-004,  1.22721531e-004,
1002                             -1.22721531e-004 , 1.22721531e-004,
1003                             -1.22721531e-004,  1.22721531e-004,
1004                             -1.22721531e-004,  -4.57969873e-005], atol=1.0e-5) or \
1005               num.allclose(domain.quantities['ymomentum'].centroid_values[:12],
1006                            [ -6.38239364e-05,   2.16943067e-05,
1007                              -2.16943067e-05,   2.16943067e-05,
1008                              -2.16943067e-05,   2.16943067e-05,
1009                              -2.16943067e-05,   2.16943067e-05,
1010                              -2.16943067e-05,   2.16943067e-05,
1011                              -2.16943067e-05,  -4.62796434e-04], atol=1.0e-5)
1012       
1013        os.remove(domain.get_name() + '.sww')
1014
1015    def test_second_order_flat_bed_moresteps(self):
1016        from mesh_factory import rectangular
1017
1018        # Create basic mesh
1019        points, vertices, boundary = rectangular(6, 6)
1020
1021        # Create shallow water domain
1022        domain = anuga.Domain(points, vertices, boundary)
1023        domain.smooth = False
1024        domain.default_order = 2
1025
1026        # Boundary conditions
1027        Br = anuga.Reflective_boundary(domain)
1028        Bd = anuga.Dirichlet_boundary([0.1, 0., 0.])
1029        domain.set_boundary({'left': Bd, 'right': Br, 'top': Br, 'bottom': Br})
1030
1031        domain.check_integrity()
1032
1033        # Evolution
1034        for t in domain.evolve(yieldstep=0.05, finaltime=0.1):
1035            pass
1036
1037        # Data from earlier version of abstract_2d_finite_volumes
1038        #assert allclose(domain.recorded_min_timestep, 0.0396825396825)
1039        #assert allclose(domain.recorded_max_timestep, 0.0396825396825)
1040        #print domain.quantities['stage'].centroid_values
1041
1042        os.remove(domain.get_name() + '.sww')
1043
1044    def test_flatbed_first_order(self):
1045        from mesh_factory import rectangular
1046
1047        # Create basic mesh
1048        N = 8
1049        points, vertices, boundary = rectangular(N, N)
1050
1051        # Create shallow water domain
1052        domain = Domain(points, vertices, boundary)
1053        domain.smooth = False
1054        domain.default_order = 1
1055        domain.H0 = 1.0e-3 # As suggested in the manual
1056
1057        # Boundary conditions
1058        Br = anuga.Reflective_boundary(domain)
1059        Bd = anuga.Dirichlet_boundary([0.2, 0., 0.])
1060
1061        domain.set_boundary({'left': Bd, 'right': Br, 'top': Br, 'bottom': Br})
1062        domain.check_integrity()
1063
1064        # Evolution
1065        for t in domain.evolve(yieldstep=0.02, finaltime=0.5):
1066            pass
1067
1068        # FIXME: These numbers were from version before 25/10
1069        #assert allclose(domain.recorded_min_timestep, 0.0140413643926)
1070        #assert allclose(domain.recorded_max_timestep, 0.0140947355753)
1071
1072        for i in range(3):
1073            #assert allclose(domain.quantities['stage'].edge_values[:4,i],
1074            #                [0.10730244,0.12337617,0.11200126,0.12605666])
1075            assert num.allclose(domain.quantities['xmomentum'].\
1076                                    edge_values[:4,i],
1077                                [0.07610894,0.06901572,0.07284461,0.06819712])
1078            #assert allclose(domain.quantities['ymomentum'].edge_values[:4,i],
1079            #                [-0.0060238, -0.00157404, -0.00309633, -0.0001637])
1080
1081        os.remove(domain.get_name() + '.sww')
1082
1083    def test_flatbed_second_order(self):
1084        from mesh_factory import rectangular
1085
1086        # Create basic mesh
1087        N = 8
1088        points, vertices, boundary = rectangular(N, N)
1089
1090        # Create shallow water domain
1091        domain = Domain(points, vertices, boundary)
1092        domain.set_store_vertices_uniquely(True)
1093        domain.set_default_order(2)
1094
1095        # Boundary conditions
1096        Br = anuga.Reflective_boundary(domain)
1097        Bd = anuga.Dirichlet_boundary([0.2, 0., 0.])
1098
1099        domain.set_boundary({'left': Bd, 'right': Br, 'top': Br, 'bottom': Br})
1100        domain.check_integrity()
1101
1102        # Evolution
1103        for t in domain.evolve(yieldstep=0.01, finaltime=0.03):
1104            pass
1105
1106
1107       
1108        msg = 'min step was %f instead of %f' % (domain.recorded_min_timestep,
1109                                                 0.0155604907816)
1110
1111        assert num.allclose(domain.recorded_min_timestep, 0.0155604907816), msg
1112        assert num.allclose(domain.recorded_max_timestep, 0.0155604907816)
1113
1114
1115        assert num.allclose(domain.quantities['stage'].vertex_values[:4,0],
1116                            [-0.009, 0.0535,  0.0, 0.0535], atol=1.0e-3) or \
1117               num.allclose(domain.quantities['stage'].vertex_values[:4,0],
1118                      [-3.54158995e-03,1.22050959e-01,-2.36227400e-05,1.21501627e-01], atol=1.0e-3)
1119       
1120       
1121        assert num.allclose(domain.quantities['xmomentum'].vertex_values[:4,0],
1122                            [-0.008, 0.0368, 0.0, 0.0368], atol=1.0e-3) or \
1123               num.allclose(domain.quantities['xmomentum'].vertex_values[:4,0],
1124                      [-2.32056226e-03,9.10618822e-02, -1.06135035e-05,9.75175956e-02], atol=1.0e-3)
1125
1126        assert num.allclose(domain.quantities['ymomentum'].vertex_values[:4,0],
1127                            [ 0.002 , 6.0e-04, 0.0, 6.0e-04],
1128                            atol=1.0e-3) or \
1129               num.allclose(domain.quantities['ymomentum'].vertex_values[:4,0],
1130                            [ 1.43500775e-03,  6.07102924e-05,   1.59329371e-06,   8.44572599e-03],
1131                            atol=1.0e-3)
1132
1133        os.remove(domain.get_name() + '.sww')
1134
1135
1136    def test_flatbed_second_order_vmax_0(self):
1137        from mesh_factory import rectangular
1138
1139        # Create basic mesh
1140        N = 8
1141        points, vertices, boundary = rectangular(N, N)
1142
1143        # Create shallow water domain
1144        domain = Domain(points, vertices, boundary)
1145
1146        domain.set_store_vertices_uniquely(True)
1147        domain.set_default_order(2)
1148
1149
1150        # Boundary conditions
1151        Br = anuga.Reflective_boundary(domain)
1152        Bd = anuga.Dirichlet_boundary([0.2, 0., 0.])
1153
1154        domain.set_boundary({'left': Bd, 'right': Br, 'top': Br, 'bottom': Br})
1155        domain.check_integrity()
1156
1157        # Evolution
1158        for t in domain.evolve(yieldstep=0.01, finaltime=0.03):
1159            pass
1160
1161
1162        assert num.allclose(domain.recorded_min_timestep, 0.0210448446782) or \
1163               num.allclose(domain.recorded_min_timestep, 0.0155604907816)
1164               
1165        assert num.allclose(domain.recorded_max_timestep, 0.0210448446782) or \
1166               num.allclose(domain.recorded_min_timestep, 0.0155604907816)
1167
1168
1169        assert num.allclose(domain.quantities['xmomentum'].vertex_values[:4,0],
1170                            [ -2.32056226e-03,   9.10618822e-02,  -1.06135035e-05,   9.75175956e-02],
1171                            atol=1.0e-3)
1172
1173        assert num.allclose(domain.quantities['ymomentum'].vertex_values[:4,0],
1174                            [  1.43500775e-03,   6.07102924e-05,   1.59329371e-06,   8.44572599e-03],
1175                            atol=1.0e-3)
1176
1177        os.remove(domain.get_name() + '.sww')
1178
1179    def test_flatbed_second_order_distribute(self):
1180        #Use real data from anuga.abstract_2d_finite_volumes 2
1181        #painfully setup and extracted.
1182
1183        from mesh_factory import rectangular
1184
1185        # Create basic mesh
1186        N = 8
1187        points, vertices, boundary = rectangular(N, N)
1188
1189        # Create shallow water domain
1190        domain = Domain(points, vertices, boundary)
1191
1192        domain.set_store_vertices_uniquely(True)
1193        domain.set_default_order(2)
1194
1195        # Boundary conditions
1196        Br = anuga.Reflective_boundary(domain)
1197        Bd = anuga.Dirichlet_boundary([0.2, 0., 0.])
1198
1199        domain.set_boundary({'left': Bd, 'right': Br, 'top': Br, 'bottom': Br})
1200        domain.check_integrity()
1201
1202        for V in [False, True]:
1203            if V:
1204                # Set centroids as if system had been evolved
1205                L = num.zeros(2*N*N, num.float)
1206                L[:32] = [7.21205592e-003, 5.35214298e-002, 1.00910824e-002,
1207                          5.35439433e-002, 1.00910824e-002, 5.35439433e-002,
1208                          1.00910824e-002, 5.35439433e-002, 1.00910824e-002,
1209                          5.35439433e-002, 1.00910824e-002, 5.35439433e-002,
1210                          1.00910824e-002, 5.35393928e-002, 1.02344264e-002,
1211                          5.59605058e-002, 0.00000000e+000, 3.31027800e-004,
1212                          0.00000000e+000, 4.37962142e-005, 0.00000000e+000,
1213                          4.37962142e-005, 0.00000000e+000, 4.37962142e-005,
1214                          0.00000000e+000, 4.37962142e-005, 0.00000000e+000,
1215                          4.37962142e-005, 0.00000000e+000, 4.37962142e-005,
1216                          0.00000000e+000, 5.57305948e-005]
1217
1218                X = num.zeros(2*N*N, num.float)
1219                X[:32] = [6.48351607e-003, 3.68571894e-002, 8.50733285e-003,
1220                          3.68731327e-002, 8.50733285e-003, 3.68731327e-002,
1221                          8.50733285e-003, 3.68731327e-002, 8.50733285e-003,
1222                          3.68731327e-002, 8.50733285e-003, 3.68731327e-002,
1223                          8.50733285e-003, 3.68693861e-002, 8.65220973e-003,
1224                          3.85055387e-002, 0.00000000e+000, 2.86060840e-004,
1225                          0.00000000e+000, 3.58905503e-005, 0.00000000e+000,
1226                          3.58905503e-005, 0.00000000e+000, 3.58905503e-005,
1227                          0.00000000e+000, 3.58905503e-005, 0.00000000e+000,
1228                          3.58905503e-005, 0.00000000e+000, 3.58905503e-005,
1229                          0.00000000e+000, 4.57662812e-005]
1230
1231                Y = num.zeros(2*N*N, num.float)
1232                Y[:32] = [-1.39463104e-003, 6.15600298e-004, -6.03637382e-004,
1233                          6.18272251e-004, -6.03637382e-004, 6.18272251e-004,
1234                          -6.03637382e-004, 6.18272251e-004, -6.03637382e-004,
1235                          6.18272251e-004, -6.03637382e-004, 6.18272251e-004,
1236                          -6.03637382e-004, 6.18599320e-004, -6.74622797e-004,
1237                          -1.48934756e-004, 0.00000000e+000, -5.35079969e-005,
1238                          0.00000000e+000, -2.57264987e-005, 0.00000000e+000,
1239                          -2.57264987e-005, 0.00000000e+000, -2.57264987e-005,
1240                          0.00000000e+000, -2.57264987e-005, 0.00000000e+000,
1241                          -2.57264987e-005, 0.00000000e+000, -2.57264987e-005,
1242                          0.00000000e+000, -2.57635178e-005]
1243
1244                domain.set_quantity('stage', L, location='centroids')
1245                domain.set_quantity('xmomentum', X, location='centroids')
1246                domain.set_quantity('ymomentum', Y, location='centroids')
1247
1248                domain.check_integrity()
1249            else:
1250                # Evolution
1251                for t in domain.evolve(yieldstep=0.01, finaltime=0.03):
1252                    pass
1253
1254               
1255                assert num.allclose(domain.recorded_min_timestep, 0.0155604907816)
1256                assert num.allclose(domain.recorded_max_timestep, 0.0155604907816)
1257
1258            #print domain.quantities['stage'].centroid_values[:4]
1259            #print domain.quantities['xmomentum'].centroid_values[:4]
1260            #print domain.quantities['ymomentum'].centroid_values[:4]                       
1261               
1262            #Centroids were correct but not vertices.
1263            #Hence the check of distribute below.
1264
1265            if not V:
1266
1267                assert num.allclose(domain.quantities['stage'].centroid_values[:4],
1268                               [0.00725574, 0.05350737, 0.01008413, 0.0535293], atol=1.0e-3) or \
1269                       num.allclose(domain.quantities['stage'].centroid_values[:4],
1270                               [0.00318259,  0.06261678,  0.00420215,  0.06285189], atol=1.0e-3)
1271               
1272                assert num.allclose(domain.quantities['xmomentum'].centroid_values[:4],
1273                               [0.00654964, 0.03684904, 0.00852561, 0.03686323],atol=1.0e-3) or \
1274                       num.allclose(domain.quantities['xmomentum'].centroid_values[:4],
1275                               [0.00218173, 0.04482164, 0.0026334,  0.04491656],atol=1.0e-3)       
1276
1277                assert num.allclose(domain.quantities['ymomentum'].centroid_values[:4],
1278                               [-0.00143169, 0.00061027, -0.00062325, 0.00061408],atol=1.0e-3) or \
1279                       num.allclose(domain.quantities['ymomentum'].centroid_values[:4],
1280                [-6.46340592e-04,-6.16702557e-05,-2.83424134e-04, 6.48556590e-05],atol=1.0e-3)
1281
1282                                   
1283                assert num.allclose(domain.quantities['xmomentum'].centroid_values[17], 0.0,
1284                                    atol=3.0e-4)               
1285            else:
1286                assert num.allclose(domain.quantities['xmomentum'].\
1287                                        centroid_values[17],
1288                                    0.00028606084)
1289                return #FIXME - Bailout for V True
1290
1291            import copy
1292
1293            XX = copy.copy(domain.quantities['xmomentum'].centroid_values)
1294            assert num.allclose(domain.quantities['xmomentum'].centroid_values,
1295                                XX)
1296
1297            domain.distribute_to_vertices_and_edges()
1298
1299            assert num.allclose(domain.quantities['xmomentum'].centroid_values[17], 0.0, atol=3.0e-4)
1300
1301            assert num.allclose(domain.quantities['ymomentum'].vertex_values[:4,0],
1302                                [ 1.84104149e-03, 6.05658846e-04, 1.77092716e-07, 6.10687334e-04],
1303                                atol=1.0e-4) or \
1304                   num.allclose(domain.quantities['ymomentum'].vertex_values[:4,0],
1305                                [1.43500775e-03, 6.07102924e-05, 1.59329371e-06, 8.44572599e-03],
1306                                atol=1.0e-4)             
1307
1308            assert num.allclose(domain.quantities['xmomentum'].vertex_values[:4,0],
1309                                [ -8.31184293e-03, 3.68841505e-02, -2.42843889e-06, 3.68900189e-02],
1310                                atol=1.0e-4) or \
1311                   num.allclose(domain.quantities['xmomentum'].vertex_values[:4,0],
1312                                [-2.32056226e-03, 9.10618822e-02, -1.06135035e-05, 9.75175956e-02],
1313                                rtol=1.0e-2)             
1314
1315
1316        os.remove(domain.get_name() + '.sww')
1317
1318
1319    def test_bedslope_problem_second_order_more_steps(self):
1320        """test_bedslope_problem_second_order_more_step
1321
1322        Test shallow water balanced finite volumes
1323        """
1324
1325        from mesh_factory import rectangular
1326
1327        # Create basic mesh
1328        points, vertices, boundary = rectangular(6, 6)
1329
1330        # Create shallow water domain
1331        domain = Domain(points, vertices, boundary)
1332
1333        domain.set_store_vertices_uniquely(True)
1334        domain.set_default_order(2)
1335
1336        # Bed-slope and friction at vertices (and interpolated elsewhere)
1337        def x_slope(x, y):
1338            return -x/3
1339
1340        domain.set_quantity('elevation', x_slope)
1341
1342        # Boundary conditions
1343        Br = anuga.Reflective_boundary(domain)
1344        domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
1345
1346        # Initial condition
1347        domain.set_quantity('stage', expression='elevation+0.05')
1348        domain.check_integrity()
1349
1350        assert num.allclose(domain.quantities['stage'].centroid_values,
1351                            [ 0.01296296,  0.03148148,  0.01296296,
1352                              0.03148148,  0.01296296,  0.03148148,
1353                              0.01296296,  0.03148148,  0.01296296,
1354                              0.03148148,  0.01296296,  0.03148148,
1355                             -0.04259259, -0.02407407, -0.04259259,
1356                             -0.02407407, -0.04259259, -0.02407407,
1357                             -0.04259259, -0.02407407, -0.04259259,
1358                             -0.02407407, -0.04259259, -0.02407407,
1359                             -0.09814815, -0.07962963, -0.09814815,
1360                             -0.07962963, -0.09814815, -0.07962963,
1361                             -0.09814815, -0.07962963, -0.09814815,
1362                             -0.07962963, -0.09814815, -0.07962963,
1363                             -0.1537037 , -0.13518519, -0.1537037,
1364                             -0.13518519, -0.1537037,  -0.13518519,
1365                             -0.1537037 , -0.13518519, -0.1537037,
1366                             -0.13518519, -0.1537037,  -0.13518519,
1367                             -0.20925926, -0.19074074, -0.20925926,
1368                             -0.19074074, -0.20925926, -0.19074074,
1369                             -0.20925926, -0.19074074, -0.20925926,
1370                             -0.19074074, -0.20925926, -0.19074074,
1371                             -0.26481481, -0.2462963,  -0.26481481,
1372                             -0.2462963,  -0.26481481, -0.2462963,
1373                             -0.26481481, -0.2462963,  -0.26481481,
1374                             -0.2462963,  -0.26481481, -0.2462963])
1375
1376        # Evolution
1377        for t in domain.evolve(yieldstep = 0.05, finaltime = 0.5):
1378            pass
1379
1380
1381        assert num.allclose(domain.quantities['stage'].centroid_values,
1382     [-0.02901283, -0.01619385, -0.03040423, -0.01564474, -0.02936756, -0.01507953,
1383      -0.02858108, -0.01491531, -0.02793549, -0.0147037,  -0.02792804, -0.014363,
1384      -0.07794301, -0.05951952, -0.07675098, -0.06182336, -0.07736607, -0.06079504,
1385      -0.07696764, -0.06039043, -0.07708793, -0.0601453,  -0.07669911, -0.06020719,
1386      -0.12223185, -0.10857309, -0.12286676, -0.10837591, -0.12386938, -0.10842744,
1387      -0.12363769, -0.10790002, -0.12304837, -0.10871278, -0.12543768, -0.10961026,
1388      -0.15664473, -0.14630207, -0.15838364, -0.14910271, -0.15804002, -0.15029627,
1389      -0.15829717, -0.1503869,  -0.15852604, -0.14971109, -0.15856346, -0.15205092,
1390      -0.20900931, -0.19658843, -0.20669607, -0.19558708, -0.20654186, -0.19492423,
1391      -0.20438765, -0.19492931, -0.20644142, -0.19423147, -0.20237449, -0.19198454,
1392      -0.13699658, -0.14209126, -0.13600697, -0.14334968, -0.1347657,  -0.14224247,
1393      -0.13442376, -0.14136926, -0.13501004, -0.14339389, -0.13479263, -0.14304073], atol=1.0e-2) or \
1394              num.allclose(domain.quantities['stage'].centroid_values,     
1395     [-0.03393968, -0.0166423,  -0.03253538, -0.01722023, -0.03270405, -0.01728606,
1396      -0.03277786, -0.0173903,  -0.03333736, -0.01743236, -0.03189526, -0.01463918,
1397      -0.07951756, -0.06410763, -0.07847973, -0.06350794, -0.07842429, -0.06240852,
1398      -0.07808697, -0.06255924, -0.07854662, -0.06322442, -0.07867314, -0.06287121,
1399      -0.11533356, -0.10559238, -0.11971301, -0.10742123, -0.1215759 , -0.10830046,
1400      -0.12202867, -0.10831703, -0.122214,   -0.10854099, -0.12343779, -0.11035803,
1401      -0.15725714, -0.14300757, -0.15559898, -0.1447275 , -0.15478568, -0.14483551,
1402      -0.15461918, -0.14489704, -0.15462074, -0.14516256, -0.15522298, -0.1452902,
1403      -0.22637615, -0.19192974, -0.20922654, -0.1907441 , -0.20900039, -0.19074809,
1404      -0.20897969, -0.19073365, -0.209195,   -0.19071396, -0.20922513, -0.19067714,
1405      -0.11357515, -0.14185801, -0.13224763, -0.14395805, -0.13379438, -0.14497114,
1406      -0.13437773, -0.14536013, -0.13607796, -0.14799629, -0.13148351, -0.15568502], atol=1.0e-1)
1407
1408
1409
1410        assert num.allclose(domain.quantities['xmomentum'].centroid_values,
1411               [ 0.00478273,  0.003297,    0.00471129,  0.00320957,  0.00462171,  0.00320135,
1412                 0.00458295,  0.00317193,  0.00451704,  0.00314308,  0.00442684,  0.00320466,
1413                 0.01512907,  0.01150756,  0.01604672,  0.01156605,  0.01583911,  0.01135809,
1414                 0.01578499,  0.01132479,  0.01543668,  0.01100614,  0.01570445,  0.0120152,
1415                 0.04019477,  0.02721469,  0.03509982,  0.02735229,  0.03369315,  0.02727871,
1416                 0.03317931,  0.02706421,  0.03332704,  0.02722779,  0.03170258,  0.02556134,
1417                 0.07157025,  0.06074271,  0.07249738,  0.05570979,  0.07311261,  0.05428175,
1418                 0.07316986,  0.05379702,  0.0719581,   0.05230996,  0.07034837,  0.05468702,
1419                 0.08145001,  0.07753479,  0.08148804,  0.08119069,  0.08247295,  0.08134969,
1420                 0.0823216,   0.081411,    0.08190964,  0.08151441,  0.08163076,  0.08166174,
1421                 0.03680205,  0.0768216,   0.03943625,  0.07791183,  0.03930529,  0.07760588,
1422                 0.03949756,  0.07839929,  0.03992892,  0.08001416,  0.04444335,  0.08628738],
1423                            atol=1.0e-2) or \
1424               num.allclose(domain.quantities['xmomentum'].centroid_values,
1425               [ 0.00178414,  0.00147791,  0.00373636,  0.00169124,  0.00395649,  0.0014468,
1426                 0.00387617,  0.00135572,  0.00338418,  0.00134554,  0.00404961,  0.00252769,
1427                 0.01365204,  0.00890416,  0.01381613,  0.00986246,  0.01419385,  0.01145017,
1428                 0.01465116,  0.01125933,  0.01407359,  0.01055426,  0.01403563,  0.01095544,
1429                 0.04653827,  0.03018236,  0.03709973,  0.0265533 ,  0.0337694 ,  0.02541724,
1430                 0.03304266,  0.02535335,  0.03264548,  0.02484769,  0.03047682,  0.02205757,
1431                 0.07400338,  0.06470583,  0.07756503,  0.06098108,  0.07942593,  0.06086531,
1432                 0.07977427,  0.06074404,  0.07979513,  0.06019911,  0.07806395,  0.06011152,
1433                 0.07305045,  0.07883894,  0.08120393,  0.08166623,  0.08180501,  0.08166251,
1434                 0.0818353 ,  0.08169641,  0.08173762,  0.08174118,  0.08176467,  0.08181817,
1435                 0.01549926,  0.08259719,  0.01835423,  0.07302656,  0.01672924,  0.07198839,
1436                 0.01676006,  0.07223233,  0.01775672,  0.07362164,  0.01955846,  0.09361223],
1437                            atol=1.0e-2)
1438
1439
1440        assert num.allclose(domain.quantities['ymomentum'].centroid_values,
1441                            [ -1.09771684e-05,  -2.60328801e-05,  -1.03481959e-05,  -7.75907380e-05,
1442                              -5.00409090e-05,  -7.83807512e-05,  -3.60509918e-05,  -6.19321031e-05,
1443                              -1.40041903e-05,  -2.95707259e-05,   3.90296618e-06,   1.87556544e-05,
1444                              9.27848053e-05,   6.66937557e-07,   1.00653468e-04,   8.24734209e-06,
1445                              -1.04548672e-05,  -4.40402988e-05,  -2.95549946e-05,  -1.86360736e-05,
1446                              1.12527016e-04,   1.27240681e-04,   2.02147284e-04,   9.18457482e-05,
1447                              1.41781748e-03,   7.23407624e-04,   5.09160779e-04,   1.29136939e-04,
1448                              -4.70131286e-05,  -1.00180290e-04,  -1.76806614e-05,  -4.19421384e-06,
1449                              -6.17759681e-05,  -3.02124967e-05,   4.32689360e-04,   5.49717934e-04,
1450                              1.15031101e-03,   1.02737170e-03,   5.77937840e-04,   3.36230967e-04,
1451                              5.44877516e-04,  -7.28594977e-05,   4.60064858e-04,  -3.94125434e-05,
1452                              7.48242964e-04,   2.88528341e-04,   6.25148041e-05,  -1.74477175e-04,
1453                              -5.06603166e-05,   7.07720999e-04,  -2.04937748e-04,   3.38595573e-05,
1454                              -4.64116229e-05,   1.49325340e-04,  -2.41342281e-05,   1.83817970e-04,
1455                              -1.44417277e-05,   2.47823834e-04,   7.91185571e-05,   1.71615793e-04,
1456                              1.56883043e-03,   8.39352974e-04,   3.23353846e-03,   1.70597880e-03,
1457                              2.27789107e-03,   1.48928169e-03,   2.09854126e-03,   1.50248643e-03,
1458                              2.83029467e-03,   1.09151499e-03,   6.52455118e-03,  -2.04468968e-03],
1459                            atol=1.0e-3) or \
1460               num.allclose(domain.quantities['ymomentum'].centroid_values,
1461                             [ -1.24810991e-04,  -3.08228767e-04,  -1.56701128e-04,  -1.01904208e-04,
1462                               -3.36282053e-05,  -1.17956840e-04,  -3.55986664e-05,  -9.38578996e-05,
1463                               7.13704069e-05,   2.47022380e-05,   1.71121489e-04,   2.65941677e-04,
1464                               6.90055205e-04,   1.99195585e-04,   1.33804448e-04,  -1.66563316e-04,
1465                               -2.00962830e-04,  -3.81664130e-05,  -9.50456053e-05,  -3.14620186e-06,
1466                               1.29388102e-04,   3.16945980e-04,   4.77556581e-04,   2.57217342e-04,
1467                               1.42300612e-03,   9.60776359e-04,   5.08941026e-04,   1.06939990e-04,
1468                               6.37673950e-05,  -2.69783047e-04,  -8.55760509e-05,  -2.12987309e-04,
1469                               -5.86840949e-06,  -9.75751293e-05,   8.25447727e-04,   1.14139065e-03,
1470                               8.56206468e-04,   3.83113329e-04,   1.75041847e-04,   4.39999200e-04,
1471                               3.75156469e-04,   2.48774698e-04,   4.09671654e-04,   2.07125615e-04,
1472                               4.59587647e-04,   2.70581830e-04,  -1.24082302e-06,  -4.29155678e-04,
1473                               -9.66841218e-03,   4.93278794e-04,  -5.25778806e-06,  -4.90396857e-05,
1474                               -9.75373988e-06,   7.28023591e-06,  -5.20499868e-06,   3.61013683e-05,
1475                               -7.54919544e-06,   4.14115771e-05,  -1.35778834e-05,  -2.23991903e-05,
1476                               3.63635844e-02,   5.29865244e-04,   5.13015379e-03,   1.19233296e-03,
1477                               4.70681275e-04,   2.62292296e-04,  -1.28084045e-04,   7.04826916e-04,
1478                               1.50377987e-04,   1.35053814e-03,   1.30710492e-02,   1.93011958e-03],
1479                            atol=1.0e-1)
1480                           
1481
1482
1483        os.remove(domain.get_name() + '.sww')
1484
1485
1486    def test_temp_play(self):
1487        from mesh_factory import rectangular
1488
1489        # Create basic mesh
1490        points, vertices, boundary = rectangular(5, 5)
1491
1492        # Create shallow water domain
1493        domain = Domain(points, vertices, boundary)
1494
1495        domain.set_store_vertices_uniquely(True)
1496        domain.set_default_order(2)
1497
1498
1499
1500        # Bed-slope and friction at vertices (and interpolated elsewhere)
1501        def x_slope(x, y):
1502            return -x/3
1503
1504        domain.set_quantity('elevation', x_slope)
1505
1506        # Boundary conditions
1507        Br = anuga.Reflective_boundary(domain)
1508        domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
1509
1510        # Initial condition
1511        domain.set_quantity('stage', expression='elevation+0.05')
1512        domain.check_integrity()
1513
1514        # Evolution
1515        for t in domain.evolve(yieldstep=0.05, finaltime=0.1):
1516            pass
1517
1518
1519        assert num.allclose(domain.quantities['stage'].centroid_values[:4],
1520                            [ 0.01,   0.015,  0.01,  0.015], atol=1.0e-2)
1521                           
1522        assert num.allclose(domain.quantities['xmomentum'].centroid_values[:4],
1523                            [ 0.015,  0.01,  0.015,  0.01], atol=1.0e-2)
1524       
1525        assert num.allclose(domain.quantities['ymomentum'].centroid_values[:4],
1526                            [  0.0,  0.0,  0.0,   0.0]
1527                            , atol=1.0e-3)
1528
1529        os.remove(domain.get_name() + '.sww')
1530
1531    def test_complex_bed(self):
1532        # No friction is tested here
1533
1534        from mesh_factory import rectangular
1535
1536        N = 12
1537        points, vertices, boundary = rectangular(N, N/2, len1=1.2, len2=0.6,
1538                                                 origin=(-0.07, 0))
1539
1540
1541        domain = Domain(points, vertices, boundary)
1542        domain.smooth = False
1543        domain.set_default_order(2)
1544        domain.set_timestepping_method('rk2')
1545        domain.set_beta(1.0)
1546
1547        inflow_stage = 0.1
1548        Z = Weir(inflow_stage)
1549        domain.set_quantity('elevation', Z)
1550
1551        Br = anuga.Reflective_boundary(domain)
1552        Bd = anuga.Dirichlet_boundary([inflow_stage, 0.0, 0.0])
1553        domain.set_boundary({'left': Bd, 'right': Br, 'bottom': Br, 'top': Br})
1554
1555        domain.set_quantity('stage', expression='elevation')
1556
1557        for t in domain.evolve(yieldstep=0.02, finaltime=0.2):
1558            pass
1559
1560        #FIXME: These numbers were from version before 25/10
1561        #assert allclose(domain.quantities['stage'].centroid_values,
1562# [3.95822638e-002,  5.61022588e-002,  4.66437868e-002,  5.73081011e-002,
1563#  4.72394613e-002,  5.74684939e-002,  4.74309483e-002,  5.77458084e-002,
1564#  4.80628177e-002,  5.85656225e-002,  4.90498542e-002,  6.02609831e-002,
1565#  1.18470315e-002,  1.75136443e-002,  1.18035266e-002,  2.15565695e-002,
1566#  1.31620268e-002,  2.14351640e-002,  1.32351076e-002,  2.15450687e-002,
1567#  1.36414028e-002,  2.24274619e-002,  1.51689511e-002,  2.21789655e-002,
1568# -7.54337535e-003, -6.86362021e-004, -7.74146760e-003, -1.83756530e-003,
1569# -8.16773628e-003, -4.49916813e-004, -8.08202599e-003, -3.91118720e-004,
1570# -8.10292716e-003, -3.88584984e-004, -7.35226124e-003,  2.73985295e-004,
1571#  1.86166683e-001,  8.74070369e-002,  1.86166712e-001,  8.74035875e-002,
1572#  6.11666935e-002, -3.76173225e-002, -6.38333276e-002, -3.76147365e-002,
1573#  6.11666725e-002,  8.73846774e-002,  1.86166697e-001,  8.74171550e-002,
1574# -4.83333333e-002,  1.18333333e-001, -4.83333333e-002,  1.18333333e-001,
1575# -4.83333333e-002, -6.66666667e-003, -1.73333333e-001, -1.31666667e-001,
1576# -1.73333333e-001, -6.66666667e-003, -4.83333333e-002,  1.18333333e-001,
1577# -2.48333333e-001, -2.31666667e-001, -2.48333333e-001, -2.31666667e-001,
1578# -2.48333333e-001, -2.31666667e-001, -2.48333333e-001, -2.31666667e-001,
1579# -2.48333333e-001, -2.31666667e-001, -2.48333333e-001, -2.31666667e-001,
1580# -4.65000000e-001, -3.65000000e-001, -4.65000000e-001, -3.65000000e-001,
1581# -4.65000000e-001, -3.65000000e-001, -4.65000000e-001, -3.65000000e-001,
1582# -4.65000000e-001, -3.65000000e-001, -4.65000000e-001, -3.65000000e-001,
1583# -5.98333333e-001, -5.81666667e-001, -5.98333333e-001, -5.81666667e-001,
1584# -5.98333333e-001, -5.81666667e-001, -5.98333333e-001, -5.81666667e-001,
1585# -5.98333333e-001, -5.81666667e-001, -5.98333333e-001, -5.81666667e-001,
1586# -6.48333333e-001, -6.31666667e-001, -6.48333333e-001, -6.31666667e-001,
1587# -6.48333333e-001, -6.31666667e-001, -6.48333333e-001, -6.31666667e-001,
1588# -6.48333333e-001, -6.31666667e-001, -6.48333333e-001, -6.31666667e-001,
1589# -5.31666667e-001, -5.98333333e-001, -5.31666667e-001, -5.98333333e-001,
1590# -5.31666667e-001, -5.98333333e-001, -5.31666667e-001, -5.98333333e-001,
1591# -5.31666667e-001, -5.98333333e-001, -5.31666667e-001, -5.98333333e-001,
1592# -4.98333333e-001, -4.81666667e-001, -4.98333333e-001, -4.81666667e-001,
1593# -4.98333333e-001, -4.81666667e-001, -4.98333333e-001, -4.81666667e-001,
1594# -4.98333333e-001, -4.81666667e-001, -4.98333333e-001, -4.81666667e-001,
1595# -5.48333333e-001, -5.31666667e-001, -5.48333333e-001, -5.31666667e-001,
1596# -5.48333333e-001, -5.31666667e-001, -5.48333333e-001, -5.31666667e-001,
1597# -5.48333333e-001, -5.31666667e-001, -5.48333333e-001, -5.31666667e-001])
1598
1599        os.remove(domain.get_name() + '.sww')
1600
1601
1602    def test_tight_slope_limiters(self):
1603        """Test that new slope limiters (Feb 2007) don't induce extremely
1604        small timesteps. This test actually reveals the problem as it
1605        was in March-April 2007
1606        """
1607        import time, os
1608        from Scientific.IO.NetCDF import NetCDFFile
1609        from mesh_factory import rectangular_cross
1610
1611        # Create basic mesh
1612        points, vertices, boundary = rectangular_cross(2, 2)
1613
1614        # Create shallow water domain
1615        domain = Domain(points, vertices, boundary)
1616        domain.set_default_order(2)
1617        domain.set_beta(1.0)
1618        domain.set_timestepping_method('euler')
1619        #domain.set_CFL(0.5)
1620       
1621
1622        # This will pass
1623        #domain.tight_slope_limiters = 1
1624        #domain.H0 = 0.01
1625
1626        # This will fail
1627        #domain.tight_slope_limiters = 1
1628        #domain.H0 = 0.001
1629
1630        # This will pass provided C extension implements limiting of
1631        # momentum in _compute_speeds
1632        #domain.tight_slope_limiters = 1
1633        #domain.H0 = 0.001
1634        #domain.protect_against_isolated_degenerate_timesteps = True
1635
1636        # Set some field values
1637        domain.set_quantity('elevation', lambda x,y: -x)
1638        domain.set_quantity('friction', 0.03)
1639
1640        # Boundary conditions
1641        B = Transmissive_boundary(domain)
1642        domain.set_boundary({'left': B, 'right': B, 'top': B, 'bottom': B})
1643
1644        # Initial condition - with jumps
1645        bed = domain.quantities['elevation'].vertex_values
1646        stage = num.zeros(bed.shape, num.float)
1647
1648        h = 0.3
1649        for i in range(stage.shape[0]):
1650            if i % 2 == 1:
1651                stage[i,:] = bed[i,:] + h
1652            else:
1653                stage[i,:] = bed[i,:]
1654
1655        domain.set_quantity('stage', stage)
1656
1657        domain.distribute_to_vertices_and_edges()
1658
1659        domain.set_name('tight_limiters')
1660        domain.smooth = True
1661        domain.reduction = mean
1662        domain.set_datadir('.')
1663        domain.smooth = False
1664        domain.store = True
1665
1666        # Evolution
1667        for t in domain.evolve(yieldstep=0.1, finaltime=0.3):
1668            #domain.write_time(track_speeds=True)
1669            stage = domain.quantities['stage'].vertex_values
1670
1671            # Get NetCDF
1672            #fid = NetCDFFile(domain.writer.filename, netcdf_mode_r)
1673            #stage_file = fid.variables['stage']
1674
1675            #fid.close()
1676
1677        os.remove(domain.writer.filename)
1678
1679
1680
1681    def test_pmesh2Domain(self):
1682         import os
1683         import tempfile
1684
1685         fileName = tempfile.mktemp(".tsh")
1686         file = open(fileName, "w")
1687         file.write("4 3 # <vertex #> <x> <y> [attributes]\n \
16880 0.0 0.0 0.0 0.0 0.01 \n \
16891 1.0 0.0 10.0 10.0 0.02  \n \
16902 0.0 1.0 0.0 10.0 0.03  \n \
16913 0.5 0.25 8.0 12.0 0.04  \n \
1692# Vert att title  \n \
1693elevation  \n \
1694stage  \n \
1695friction  \n \
16962 # <triangle #> [<vertex #>] [<neigbouring triangle #>]  \n\
16970 0 3 2 -1  -1  1 dsg\n\
16981 0 1 3 -1  0 -1   ole nielsen\n\
16994 # <segment #> <vertex #>  <vertex #> [boundary tag] \n\
17000 1 0 2 \n\
17011 0 2 3 \n\
17022 2 3 \n\
17033 3 1 1 \n\
17043 0 # <x> <y> [attributes] ...Mesh Vertices... \n \
17050 216.0 -86.0 \n \
17061 160.0 -167.0 \n \
17072 114.0 -91.0 \n \
17083 # <vertex #>  <vertex #> [boundary tag] ...Mesh Segments... \n \
17090 0 1 0 \n \
17101 1 2 0 \n \
17112 2 0 0 \n \
17120 # <x> <y> ...Mesh Holes... \n \
17130 # <x> <y> <attribute>...Mesh Regions... \n \
17140 # <x> <y> <attribute>...Mesh Regions, area... \n\
1715#Geo reference \n \
171656 \n \
1717140 \n \
1718120 \n")
1719         file.close()
1720
1721         tags = {}
1722         b1 =  anuga.Dirichlet_boundary(conserved_quantities = num.array([0.0]))
1723         b2 =  anuga.Dirichlet_boundary(conserved_quantities = num.array([1.0]))
1724         b3 =  anuga.Dirichlet_boundary(conserved_quantities = num.array([2.0]))
1725         tags["1"] = b1
1726         tags["2"] = b2
1727         tags["3"] = b3
1728
1729         domain = Domain(mesh_filename=fileName)
1730                         # verbose=True, use_cache=True)
1731
1732         ## check the quantities
1733         answer = [[0., 8., 0.],
1734                   [0., 10., 8.]]
1735         assert num.allclose(domain.quantities['elevation'].vertex_values,
1736                             answer)
1737
1738         answer = [[0., 12., 10.],
1739                   [0., 10., 12.]]
1740         assert num.allclose(domain.quantities['stage'].vertex_values,
1741                             answer)
1742
1743         answer = [[0.01, 0.04, 0.03],
1744                   [0.01, 0.02, 0.04]]
1745         assert num.allclose(domain.quantities['friction'].vertex_values,
1746                             answer)
1747
1748         tagged_elements = domain.get_tagged_elements()
1749         assert num.allclose(tagged_elements['dsg'][0], 0)
1750         assert num.allclose(tagged_elements['ole nielsen'][0], 1)
1751
1752         msg = "test_tags_to_boundaries failed. Single boundary wasn't added."
1753         self.failUnless( domain.boundary[(1, 0)]  == '1', msg)
1754         self.failUnless( domain.boundary[(1, 2)]  == '2', msg)
1755         self.failUnless( domain.boundary[(0, 1)]  == '3', msg)
1756         self.failUnless( domain.boundary[(0, 0)]  == 'exterior', msg)
1757         msg = "test_pmesh2Domain Too many boundaries"
1758         self.failUnless( len(domain.boundary)  == 4, msg)
1759
1760         # FIXME change to use get_xllcorner
1761         msg = 'Bad geo-reference'
1762         self.failUnless(domain.geo_reference.xllcorner  == 140.0, msg)
1763
1764         domain = Domain(fileName)
1765
1766         answer = [[0., 8., 0.],
1767                   [0., 10., 8.]]
1768         assert num.allclose(domain.quantities['elevation'].vertex_values,
1769                             answer)
1770
1771         answer = [[0., 12., 10.],
1772                   [0., 10., 12.]]
1773         assert num.allclose(domain.quantities['stage'].vertex_values,
1774                             answer)
1775
1776         answer = [[0.01, 0.04, 0.03],
1777                   [0.01, 0.02, 0.04]]
1778         assert num.allclose(domain.quantities['friction'].vertex_values,
1779                             answer)
1780
1781         tagged_elements = domain.get_tagged_elements()
1782         assert num.allclose(tagged_elements['dsg'][0], 0)
1783         assert num.allclose(tagged_elements['ole nielsen'][0], 1)
1784
1785         msg = "test_tags_to_boundaries failed. Single boundary wasn't added."
1786         self.failUnless(domain.boundary[(1, 0)]  == '1', msg)
1787         self.failUnless(domain.boundary[(1, 2)]  == '2', msg)
1788         self.failUnless(domain.boundary[(0, 1)]  == '3', msg)
1789         self.failUnless(domain.boundary[(0, 0)]  == 'exterior', msg)
1790         msg = "test_pmesh2Domain Too many boundaries"
1791         self.failUnless(len(domain.boundary)  == 4, msg)
1792
1793         # FIXME change to use get_xllcorner
1794         msg = 'Bad geo_reference'
1795         self.failUnless(domain.geo_reference.xllcorner  == 140.0, msg)
1796
1797         os.remove(fileName)
1798
1799    def test_get_lone_vertices(self):
1800        a = [0.0, 0.0]
1801        b = [0.0, 2.0]
1802        c = [2.0, 0.0]
1803        d = [0.0, 4.0]
1804        e = [2.0, 2.0]
1805        f = [4.0, 0.0]
1806
1807        points = [a, b, c, d, e, f]
1808        #             bac,     bce,     ecf,     dbe
1809        vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4] ]
1810        boundary = {(0, 0): 'Third',
1811                    (0, 2): 'First',
1812                    (2, 0): 'Second',
1813                    (2, 1): 'Second',
1814                    (3, 1): 'Second',
1815                    (3, 2): 'Third'}
1816
1817        domain = Domain(points, vertices, boundary)
1818        domain.get_lone_vertices()
1819
1820    def test_fitting_using_shallow_water_domain(self):
1821        #Mesh in zone 56 (absolute coords)
1822
1823        x0 = 314036.58727982
1824        y0 = 6224951.2960092
1825
1826        a = [x0+0.0, y0+0.0]
1827        b = [x0+0.0, y0+2.0]
1828        c = [x0+2.0, y0+0.0]
1829        d = [x0+0.0, y0+4.0]
1830        e = [x0+2.0, y0+2.0]
1831        f = [x0+4.0, y0+0.0]
1832
1833        points = [a, b, c, d, e, f]
1834
1835        #             bac,     bce,     ecf,     dbe
1836        elements = [[1,0,2], [1,2,4], [4,2,5], [3,1,4] ]
1837
1838        # absolute going in ..
1839        mesh4 = Domain(points, elements, geo_reference=Geo_reference(56, 0, 0))
1840        mesh4.check_integrity()
1841        quantity = Quantity(mesh4)
1842
1843        # Get (enough) datapoints (relative to georef)
1844        data_points_rel = [[ 0.66666667, 0.66666667],
1845                           [ 1.33333333, 1.33333333],
1846                           [ 2.66666667, 0.66666667],
1847                           [ 0.66666667, 2.66666667],
1848                           [ 0.0,        1.0],
1849                           [ 0.0,        3.0],
1850                           [ 1.0,        0.0],
1851                           [ 1.0,        1.0],
1852                           [ 1.0,        2.0],
1853                           [ 1.0,        3.0],
1854                           [ 2.0,        1.0],
1855                           [ 3.0,        0.0],
1856                           [ 3.0,        1.0]]
1857
1858        data_geo_spatial = Geospatial_data(data_points_rel,
1859                                           geo_reference=Geo_reference(56,
1860                                                                       x0,
1861                                                                       y0))
1862        data_points_absolute = data_geo_spatial.get_data_points(absolute=True)
1863        attributes = linear_function(data_points_absolute)
1864        att = 'spam_and_eggs'
1865
1866        # Create .txt file
1867        ptsfile = tempfile.mktemp(".txt")
1868        file = open(ptsfile, "w")
1869        file.write(" x,y," + att + " \n")
1870        for data_point, attribute in map(None, data_points_absolute, attributes):
1871            row = (str(data_point[0]) + ',' +
1872                   str(data_point[1]) + ',' +
1873                   str(attribute))
1874            file.write(row + "\n")
1875        file.close()
1876
1877        # Check that values can be set from file
1878        quantity.set_values(filename=ptsfile, attribute_name=att, alpha=0)
1879        answer = linear_function(quantity.domain.get_vertex_coordinates())
1880
1881        assert num.allclose(quantity.vertex_values.flat, answer)
1882
1883        # Check that values can be set from file using default attribute
1884        quantity.set_values(filename = ptsfile, alpha = 0)
1885        assert num.allclose(quantity.vertex_values.flat, answer)
1886
1887        # Cleanup
1888        import os
1889        os.remove(ptsfile)
1890
1891    def test_fitting_example_that_crashed(self):
1892        """This unit test has been derived from a real world example
1893        (the Towradgi '98 rainstorm simulation).
1894
1895        It shows a condition where fitting as called from set_quantity crashes
1896        when ANUGA mesh is reused. The test passes in the case where a new mesh
1897        is created.
1898
1899        See ticket:314
1900        """
1901
1902        verbose = False
1903       
1904        # Get path where this test is run
1905        path = get_pathname_from_package('anuga.shallow_water')       
1906       
1907       
1908        #----------------------------------------------------------------------
1909        # Create domain
1910        #--------------------------------------------------------------------
1911        W = 303400
1912        N = 6195800
1913        E = 308640
1914        S = 6193120
1915        bounding_polygon = [[W, S], [E, S], [E, N], [W, N]]
1916
1917        offending_regions = []
1918
1919        # From culvert 8
1920        offending_regions.append([[307611.43896231, 6193631.6894806],
1921                                  [307600.11394969, 6193608.2855474],
1922                                  [307597.41349586, 6193609.59227963],
1923                                  [307608.73850848, 6193632.99621282]])
1924        offending_regions.append([[307633.69143231, 6193620.9216536],
1925                                  [307622.36641969, 6193597.5177204],
1926                                  [307625.06687352, 6193596.21098818],
1927                                  [307636.39188614, 6193619.61492137]])
1928
1929        # From culvert 9
1930        offending_regions.append([[306326.69660524, 6194818.62900522],
1931                                  [306324.67939476, 6194804.37099478],
1932                                  [306323.75856492, 6194804.50127295],
1933                                  [306325.7757754, 6194818.7592834]])
1934        offending_regions.append([[306365.57160524, 6194813.12900522],
1935                                  [306363.55439476, 6194798.87099478],
1936                                  [306364.4752246, 6194798.7407166],
1937                                  [306366.49243508, 6194812.99872705]])
1938
1939        # From culvert 10
1940        offending_regions.append([[306955.071019428608, 6194465.704096679576],
1941                                  [306951.616980571358, 6194457.295903320424],
1942                                  [306950.044491164153, 6194457.941873183474],
1943                                  [306953.498530021403, 6194466.350066542625]])
1944        offending_regions.append([[307002.540019428649, 6194446.204096679576],
1945                                  [306999.085980571399, 6194437.795903320424],
1946                                  [307000.658469978604, 6194437.149933457375],
1947                                  [307004.112508835853, 6194445.558126816526]])
1948
1949        interior_regions = []
1950        for polygon in offending_regions:
1951            interior_regions.append( [polygon, 100] ) 
1952
1953        meshname = os.path.join(path, 'offending_mesh.msh')
1954        anuga.create_mesh_from_regions(bounding_polygon,
1955                                 boundary_tags={'south': [0], 'east': [1],
1956                                                'north': [2], 'west': [3]},
1957                                 maximum_triangle_area=1000000,
1958                                 interior_regions=interior_regions,
1959                                 filename=meshname,
1960                                 use_cache=False,
1961                                 verbose=verbose)
1962
1963        domain = anuga.Domain(meshname, use_cache=False, verbose=verbose)
1964
1965        #----------------------------------------------------------------------
1966        # Fit data point to mesh
1967        #----------------------------------------------------------------------
1968
1969        points_file = os.path.join(path, 'offending_point.pts')
1970
1971        # Offending point
1972        G = Geospatial_data(data_points=[[306953.344, 6194461.5]],
1973                            attributes=[1])
1974        G.export_points_file(points_file)
1975       
1976        try:
1977            domain.set_quantity('elevation', filename=points_file,
1978                                use_cache=False, verbose=verbose, alpha=0.01)
1979        except RuntimeError, e:
1980            msg = 'Test failed: %s' % str(e)
1981            raise Exception, msg
1982            # clean up in case raise fails
1983            os.remove(meshname)
1984            os.remove(points_file)
1985        else:
1986            os.remove(meshname)
1987            os.remove(points_file)           
1988       
1989
1990    def test_fitting_example_that_crashed_2(self):
1991        """test_fitting_example_that_crashed_2
1992       
1993        This unit test has been derived from a real world example
1994        (the JJKelly study, by Petar Milevski).
1995       
1996        It shows a condition where set_quantity crashes due to AtA
1997        not being built properly
1998       
1999        See ticket:314
2000        """
2001
2002        verbose = False       
2003       
2004        # Get path where this test is run
2005        path = get_pathname_from_package('anuga.shallow_water')       
2006
2007        meshname = os.path.join(path, 'test_mesh.msh')
2008        W = 304180
2009        S = 6185270
2010        E = 307650
2011        N = 6189040
2012        maximum_triangle_area = 1000000
2013
2014        bounding_polygon = [[W, S], [E, S], [E, N], [W, N]]
2015
2016        anuga.create_mesh_from_regions(bounding_polygon,
2017                                 boundary_tags={'south': [0], 
2018                                                'east': [1], 
2019                                                'north': [2], 
2020                                                'west': [3]},
2021                                 maximum_triangle_area=maximum_triangle_area,
2022                                 filename=meshname,
2023                                 use_cache=False,
2024                                 verbose=verbose)
2025
2026        domain = anuga.Domain(meshname, use_cache=True, verbose=verbose)
2027       
2028        # Large test set revealed one problem
2029        points_file = os.path.join(path, 'test_points_large.csv')
2030        try:
2031            domain.set_quantity('elevation', filename=points_file,
2032                                use_cache=False, verbose=verbose)
2033        except AssertionError, e:
2034            msg = 'Test failed: %s' % str(e)
2035            raise Exception, msg
2036            # Cleanup in case this failed
2037            os.remove(meshname)
2038
2039        # Small test set revealed another problem
2040        points_file = os.path.join(path, 'test_points_small.csv')
2041        try:   
2042            domain.set_quantity('elevation', filename=points_file,
2043                                use_cache=False, verbose=verbose)                           
2044        except AssertionError, e:
2045            msg = 'Test failed: %s' % str(e)
2046            raise Exception, msg
2047            # Cleanup in case this failed
2048            os.remove(meshname)
2049        else:
2050            os.remove(meshname)
2051
2052
2053
2054
2055    def test_variable_elevation(self):           
2056        """test_variable_elevation
2057
2058        This will test that elevagtion van be stored in sww files
2059        as a time dependent quantity.
2060       
2061        It will also chck that storage of other quantities
2062        can be controlled this way.
2063        """
2064       
2065        #---------------------------------------------------------------------
2066        # Setup computational domain
2067        #---------------------------------------------------------------------
2068        length = 8.
2069        width = 6.
2070        dx = dy = 1    # Resolution: Length of subdivisions on both axes
2071       
2072        inc = 0.05 # Elevation increment
2073
2074        points, vertices, boundary = rectangular_cross(int(length/dx), 
2075                                                       int(width/dy),
2076                                                       len1=length, 
2077                                                       len2=width)
2078        domain = anuga.Domain(points, vertices, boundary)
2079        domain.set_name('channel_variable_test')  # Output name
2080        domain.set_quantities_to_be_stored({'elevation': 2,
2081                                            'stage': 2})
2082
2083        #---------------------------------------------------------------------
2084        # Setup initial conditions
2085        #---------------------------------------------------------------------
2086
2087        def pole_increment(x,y):
2088            """This provides a small increment to a pole located mid stream
2089            For use with variable elevation data
2090            """
2091           
2092            z = 0.0*x
2093
2094            N = len(x)
2095            for i in range(N):
2096                # Pole
2097                if (x[i] - 4)**2 + (y[i] - 2)**2 < 1.0**2:
2098                    z[i] += inc
2099            return z
2100           
2101        domain.set_quantity('elevation', 0.0)    # Flat bed initially
2102        domain.set_quantity('friction', 0.01)    # Constant friction
2103        domain.set_quantity('stage', 0.0)        # Dry initial condition
2104
2105        #------------------------------------------------------------------
2106        # Setup boundary conditions
2107        #------------------------------------------------------------------
2108        Bi = anuga.Dirichlet_boundary([0.4, 0, 0])          # Inflow
2109        Br = anuga.Reflective_boundary(domain)              # Solid reflective wall
2110        Bo = anuga.Dirichlet_boundary([-5, 0, 0])           # Outflow
2111
2112        domain.set_boundary({'left': Bi, 'right': Bo, 'top': Br, 'bottom': Br})
2113
2114        #-------------------------------------------------------------------
2115        # Evolve system through time
2116        #-------------------------------------------------------------------
2117
2118        for t in domain.evolve(yieldstep=1, finaltime=3.0):
2119            #print domain.timestepping_statistics()
2120
2121            domain.add_quantity('elevation', pole_increment)
2122       
2123           
2124        # Check that quantities have been stored correctly   
2125        from Scientific.IO.NetCDF import NetCDFFile
2126        sww_file = domain.get_name() + '.sww'
2127        fid = NetCDFFile(sww_file)
2128
2129        x = fid.variables['x'][:]
2130        y = fid.variables['y'][:]
2131        stage = fid.variables['stage'][:]
2132        elevation = fid.variables['elevation'][:]
2133        fid.close()
2134
2135        os.remove(sww_file)
2136       
2137                   
2138        assert len(stage.shape) == 2
2139        assert len(elevation.shape) == 2       
2140       
2141        M, N = stage.shape
2142               
2143        for i in range(M): 
2144            # For each timestep
2145            assert num.allclose(max(elevation[i,:]), i * inc) 
2146
2147
2148
2149#################################################################################
2150
2151if __name__ == "__main__":
2152    suite = unittest.makeSuite(Test_swb_basic, 'test')
2153    runner = unittest.TextTestRunner(verbosity=1)
2154    runner.run(suite)
Note: See TracBrowser for help on using the repository browser.