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

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

More swb tests passing. Cleaned up some pylint errors.

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