source: anuga_core/source/anuga/shallow_water_balanced/test_swb_basic.py @ 7559

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

Broke up test_swb_domain so that we can concentrate on new functionality

File size: 85.2 KB
Line 
1#!/usr/bin/env python
2
3import unittest, os
4import os.path
5from math import pi, sqrt
6import tempfile
7
8from anuga.config import g, epsilon
9from anuga.config import netcdf_mode_r, netcdf_mode_w, netcdf_mode_a
10from anuga.utilities.numerical_tools import mean
11from anuga.utilities.polygon import is_inside_polygon
12from anuga.coordinate_transforms.geo_reference import Geo_reference
13from anuga.abstract_2d_finite_volumes.quantity import Quantity
14from anuga.geospatial_data.geospatial_data import Geospatial_data
15from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular_cross
16
17from anuga.utilities.system_tools import get_pathname_from_package
18from swb_domain import *
19
20import numpy as num
21
22# Get gateway to C implementation of flux function for direct testing
23from shallow_water_ext import flux_function_central as flux_function
24
25
26# For test_fitting_using_shallow_water_domain example
27def linear_function(point):
28    point = num.array(point)
29    return point[:,0]+point[:,1]
30
31class Weir:
32    """Set a bathymetry for weir with a hole and a downstream gutter
33    x,y are assumed to be in the unit square
34    """
35
36    def __init__(self, stage):
37        self.inflow_stage = stage
38
39    def __call__(self, x, y):
40        N = len(x)
41        assert N == len(y)
42
43        z = num.zeros(N, num.float)
44        for i in range(N):
45            z[i] = -x[i]/2  #General slope
46
47            #Flattish bit to the left
48            if x[i] < 0.3:
49                z[i] = -x[i]/10
50
51            #Weir
52            if x[i] >= 0.3 and x[i] < 0.4:
53                z[i] = -x[i]+0.9
54
55            #Dip
56            x0 = 0.6
57            depth = -1.0
58            plateaux = -0.6
59            if y[i] < 0.7:
60                if x[i] > x0 and x[i] < 0.9:
61                    z[i] = depth
62                #RHS plateaux
63                if x[i] >= 0.9:
64                    z[i] = plateaux
65            elif y[i] >= 0.7 and y[i] < 1.5:
66                #Restrict and deepen
67                if x[i] >= x0 and x[i] < 0.8:
68                    z[i] = depth - (y[i]/3 - 0.3)
69                elif x[i] >= 0.8:
70                    #RHS plateaux
71                    z[i] = plateaux
72            elif y[i] >= 1.5:
73                if x[i] >= x0 and x[i] < 0.8 + (y[i]-1.5)/1.2:
74                    #Widen up and stay at constant depth
75                    z[i] = depth-1.5/5
76                elif x[i] >= 0.8 + (y[i]-1.5)/1.2:
77                    #RHS plateaux
78                    z[i] = plateaux
79
80            #Hole in weir (slightly higher than inflow condition)
81            if x[i] >= 0.3 and x[i] < 0.4 and y[i] > 0.2 and y[i] < 0.4:
82                z[i] = -x[i]+self.inflow_stage + 0.02
83
84            #Channel behind weir
85            x0 = 0.5
86            if x[i] >= 0.4 and x[i] < x0 and y[i] > 0.2 and y[i] < 0.4:
87                z[i] = -x[i]+self.inflow_stage + 0.02
88
89            if x[i] >= x0 and x[i] < 0.6 and y[i] > 0.2 and y[i] < 0.4:
90                #Flatten it out towards the end
91                z[i] = -x0+self.inflow_stage + 0.02 + (x0-x[i])/5
92
93            # Hole to the east
94            x0 = 1.1
95            y0 = 0.35
96            if num.sqrt((2*(x[i]-x0))**2 + (2*(y[i]-y0))**2) < 0.2:
97                z[i] = num.sqrt(((x[i]-x0))**2 + ((y[i]-y0))**2)-1.0
98
99            #Tiny channel draining hole
100            if x[i] >= 1.14 and x[i] < 1.2 and y[i] >= 0.4 and y[i] < 0.6:
101                z[i] = -0.9 #North south
102
103            if x[i] >= 0.9 and x[i] < 1.18 and y[i] >= 0.58 and y[i] < 0.65:
104                z[i] = -1.0 + (x[i]-0.9)/3 #East west
105
106            # Stuff not in use
107
108            # Upward slope at inlet to the north west
109            # if x[i] < 0.0: # and y[i] > 0.5:
110            #    #z[i] = -y[i]+0.5  #-x[i]/2
111            #    z[i] = x[i]/4 - y[i]**2 + 0.5
112
113            # Hole to the west
114            # x0 = -0.4; y0 = 0.35 # center
115            # if sqrt((2*(x[i]-x0))**2 + (2*(y[i]-y0))**2) < 0.2:
116            #    z[i] = sqrt(((x[i]-x0))**2 + ((y[i]-y0))**2)-0.2
117
118        return z/2
119
120
121
122#########################################################
123
124class Test_swb_basic(unittest.TestCase):
125    def setUp(self):
126        pass
127
128    def tearDown(self):
129        pass
130
131
132    def test_rotate(self):
133        normal = num.array([0.0, -1.0])
134
135        q = num.array([1.0, 2.0, 3.0])
136
137        r = rotate(q, normal, direction = 1)
138        assert r[0] == 1
139        assert r[1] == -3
140        assert r[2] == 2
141
142        w = rotate(r, normal, direction = -1)
143        assert num.allclose(w, q)
144
145        # Check error check
146        try:
147            rotate(r, num.array([1, 1, 1]))
148        except:
149            pass
150        else:
151            raise Exception, 'Should have raised an exception'
152
153    # Individual flux tests
154    def test_flux_zero_case(self):
155        ql = num.zeros(3, num.float)
156        qr = num.zeros(3, num.float)
157        normal = num.zeros(2, num.float)
158        edgeflux = num.zeros(3, num.float)
159        zl = zr = 0.
160        H0 = 1.0e-3 # As suggested in the manual
161       
162        max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux, epsilon, g, H0)
163
164        assert num.allclose(edgeflux, [0,0,0])
165        assert max_speed == 0.
166
167    def test_flux_constants(self):
168        w = 2.0
169
170        normal = num.array([1.,0])
171        ql = num.array([w, 0, 0])
172        qr = num.array([w, 0, 0])
173        edgeflux = num.zeros(3, num.float)
174        zl = zr = 0.
175        h = w - (zl+zr)/2
176        H0 = 0.0
177
178        max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux, epsilon, g, H0)       
179        assert num.allclose(edgeflux, [0., 0.5*g*h**2, 0.])
180        assert max_speed == num.sqrt(g*h)
181
182    #def test_flux_slope(self):
183    #    #FIXME: TODO
184    #    w = 2.0
185    #
186    #    normal = array([1.,0])
187    #    ql = array([w, 0, 0])
188    #    qr = array([w, 0, 0])
189    #    zl = zr = 0.
190    #    h = w - (zl+zr)/2
191    #
192    #    flux, max_speed = flux_function(normal, ql, qr, zl, zr)
193    #
194    #    assert allclose(flux, [0., 0.5*g*h**2, 0.])
195    #    assert max_speed == sqrt(g*h)
196
197    def test_flux1(self):
198        # Use data from previous version of abstract_2d_finite_volumes
199        normal = num.array([1., 0])
200        ql = num.array([-0.2, 2, 3])
201        qr = num.array([-0.2, 2, 3])
202        zl = zr = -0.5
203        edgeflux = num.zeros(3, num.float)
204
205        H0 = 0.0
206
207        max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux,
208                                  epsilon, g, H0)
209
210        assert num.allclose(edgeflux, [2., 13.77433333, 20.])
211        assert num.allclose(max_speed, 8.38130948661)
212
213    def test_flux2(self):
214        # Use data from previous version of abstract_2d_finite_volumes
215        normal = num.array([0., -1.])
216        ql = num.array([-0.075, 2, 3])
217        qr = num.array([-0.075, 2, 3])
218        zl = zr = -0.375
219
220        edgeflux = num.zeros(3, num.float)
221        H0 = 0.0
222        max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux,
223                                  epsilon, g, H0)
224
225        assert num.allclose(edgeflux, [-3., -20.0, -30.441])
226        assert num.allclose(max_speed, 11.7146428199)
227
228    def test_flux3(self):
229        # Use data from previous version of abstract_2d_finite_volumes
230        normal = num.array([-sqrt(2)/2, sqrt(2)/2])
231        ql = num.array([-0.075, 2, 3])
232        qr = num.array([-0.075, 2, 3])
233        zl = zr = -0.375
234
235        edgeflux = num.zeros(3, num.float)
236        H0 = 0.0
237        max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux,
238                                  epsilon, g, H0)
239
240        assert num.allclose(edgeflux, [sqrt(2)/2, 4.40221112, 7.3829019])
241        assert num.allclose(max_speed, 4.0716654239)
242
243    def test_flux4(self):
244        # Use data from previous version of abstract_2d_finite_volumes
245        normal = num.array([-sqrt(2)/2, sqrt(2)/2])
246        ql = num.array([-0.34319278, 0.10254161, 0.07273855])
247        qr = num.array([-0.30683287, 0.1071986, 0.05930515])
248        zl = zr = -0.375
249
250        edgeflux = num.zeros(3, num.float)
251        H0 = 0.0
252        max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux,
253                                  epsilon, g, H0)
254
255        assert num.allclose(edgeflux, [-0.04072676, -0.07096636, -0.01604364])
256        assert num.allclose(max_speed, 1.31414103233)
257
258    def test_flux_computation(self):
259        """test flux calculation (actual C implementation)
260
261        This one tests the constant case where only the pressure term
262        contributes to each edge and cancels out once the total flux has
263        been summed up.
264        """
265
266        a = [0.0, 0.0]
267        b = [0.0, 2.0]
268        c = [2.0, 0.0]
269        d = [0.0, 4.0]
270        e = [2.0, 2.0]
271        f = [4.0, 0.0]
272
273        points = [a, b, c, d, e, f]
274        #              bac,     bce,     ecf,     dbe
275        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
276
277        domain = Domain(points, vertices)
278        domain.check_integrity()
279
280        # The constant case
281        domain.set_quantity('elevation', -1)
282        domain.set_quantity('stage', 1)
283
284        domain.compute_fluxes()
285        # Central triangle
286        assert num.allclose(domain.get_quantity('stage').explicit_update[1], 0)
287
288        # The more general case
289        def surface(x, y):
290            return -x/2
291
292        domain.set_quantity('elevation', -10)
293        domain.set_quantity('stage', surface)
294        domain.set_quantity('xmomentum', 1)
295
296        domain.compute_fluxes()
297
298        #print domain.get_quantity('stage').explicit_update
299        # FIXME (Ole): TODO the general case
300        #assert allclose(domain.get_quantity('stage').explicit_update[1], ...??)
301
302    def test_sw_domain_simple(self):
303        a = [0.0, 0.0]
304        b = [0.0, 2.0]
305        c = [2.0, 0.0]
306        d = [0.0, 4.0]
307        e = [2.0, 2.0]
308        f = [4.0, 0.0]
309
310        points = [a, b, c, d, e, f]
311        #             bac,     bce,     ecf,     dbe
312        vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
313
314        #from anuga.abstract_2d_finite_volumes.domain import Domain as Generic_domain
315        #msg = 'The class %s is not a subclass of the generic domain class %s'\
316        #      %(DomainClass, Domain)
317        #assert issubclass(DomainClass, Domain), msg
318
319        domain = Domain(points, vertices)
320        domain.check_integrity()
321
322        for name in ['stage', 'xmomentum', 'ymomentum',
323                     'elevation', 'friction']:
324            assert domain.quantities.has_key(name)
325
326        assert num.alltrue(domain.get_conserved_quantities(0, edge=1) == 0.)
327
328    def xtest_catching_negative_heights(self):
329        #OBSOLETE
330
331        a = [0.0, 0.0]
332        b = [0.0, 2.0]
333        c = [2.0, 0.0]
334        d = [0.0, 4.0]
335        e = [2.0, 2.0]
336        f = [4.0, 0.0]
337
338        points = [a, b, c, d, e, f]
339        #             bac,     bce,     ecf,     dbe
340        vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
341
342        domain = Domain(points, vertices)
343        val0 = 2. + 2.0/3
344        val1 = 4. + 4.0/3
345        val2 = 8. + 2.0/3
346        val3 = 2. + 8.0/3
347
348        zl = zr = 4    # Too large
349        domain.set_quantity('elevation', zl*num.ones((4, 3), num.int)) #array default#
350        domain.set_quantity('stage', [[val0, val0-1, val0-2],
351                                      [val1, val1+1, val1],
352                                      [val2, val2-2, val2],
353                                      [val3-0.5, val3, val3]])
354
355        #Should fail
356        try:
357            domain.check_integrity()
358        except:
359            pass
360
361    def test_get_wet_elements(self):
362        a = [0.0, 0.0]
363        b = [0.0, 2.0]
364        c = [2.0, 0.0]
365        d = [0.0, 4.0]
366        e = [2.0, 2.0]
367        f = [4.0, 0.0]
368
369        points = [a, b, c, d, e, f]
370        #             bac,     bce,     ecf,     dbe
371        vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
372
373        domain = Domain(points, vertices)
374
375        val0 = 2. + 2.0/3
376        val1 = 4. + 4.0/3
377        val2 = 8. + 2.0/3
378        val3 = 2. + 8.0/3
379
380        zl = zr = 5
381        domain.set_quantity('elevation', zl*num.ones((4, 3), num.int)) #array default#
382        domain.set_quantity('stage', [[val0, val0-1, val0-2],
383                                      [val1, val1+1, val1],
384                                      [val2, val2-2, val2],
385                                      [val3-0.5, val3, val3]])
386
387        domain.check_integrity()
388
389        indices = domain.get_wet_elements()
390        assert num.allclose(indices, [1, 2])
391
392        indices = domain.get_wet_elements(indices=[0, 1, 3])
393        assert num.allclose(indices, [1])
394
395    def test_get_maximum_inundation_1(self):
396        a = [0.0, 0.0]
397        b = [0.0, 2.0]
398        c = [2.0, 0.0]
399        d = [0.0, 4.0]
400        e = [2.0, 2.0]
401        f = [4.0, 0.0]
402
403        points = [a, b, c, d, e, f]
404        #             bac,     bce,     ecf,     dbe
405        vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
406
407        domain = Domain(points, vertices)
408
409        domain.set_quantity('elevation', lambda x, y: x+2*y)    # 2 4 4 6
410        domain.set_quantity('stage', 3)
411
412        domain.check_integrity()
413
414        indices = domain.get_wet_elements()
415        assert num.allclose(indices, [0])
416
417        q = domain.get_maximum_inundation_elevation()
418        assert num.allclose(q, domain.get_quantity('elevation').\
419                                   get_values(location='centroids')[0])
420
421        x, y = domain.get_maximum_inundation_location()
422        assert num.allclose([x, y], domain.get_centroid_coordinates()[0])
423
424    def test_get_maximum_inundation_2(self):
425        """test_get_maximum_inundation_2(self)
426
427        Test multiple wet cells with same elevation
428        """
429
430        a = [0.0, 0.0]
431        b = [0.0, 2.0]
432        c = [2.0, 0.0]
433        d = [0.0, 4.0]
434        e = [2.0, 2.0]
435        f = [4.0, 0.0]
436
437        points = [a, b, c, d, e, f]
438        #             bac,     bce,     ecf,     dbe
439        vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
440
441        domain = Domain(points, vertices)
442
443        domain.set_quantity('elevation', lambda x, y: x+2*y)    # 2 4 4 6
444        domain.set_quantity('stage', 4.1)
445
446        domain.check_integrity()
447
448        indices = domain.get_wet_elements()
449        assert num.allclose(indices, [0, 1, 2])
450
451        q = domain.get_maximum_inundation_elevation()
452        assert num.allclose(q, 4)
453
454        x, y = domain.get_maximum_inundation_location()
455        assert num.allclose([x, y], domain.get_centroid_coordinates()[1])
456
457    def test_get_maximum_inundation_3(self):
458        """test_get_maximum_inundation_3(self)
459
460        Test of real runup example:
461        """
462
463        from anuga.abstract_2d_finite_volumes.mesh_factory \
464                import rectangular_cross
465
466        initial_runup_height = -0.4
467        final_runup_height = -0.3
468
469        #--------------------------------------------------------------
470        # Setup computational domain
471        #--------------------------------------------------------------
472        N = 5
473        points, vertices, boundary = rectangular_cross(N, N)
474        domain = Domain(points, vertices, boundary)
475        domain.set_maximum_allowed_speed(1.0)
476
477        #--------------------------------------------------------------
478        # Setup initial conditions
479        #--------------------------------------------------------------
480        def topography(x, y):
481            return -x/2                             # linear bed slope
482
483        # Use function for elevation
484        domain.set_quantity('elevation', topography)
485        domain.set_quantity('friction', 0.)                # Zero friction
486        # Constant negative initial stage
487        domain.set_quantity('stage', initial_runup_height)
488
489        #--------------------------------------------------------------
490        # Setup boundary conditions
491        #--------------------------------------------------------------
492        Br = Reflective_boundary(domain)                       # Reflective wall
493        Bd = Dirichlet_boundary([final_runup_height, 0, 0])    # Constant inflow
494
495        # All reflective to begin with (still water)
496        domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
497
498        #--------------------------------------------------------------
499        # Test initial inundation height
500        #--------------------------------------------------------------
501
502        indices = domain.get_wet_elements()
503        z = domain.get_quantity('elevation').get_values(location='centroids',
504                                                        indices=indices)
505        assert num.alltrue(z < initial_runup_height)
506
507        q = domain.get_maximum_inundation_elevation()
508        # First order accuracy
509        assert num.allclose(q, initial_runup_height, rtol=1.0/N)
510
511        x, y = domain.get_maximum_inundation_location()
512
513        qref = domain.get_quantity('elevation').\
514                     get_values(interpolation_points=[[x, y]])
515        assert num.allclose(q, qref)
516
517        wet_elements = domain.get_wet_elements()
518        wet_elevations = domain.get_quantity('elevation').\
519                                    get_values(location='centroids',
520                                               indices=wet_elements)
521        assert num.alltrue(wet_elevations < initial_runup_height)
522        assert num.allclose(wet_elevations, z)
523
524        #--------------------------------------------------------------
525        # Let triangles adjust
526        #--------------------------------------------------------------
527        for t in domain.evolve(yieldstep = 0.1, finaltime = 1.0):
528            pass
529
530        #--------------------------------------------------------------
531        # Test inundation height again
532        #--------------------------------------------------------------
533        indices = domain.get_wet_elements()
534        z = domain.get_quantity('elevation').get_values(location='centroids',
535                                                        indices=indices)
536
537        assert num.alltrue(z < initial_runup_height)
538
539        q = domain.get_maximum_inundation_elevation()
540        # First order accuracy
541        assert num.allclose(q, initial_runup_height, rtol=1.0/N)
542
543        x, y = domain.get_maximum_inundation_location()
544        qref = domain.get_quantity('elevation').\
545                        get_values(interpolation_points=[[x, y]])
546        assert num.allclose(q, qref)
547
548        #--------------------------------------------------------------
549        # Update boundary to allow inflow
550        #--------------------------------------------------------------
551        domain.set_boundary({'right': Bd})
552
553        #--------------------------------------------------------------
554        # Evolve system through time
555        #--------------------------------------------------------------
556        for t in domain.evolve(yieldstep = 0.1, finaltime = 3.0):
557            pass
558
559        #--------------------------------------------------------------
560        # Test inundation height again
561        #--------------------------------------------------------------
562        indices = domain.get_wet_elements()
563        z = domain.get_quantity('elevation').\
564                    get_values(location='centroids', indices=indices)
565
566        assert num.alltrue(z < final_runup_height)
567
568        q = domain.get_maximum_inundation_elevation()
569        # First order accuracy
570        assert num.allclose(q, final_runup_height, rtol=1.0/N)
571
572        x, y = domain.get_maximum_inundation_location()
573        qref = domain.get_quantity('elevation').\
574                        get_values(interpolation_points=[[x, y]])
575        assert num.allclose(q, qref)
576
577        wet_elements = domain.get_wet_elements()
578        wet_elevations = domain.get_quantity('elevation').\
579                             get_values(location='centroids',
580                                        indices=wet_elements)
581        assert num.alltrue(wet_elevations < final_runup_height)
582        assert num.allclose(wet_elevations, z)
583
584    def test_get_maximum_inundation_from_sww(self):
585        """test_get_maximum_inundation_from_sww(self)
586
587        Test of get_maximum_inundation_elevation()
588        and get_maximum_inundation_location() from data_manager.py
589
590        This is based on test_get_maximum_inundation_3(self) but works with the
591        stored results instead of with the internal data structure.
592
593        This test uses the underlying get_maximum_inundation_data for tests
594        """
595
596        from anuga.abstract_2d_finite_volumes.mesh_factory \
597                import rectangular_cross
598        from data_manager import get_maximum_inundation_elevation
599        from data_manager import get_maximum_inundation_location
600        from data_manager import get_maximum_inundation_data
601
602        initial_runup_height = -0.4
603        final_runup_height = -0.3
604
605        #--------------------------------------------------------------
606        # Setup computational domain
607        #--------------------------------------------------------------
608        N = 10
609        points, vertices, boundary = rectangular_cross(N, N)
610        domain = Domain(points, vertices, boundary)
611        domain.set_name('runup_test')
612        domain.set_maximum_allowed_speed(1.0)
613
614        # FIXME: This works better with old limiters so far
615        domain.tight_slope_limiters = 0
616
617        #--------------------------------------------------------------
618        # Setup initial conditions
619        #--------------------------------------------------------------
620        def topography(x, y):
621            return -x/2                             # linear bed slope
622
623        # Use function for elevation
624        domain.set_quantity('elevation', topography)
625        domain.set_quantity('friction', 0.)                # Zero friction
626        # Constant negative initial stage
627        domain.set_quantity('stage', initial_runup_height)
628
629        #--------------------------------------------------------------
630        # Setup boundary conditions
631        #--------------------------------------------------------------
632        Br = Reflective_boundary(domain)                       # Reflective wall
633        Bd = Dirichlet_boundary([final_runup_height, 0, 0])    # Constant inflow
634
635        # All reflective to begin with (still water)
636        domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
637
638        #--------------------------------------------------------------
639        # Test initial inundation height
640        #--------------------------------------------------------------
641        indices = domain.get_wet_elements()
642        z = domain.get_quantity('elevation').\
643                get_values(location='centroids', indices=indices)
644        assert num.alltrue(z < initial_runup_height)
645
646        q_ref = domain.get_maximum_inundation_elevation()
647        # First order accuracy
648        assert num.allclose(q_ref, initial_runup_height, rtol=1.0/N)
649
650        #--------------------------------------------------------------
651        # Let triangles adjust
652        #--------------------------------------------------------------
653        for t in domain.evolve(yieldstep = 0.1, finaltime = 1.0):
654            pass
655
656        #--------------------------------------------------------------
657        # Test inundation height again
658        #--------------------------------------------------------------
659        q_ref = domain.get_maximum_inundation_elevation()
660        q = get_maximum_inundation_elevation('runup_test.sww')
661        msg = 'We got %f, should have been %f' % (q, q_ref)
662        assert num.allclose(q, q_ref, rtol=1.0/N), msg
663
664        q = get_maximum_inundation_elevation('runup_test.sww')
665        msg = 'We got %f, should have been %f' % (q, initial_runup_height)
666        assert num.allclose(q, initial_runup_height, rtol = 1.0/N), msg
667
668        # Test error condition if time interval is out
669        try:
670            q = get_maximum_inundation_elevation('runup_test.sww',
671                                                 time_interval=[2.0, 3.0])
672        except ValueError:
673            pass
674        else:
675            msg = 'should have caught wrong time interval'
676            raise Exception, msg
677
678        # Check correct time interval
679        q, loc = get_maximum_inundation_data('runup_test.sww',
680                                             time_interval=[0.0, 3.0])
681        msg = 'We got %f, should have been %f' % (q, initial_runup_height)
682        assert num.allclose(q, initial_runup_height, rtol = 1.0/N), msg
683        assert num.allclose(-loc[0]/2, q)    # From topography formula
684
685        #--------------------------------------------------------------
686        # Update boundary to allow inflow
687        #--------------------------------------------------------------
688        domain.set_boundary({'right': Bd})
689
690        #--------------------------------------------------------------
691        # Evolve system through time
692        #--------------------------------------------------------------
693        q_max = None
694        for t in domain.evolve(yieldstep = 0.1, finaltime = 3.0,
695                               skip_initial_step = True):
696            q = domain.get_maximum_inundation_elevation()
697            if q > q_max:
698                q_max = q
699
700        #--------------------------------------------------------------
701        # Test inundation height again
702        #--------------------------------------------------------------
703        indices = domain.get_wet_elements()
704        z = domain.get_quantity('elevation').\
705                get_values(location='centroids', indices=indices)
706
707        assert num.alltrue(z < final_runup_height)
708
709        q = domain.get_maximum_inundation_elevation()
710        # First order accuracy
711        assert num.allclose(q, final_runup_height, rtol=1.0/N)
712
713        q, loc = get_maximum_inundation_data('runup_test.sww',
714                                             time_interval=[3.0, 3.0])
715        msg = 'We got %f, should have been %f' % (q, final_runup_height)
716        assert num.allclose(q, final_runup_height, rtol=1.0/N), msg
717        assert num.allclose(-loc[0]/2, q)    # From topography formula
718
719        q = get_maximum_inundation_elevation('runup_test.sww')
720        loc = get_maximum_inundation_location('runup_test.sww')
721        msg = 'We got %f, should have been %f' % (q, q_max)
722        assert num.allclose(q, q_max, rtol=1.0/N), msg
723        assert num.allclose(-loc[0]/2, q)    # From topography formula
724
725        q = get_maximum_inundation_elevation('runup_test.sww',
726                                             time_interval=[0, 3])
727        msg = 'We got %f, should have been %f' % (q, q_max)
728        assert num.allclose(q, q_max, rtol=1.0/N), msg
729
730        # Check polygon mode
731        # Runup region
732        polygon = [[0.3, 0.0], [0.9, 0.0], [0.9, 1.0], [0.3, 1.0]]
733        q = get_maximum_inundation_elevation('runup_test.sww',
734                                             polygon = polygon,
735                                             time_interval=[0, 3])
736        msg = 'We got %f, should have been %f' % (q, q_max)
737        assert num.allclose(q, q_max, rtol=1.0/N), msg
738
739        # Offshore region
740        polygon = [[0.9, 0.0], [1.0, 0.0], [1.0, 1.0], [0.9, 1.0]]
741        q, loc = get_maximum_inundation_data('runup_test.sww',
742                                             polygon = polygon,
743                                             time_interval=[0, 3])
744        msg = 'We got %f, should have been %f' % (q, -0.475)
745        assert num.allclose(q, -0.475, rtol=1.0/N), msg
746        assert is_inside_polygon(loc, polygon)
747        assert num.allclose(-loc[0]/2, q)    # From topography formula
748
749        # Dry region
750        polygon = [[0.0, 0.0], [0.4, 0.0], [0.4, 1.0], [0.0, 1.0]]
751        q, loc = get_maximum_inundation_data('runup_test.sww',
752                                             polygon = polygon,
753                                             time_interval=[0, 3])
754        msg = 'We got %s, should have been None' % (q)
755        assert q is None, msg
756        msg = 'We got %s, should have been None' % (loc)
757        assert loc is None, msg
758
759        # Check what happens if no time point is within interval
760        try:
761            q = get_maximum_inundation_elevation('runup_test.sww',
762                                                 time_interval=[2.75, 2.75])
763        except AssertionError:
764            pass
765        else:
766            msg = 'Time interval should have raised an exception'
767            raise Exception, msg
768
769        # Cleanup
770        try:
771            os.remove(domain.get_name() + '.sww')
772        except:
773            pass
774            #FIXME(Ole): Windows won't allow removal of this
775
776
777
778
779
780    def test_another_runup_example(self):
781        """test_another_runup_example
782
783        Test runup example where actual timeseries at interpolated
784        points are tested.
785        """
786
787        from anuga.pmesh.mesh_interface import create_mesh_from_regions
788        from anuga.abstract_2d_finite_volumes.mesh_factory \
789                import rectangular_cross
790        from anuga.shallow_water import Domain
791        from anuga.shallow_water import Reflective_boundary
792        from anuga.shallow_water import Dirichlet_boundary
793
794        #-----------------------------------------------------------------
795        # Setup computational domain
796        #-----------------------------------------------------------------
797        points, vertices, boundary = rectangular_cross(10, 10) # Basic mesh
798        domain = Domain(points, vertices, boundary) # Create domain
799        domain.set_default_order(2)
800        domain.set_quantities_to_be_stored(None)
801        domain.H0 = 1.0e-3
802
803        #-----------------------------------------------------------------
804        # Setup initial conditions
805        #-----------------------------------------------------------------
806        def topography(x, y):
807            return -x/2                              # linear bed slope
808
809        domain.set_quantity('elevation', topography)
810        domain.set_quantity('friction', 0.0)
811        domain.set_quantity('stage', expression='elevation')
812
813        #----------------------------------------------------------------
814        # Setup boundary conditions
815        #----------------------------------------------------------------
816        Br = Reflective_boundary(domain)           # Solid reflective wall
817        Bd = Dirichlet_boundary([-0.2, 0., 0.])    # Constant boundary values
818        domain.set_boundary({'left': Br, 'right': Bd, 'top': Br, 'bottom': Br})
819
820        #----------------------------------------------------------------
821        # Evolve system through time
822        #----------------------------------------------------------------
823        interpolation_points = [[0.4,0.5], [0.6,0.5], [0.8,0.5], [0.9,0.5]]
824        gauge_values = []
825        for _ in interpolation_points:
826            gauge_values.append([])
827
828        time = []
829        for t in domain.evolve(yieldstep=0.1, finaltime=5.0):
830            # Record time series at known points
831            time.append(domain.get_time())
832
833            stage = domain.get_quantity('stage')
834            w = stage.get_values(interpolation_points=interpolation_points)
835
836            for i, _ in enumerate(interpolation_points):
837                gauge_values[i].append(w[i])
838
839        #Reference (nautilus 26/6/2008)
840       
841        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]
842
843        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]
844       
845       
846        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]
847       
848        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]
849       
850        assert num.allclose(gauge_values[0], G0)
851        assert num.allclose(gauge_values[1], G1)
852        assert num.allclose(gauge_values[2], G2)
853        assert num.allclose(gauge_values[3], G3)
854
855    #####################################################
856
857
858    def test_initial_condition(self):
859        """test_initial_condition
860
861        Test that initial condition is output at time == 0 and that
862        computed values change as system evolves
863        """
864
865        from anuga.config import g
866        import copy
867
868        a = [0.0, 0.0]
869        b = [0.0, 2.0]
870        c = [2.0, 0.0]
871        d = [0.0, 4.0]
872        e = [2.0, 2.0]
873        f = [4.0, 0.0]
874
875        points = [a, b, c, d, e, f]
876        #             bac,     bce,     ecf,     dbe
877        vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
878
879        domain = Domain(points, vertices)
880
881        #Set up for a gradient of (3,0) at mid triangle (bce)
882        def slope(x, y):
883            return 3*x
884
885        h = 0.1
886        def stage(x, y):
887            return slope(x, y) + h
888
889        domain.set_quantity('elevation', slope)
890        domain.set_quantity('stage', stage)
891
892        # Allow slope limiters to work
893        # (FIXME (Ole): Shouldn't this be automatic in ANUGA?)
894        domain.distribute_to_vertices_and_edges()
895
896        initial_stage = copy.copy(domain.quantities['stage'].vertex_values)
897
898        domain.set_boundary({'exterior': Reflective_boundary(domain)})
899
900        domain.optimise_dry_cells = True
901
902        #Evolution
903        for t in domain.evolve(yieldstep=0.5, finaltime=2.0):
904            stage = domain.quantities['stage'].vertex_values
905
906            if t == 0.0:
907                assert num.allclose(stage, initial_stage)
908            else:
909                assert not num.allclose(stage, initial_stage)
910
911        os.remove(domain.get_name() + '.sww')
912
913    #####################################################
914
915    def test_second_order_flat_bed_onestep(self):
916        from mesh_factory import rectangular
917
918        #Create basic mesh
919        points, vertices, boundary = rectangular(6, 6)
920
921        #Create shallow water domain
922        domain = Domain(points, vertices, boundary)
923        domain.smooth = False
924        domain.default_order = 2
925        domain.beta_w = 0.9
926        domain.beta_w_dry = 0.9
927        domain.beta_uh = 0.9
928        domain.beta_uh_dry = 0.9
929        domain.beta_vh = 0.9
930        domain.beta_vh_dry = 0.9
931        domain.H0 = 1.0e-3
932
933        # Boundary conditions
934        Br = Reflective_boundary(domain)
935        Bd = Dirichlet_boundary([0.1, 0., 0.])
936        domain.set_boundary({'left': Bd, 'right': Br, 'top': Br, 'bottom': Br})
937
938        domain.check_integrity()
939
940        # Evolution
941        for t in domain.evolve(yieldstep=0.05, finaltime=0.05):
942            pass
943
944        # Data from earlier version of abstract_2d_finite_volumes
945        assert num.allclose(domain.recorded_min_timestep, 0.0396825396825)
946        assert num.allclose(domain.recorded_max_timestep, 0.0396825396825)
947
948        msg = ("domain.quantities['stage'].centroid_values[:12]=%s"
949               % str(domain.quantities['stage'].centroid_values[:12]))
950        assert num.allclose(domain.quantities['stage'].centroid_values[:12],
951                            [0.00171396, 0.02656103, 0.00241523, 0.02656103,
952                             0.00241523, 0.02656103, 0.00241523, 0.02656103,
953                             0.00241523, 0.02656103, 0.00241523, 0.0272623]), msg
954
955        domain.distribute_to_vertices_and_edges()
956
957        assert num.allclose(domain.quantities['stage'].vertex_values[:12,0],
958                            [0.001, 0.02656103, 0.001, 0.02656103, 0.001, 0.02656103,
959                             0.001, 0.02656103, 0.001, 0.02656103, 0.001, 0.0272623])     
960                             
961        assert num.allclose(domain.quantities['stage'].vertex_values[:12,1],
962                            [0.00237867, 0.02656103, 0.001, 0.02656103, 0.001,
963                             0.02656103, 0.001, 0.02656103, 0.001, 0.02656103, 
964                             0.00110647, 0.0272623])
965
966        assert num.allclose(domain.quantities['stage'].vertex_values[:12,2],
967                            [0.00176321, 0.02656103, 0.00524568, 
968                             0.02656103, 0.00524568, 0.02656103,
969                             0.00524568, 0.02656103, 0.00524568,
970                             0.02656103, 0.00513921, 0.0272623])
971
972        assert num.allclose(domain.quantities['xmomentum'].centroid_values[:12],
973                            [0.00113961, 0.01302432, 0.00148672,
974                             0.01302432, 0.00148672, 0.01302432,
975                             0.00148672, 0.01302432, 0.00148672 ,
976                             0.01302432, 0.00148672, 0.01337143])
977
978        assert num.allclose(domain.quantities['ymomentum'].centroid_values[:12],
979                            [-2.91240050e-004 , 1.22721531e-004,
980                             -1.22721531e-004,  1.22721531e-004 ,
981                             -1.22721531e-004,  1.22721531e-004,
982                             -1.22721531e-004 , 1.22721531e-004,
983                             -1.22721531e-004,  1.22721531e-004,
984                             -1.22721531e-004,  -4.57969873e-005])
985
986        os.remove(domain.get_name() + '.sww')
987
988    def test_second_order_flat_bed_moresteps(self):
989        from mesh_factory import rectangular
990
991        # Create basic mesh
992        points, vertices, boundary = rectangular(6, 6)
993
994        # Create shallow water domain
995        domain = Domain(points, vertices, boundary)
996        domain.smooth = False
997        domain.default_order = 2
998
999        # Boundary conditions
1000        Br = Reflective_boundary(domain)
1001        Bd = Dirichlet_boundary([0.1, 0., 0.])
1002        domain.set_boundary({'left': Bd, 'right': Br, 'top': Br, 'bottom': Br})
1003
1004        domain.check_integrity()
1005
1006        # Evolution
1007        for t in domain.evolve(yieldstep=0.05, finaltime=0.1):
1008            pass
1009
1010        # Data from earlier version of abstract_2d_finite_volumes
1011        #assert allclose(domain.recorded_min_timestep, 0.0396825396825)
1012        #assert allclose(domain.recorded_max_timestep, 0.0396825396825)
1013        #print domain.quantities['stage'].centroid_values
1014
1015        os.remove(domain.get_name() + '.sww')
1016
1017    def test_flatbed_first_order(self):
1018        from mesh_factory import rectangular
1019
1020        # Create basic mesh
1021        N = 8
1022        points, vertices, boundary = rectangular(N, N)
1023
1024        # Create shallow water domain
1025        domain = Domain(points, vertices, boundary)
1026        domain.smooth = False
1027        domain.default_order = 1
1028        domain.H0 = 1.0e-3 # As suggested in the manual
1029
1030        # Boundary conditions
1031        Br = Reflective_boundary(domain)
1032        Bd = Dirichlet_boundary([0.2, 0., 0.])
1033
1034        domain.set_boundary({'left': Bd, 'right': Br, 'top': Br, 'bottom': Br})
1035        domain.check_integrity()
1036
1037        # Evolution
1038        for t in domain.evolve(yieldstep=0.02, finaltime=0.5):
1039            pass
1040
1041        # FIXME: These numbers were from version before 25/10
1042        #assert allclose(domain.recorded_min_timestep, 0.0140413643926)
1043        #assert allclose(domain.recorded_max_timestep, 0.0140947355753)
1044
1045        for i in range(3):
1046            #assert allclose(domain.quantities['stage'].edge_values[:4,i],
1047            #                [0.10730244,0.12337617,0.11200126,0.12605666])
1048            assert num.allclose(domain.quantities['xmomentum'].\
1049                                    edge_values[:4,i],
1050                                [0.07610894,0.06901572,0.07284461,0.06819712])
1051            #assert allclose(domain.quantities['ymomentum'].edge_values[:4,i],
1052            #                [-0.0060238, -0.00157404, -0.00309633, -0.0001637])
1053
1054        os.remove(domain.get_name() + '.sww')
1055
1056    def test_flatbed_second_order(self):
1057        from mesh_factory import rectangular
1058
1059        # Create basic mesh
1060        N = 8
1061        points, vertices, boundary = rectangular(N, N)
1062
1063        # Create shallow water domain
1064        domain = Domain(points, vertices, boundary)
1065        domain.set_store_vertices_uniquely(True)
1066        domain.set_default_order(2)
1067
1068        # Boundary conditions
1069        Br = Reflective_boundary(domain)
1070        Bd = Dirichlet_boundary([0.2, 0., 0.])
1071
1072        domain.set_boundary({'left': Bd, 'right': Br, 'top': Br, 'bottom': Br})
1073        domain.check_integrity()
1074
1075        # Evolution
1076        for t in domain.evolve(yieldstep=0.01, finaltime=0.03):
1077            pass
1078
1079        msg = 'min step was %f instead of %f' % (domain.recorded_min_timestep,
1080                                                 0.0210448446782)
1081
1082        assert num.allclose(domain.recorded_min_timestep, 0.0210448446782), msg
1083        assert num.allclose(domain.recorded_max_timestep, 0.0210448446782)
1084
1085        # Slight change due to flux limiter optimisation 28/5/9
1086        assert num.allclose(domain.quantities['stage'].vertex_values[:4,0],
1087                            [0.001, 0.05350737,  0.001, 0.0535293])
1088       
1089        assert num.allclose(domain.quantities['xmomentum'].vertex_values[:4,0],
1090                            [0.00090268, 0.03684904, 0.00084545, 0.03686323])
1091
1092        assert num.allclose(domain.quantities['ymomentum'].vertex_values[:4,0],
1093        [ -1.97318203e-04, 6.10268320e-04, -6.18053986e-05, 6.14082609e-04])
1094
1095        os.remove(domain.get_name() + '.sww')
1096
1097
1098    def test_flatbed_second_order_vmax_0(self):
1099        from mesh_factory import rectangular
1100
1101        # Create basic mesh
1102        N = 8
1103        points, vertices, boundary = rectangular(N, N)
1104
1105        # Create shallow water domain
1106        domain = Domain(points, vertices, boundary)
1107
1108        domain.set_store_vertices_uniquely(True)
1109        domain.set_default_order(2)
1110
1111
1112        # Boundary conditions
1113        Br = Reflective_boundary(domain)
1114        Bd = Dirichlet_boundary([0.2, 0., 0.])
1115
1116        domain.set_boundary({'left': Bd, 'right': Br, 'top': Br, 'bottom': Br})
1117        domain.check_integrity()
1118
1119        # Evolution
1120        for t in domain.evolve(yieldstep=0.01, finaltime=0.03):
1121            pass
1122
1123
1124        assert num.allclose(domain.recorded_min_timestep, 0.0210448446782)
1125        assert num.allclose(domain.recorded_max_timestep, 0.0210448446782)
1126
1127
1128        assert num.allclose(domain.quantities['xmomentum'].vertex_values[:4,0],
1129                            [0.00090268, 0.03684904, 0.00084545, 0.03686323])
1130
1131        assert num.allclose(domain.quantities['ymomentum'].vertex_values[:4,0],
1132                            [ -1.97318203e-04, 6.10268320e-04, -6.18053986e-05, 6.14082609e-04])
1133
1134        os.remove(domain.get_name() + '.sww')
1135
1136    def test_flatbed_second_order_distribute(self):
1137        #Use real data from anuga.abstract_2d_finite_volumes 2
1138        #painfully setup and extracted.
1139
1140        from mesh_factory import rectangular
1141
1142        # Create basic mesh
1143        N = 8
1144        points, vertices, boundary = rectangular(N, N)
1145
1146        # Create shallow water domain
1147        domain = Domain(points, vertices, boundary)
1148
1149        domain.set_store_vertices_uniquely(True)
1150        domain.set_default_order(2)
1151
1152        # Boundary conditions
1153        Br = Reflective_boundary(domain)
1154        Bd = Dirichlet_boundary([0.2, 0., 0.])
1155
1156        domain.set_boundary({'left': Bd, 'right': Br, 'top': Br, 'bottom': Br})
1157        domain.check_integrity()
1158
1159        for V in [False, True]:
1160            if V:
1161                # Set centroids as if system had been evolved
1162                L = num.zeros(2*N*N, num.float)
1163                L[:32] = [7.21205592e-003, 5.35214298e-002, 1.00910824e-002,
1164                          5.35439433e-002, 1.00910824e-002, 5.35439433e-002,
1165                          1.00910824e-002, 5.35439433e-002, 1.00910824e-002,
1166                          5.35439433e-002, 1.00910824e-002, 5.35439433e-002,
1167                          1.00910824e-002, 5.35393928e-002, 1.02344264e-002,
1168                          5.59605058e-002, 0.00000000e+000, 3.31027800e-004,
1169                          0.00000000e+000, 4.37962142e-005, 0.00000000e+000,
1170                          4.37962142e-005, 0.00000000e+000, 4.37962142e-005,
1171                          0.00000000e+000, 4.37962142e-005, 0.00000000e+000,
1172                          4.37962142e-005, 0.00000000e+000, 4.37962142e-005,
1173                          0.00000000e+000, 5.57305948e-005]
1174
1175                X = num.zeros(2*N*N, num.float)
1176                X[:32] = [6.48351607e-003, 3.68571894e-002, 8.50733285e-003,
1177                          3.68731327e-002, 8.50733285e-003, 3.68731327e-002,
1178                          8.50733285e-003, 3.68731327e-002, 8.50733285e-003,
1179                          3.68731327e-002, 8.50733285e-003, 3.68731327e-002,
1180                          8.50733285e-003, 3.68693861e-002, 8.65220973e-003,
1181                          3.85055387e-002, 0.00000000e+000, 2.86060840e-004,
1182                          0.00000000e+000, 3.58905503e-005, 0.00000000e+000,
1183                          3.58905503e-005, 0.00000000e+000, 3.58905503e-005,
1184                          0.00000000e+000, 3.58905503e-005, 0.00000000e+000,
1185                          3.58905503e-005, 0.00000000e+000, 3.58905503e-005,
1186                          0.00000000e+000, 4.57662812e-005]
1187
1188                Y = num.zeros(2*N*N, num.float)
1189                Y[:32] = [-1.39463104e-003, 6.15600298e-004, -6.03637382e-004,
1190                          6.18272251e-004, -6.03637382e-004, 6.18272251e-004,
1191                          -6.03637382e-004, 6.18272251e-004, -6.03637382e-004,
1192                          6.18272251e-004, -6.03637382e-004, 6.18272251e-004,
1193                          -6.03637382e-004, 6.18599320e-004, -6.74622797e-004,
1194                          -1.48934756e-004, 0.00000000e+000, -5.35079969e-005,
1195                          0.00000000e+000, -2.57264987e-005, 0.00000000e+000,
1196                          -2.57264987e-005, 0.00000000e+000, -2.57264987e-005,
1197                          0.00000000e+000, -2.57264987e-005, 0.00000000e+000,
1198                          -2.57264987e-005, 0.00000000e+000, -2.57264987e-005,
1199                          0.00000000e+000, -2.57635178e-005]
1200
1201                domain.set_quantity('stage', L, location='centroids')
1202                domain.set_quantity('xmomentum', X, location='centroids')
1203                domain.set_quantity('ymomentum', Y, location='centroids')
1204
1205                domain.check_integrity()
1206            else:
1207                # Evolution
1208                for t in domain.evolve(yieldstep=0.01, finaltime=0.03):
1209                    pass
1210                assert num.allclose(domain.recorded_min_timestep, 0.0210448446782)
1211                assert num.allclose(domain.recorded_max_timestep, 0.0210448446782)
1212
1213            #print domain.quantities['stage'].centroid_values[:4]
1214            #print domain.quantities['xmomentum'].centroid_values[:4]
1215            #print domain.quantities['ymomentum'].centroid_values[:4]                       
1216               
1217            #Centroids were correct but not vertices.
1218            #Hence the check of distribute below.
1219
1220            if not V:
1221                assert num.allclose(domain.quantities['stage'].centroid_values[:4],
1222                                    [0.00725574, 0.05350737, 0.01008413, 0.0535293])           
1223                assert num.allclose(domain.quantities['xmomentum'].centroid_values[:4],
1224                                    [0.00654964, 0.03684904, 0.00852561, 0.03686323])
1225
1226                assert num.allclose(domain.quantities['ymomentum'].centroid_values[:4],
1227                                    [-0.00143169, 0.00061027, -0.00062325, 0.00061408])
1228
1229                                   
1230           
1231                assert num.allclose(domain.quantities['xmomentum'].centroid_values[17], 0.0)               
1232            else:
1233                assert num.allclose(domain.quantities['xmomentum'].\
1234                                        centroid_values[17],
1235                                    0.00028606084)
1236                return #FIXME - Bailout for V True
1237
1238            import copy
1239
1240            XX = copy.copy(domain.quantities['xmomentum'].centroid_values)
1241            assert num.allclose(domain.quantities['xmomentum'].centroid_values,
1242                                XX)
1243
1244            domain.distribute_to_vertices_and_edges()
1245
1246            assert num.allclose(domain.quantities['xmomentum'].centroid_values[17], 0.0)
1247
1248            assert num.allclose(domain.quantities['ymomentum'].vertex_values[:4,0],
1249                                [-1.97318203e-04, 6.10268320e-04, -6.18053986e-05, 6.14082609e-04])
1250
1251            assert num.allclose(domain.quantities['xmomentum'].vertex_values[:4,0],
1252                                [0.00090268,  0.03684904,  0.00084545,  0.03686323])
1253
1254
1255        os.remove(domain.get_name() + '.sww')
1256
1257
1258    def test_bedslope_problem_second_order_more_steps(self):
1259        """test_bedslope_problem_second_order_more_step
1260
1261        Test shallow water balanced finite volumes
1262        """
1263
1264        from mesh_factory import rectangular
1265
1266        # Create basic mesh
1267        points, vertices, boundary = rectangular(6, 6)
1268
1269        # Create shallow water domain
1270        domain = Domain(points, vertices, boundary)
1271
1272        domain.set_store_vertices_uniquely(True)
1273        domain.set_default_order(2)
1274
1275        # Bed-slope and friction at vertices (and interpolated elsewhere)
1276        def x_slope(x, y):
1277            return -x/3
1278
1279        domain.set_quantity('elevation', x_slope)
1280
1281        # Boundary conditions
1282        Br = Reflective_boundary(domain)
1283        domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
1284
1285        # Initial condition
1286        domain.set_quantity('stage', expression='elevation+0.05')
1287        domain.check_integrity()
1288
1289        assert num.allclose(domain.quantities['stage'].centroid_values,
1290                            [ 0.01296296,  0.03148148,  0.01296296,
1291                              0.03148148,  0.01296296,  0.03148148,
1292                              0.01296296,  0.03148148,  0.01296296,
1293                              0.03148148,  0.01296296,  0.03148148,
1294                             -0.04259259, -0.02407407, -0.04259259,
1295                             -0.02407407, -0.04259259, -0.02407407,
1296                             -0.04259259, -0.02407407, -0.04259259,
1297                             -0.02407407, -0.04259259, -0.02407407,
1298                             -0.09814815, -0.07962963, -0.09814815,
1299                             -0.07962963, -0.09814815, -0.07962963,
1300                             -0.09814815, -0.07962963, -0.09814815,
1301                             -0.07962963, -0.09814815, -0.07962963,
1302                             -0.1537037 , -0.13518519, -0.1537037,
1303                             -0.13518519, -0.1537037,  -0.13518519,
1304                             -0.1537037 , -0.13518519, -0.1537037,
1305                             -0.13518519, -0.1537037,  -0.13518519,
1306                             -0.20925926, -0.19074074, -0.20925926,
1307                             -0.19074074, -0.20925926, -0.19074074,
1308                             -0.20925926, -0.19074074, -0.20925926,
1309                             -0.19074074, -0.20925926, -0.19074074,
1310                             -0.26481481, -0.2462963,  -0.26481481,
1311                             -0.2462963,  -0.26481481, -0.2462963,
1312                             -0.26481481, -0.2462963,  -0.26481481,
1313                             -0.2462963,  -0.26481481, -0.2462963])
1314
1315        # Evolution
1316        for t in domain.evolve(yieldstep = 0.05, finaltime = 0.5):
1317            pass
1318
1319
1320       
1321        assert num.allclose(domain.quantities['stage'].centroid_values,
1322     [-0.02901283, -0.01619385, -0.03040423, -0.01564474, -0.02936756, -0.01507953,
1323      -0.02858108, -0.01491531, -0.02793549, -0.0147037,  -0.02792804, -0.014363,
1324      -0.07794301, -0.05951952, -0.07675098, -0.06182336, -0.07736607, -0.06079504,
1325      -0.07696764, -0.06039043, -0.07708793, -0.0601453,  -0.07669911, -0.06020719,
1326      -0.12223185, -0.10857309, -0.12286676, -0.10837591, -0.12386938, -0.10842744,
1327      -0.12363769, -0.10790002, -0.12304837, -0.10871278, -0.12543768, -0.10961026,
1328      -0.15664473, -0.14630207, -0.15838364, -0.14910271, -0.15804002, -0.15029627,
1329      -0.15829717, -0.1503869,  -0.15852604, -0.14971109, -0.15856346, -0.15205092,
1330      -0.20900931, -0.19658843, -0.20669607, -0.19558708, -0.20654186, -0.19492423,
1331      -0.20438765, -0.19492931, -0.20644142, -0.19423147, -0.20237449, -0.19198454,
1332      -0.13699658, -0.14209126, -0.13600697, -0.14334968, -0.1347657,  -0.14224247,
1333      -0.13442376, -0.14136926, -0.13501004, -0.14339389, -0.13479263, -0.14304073])
1334
1335       
1336
1337        assert num.allclose(domain.quantities['xmomentum'].centroid_values,
1338[0.0053808980620685537, 0.0018623678353073116, 0.0047019894631921888, 0.002140100234595967, 0.0050544722142569594, 0.002413322128452734, 0.005489085956540414, 0.0025243661256655553, 0.0060037830243309179, 0.0026391226895162612, 0.0057883734209939882, 0.0026956732546692926, 0.015940649159142187, 0.013221882156398579, 0.015888047512240901, 0.01158701848044226, 0.01546380192164616, 0.012392900653821705, 0.01581341145685258, 0.012655591326883386, 0.015589268030001307, 0.012882167778653425, 0.016378187240908375, 0.01285389097473179, 0.033554812180096157, 0.024722935925976689, 0.031394874407697747, 0.025522168163393786, 0.030489844987314309, 0.025816648228773609, 0.03093376828330165, 0.026526382431385956, 0.031779480374536206, 0.025244211480815279, 0.02777110473340063, 0.0235612830668613, 0.072805732021172603, 0.057771452382474019, 0.070995840015208006, 0.052063135807599033, 0.071682005137710475, 0.050427261692222392, 0.071198588672075042, 0.050412342621995315, 0.071083783829814381, 0.05152744137515769, 0.071295902924896015, 0.046793561462568523, 0.08896512801938701, 0.073620532919998594, 0.09050528117516124, 0.076886136136002675, 0.089887310709307736, 0.076834171935627513, 0.090202740570079903, 0.077601818401490483, 0.091197277460468809, 0.07791128140944184, 0.091598726283965259, 0.077544779639518807, 0.020200091779226687, 0.058267129156556331, 0.026187571427719752, 0.057931516481767524, 0.025402693883943676, 0.056813327712684755, 0.024916381753277369, 0.056341859717484941, 0.024736292276481896, 0.05716071583083765, 0.026274297871292408, 0.07511805936685842])
1339
1340
1341        assert num.allclose(domain.quantities['ymomentum'].centroid_values,
1342[0.00033542456924774407, -0.00020012758197726979, -0.00033105661215122639, -0.00012291548928474255, -8.6598990751306055e-05, -5.3679813150316306e-05, 2.7774382762402742e-05, -9.3519331403178185e-05, -0.00019125171262773737, -0.00022905710988083868, -0.00038249823793758749, -5.2279522092562622e-05, -0.00032248606612494252, 8.870124179963529e-05, -0.00037840179563069224, -0.00038359452748757826, -0.0003248974462485901, -0.00012753304045182617, -0.0001591017972094672, 1.8279217568228501e-05, -6.1000546782769864e-05, -1.331229915505809e-05, -6.253286589681782e-05, -0.0002488614999307059, 0.0011988796270581575, 0.00017258171683877814, 3.4021626634331032e-05, -0.00036969454859050733, -0.00033874782370460032, -0.00031089795347570575, -0.0002999150988746842, -0.00037403606927378225, -0.00048310389377791813, -0.00010570764663927565, 0.00079563172917226932, 0.00078476438788764808, 0.0018347928776038006, 0.00072944213736041619, 0.0007060579464053257, 3.394251412720066e-05, 0.00053023191377378964, -0.00038114184186005149, 0.00037324507881442268, -0.00029480847037389904, 0.00052318790374037533, -0.00065867970688702289, -0.00047558178231081052, 0.00038297067147786805, -0.00010701572465258302, 0.0016609093113296653, 0.00072099989450924115, 0.00083723255250903368, 0.0004402978878512923, 0.00071527026447847206, 0.00061764112501907131, 0.0009682410892776223, 8.8194340495455049e-05, 0.00032386823106466095, -0.0014131220131695192, 0.00034752669349133577, -0.0017518907583968107, -0.0032179499168180402, -0.0030608351073009841, -0.0019003818511794108, -0.0019268160268303125, -0.0016355558234909565, -0.0018559374675595419, -0.0012213557447637634, -0.00195465742442113, 0.00016169839310254064, 0.0031989528671111625, -0.0018028271632022301])
1343
1344        os.remove(domain.get_name() + '.sww')
1345
1346
1347    def test_temp_play(self):
1348        from mesh_factory import rectangular
1349
1350        # Create basic mesh
1351        points, vertices, boundary = rectangular(5, 5)
1352
1353        # Create shallow water domain
1354        domain = Domain(points, vertices, boundary)
1355        domain.smooth = False
1356        domain.default_order = 2
1357        domain.beta_w = 0.9
1358        domain.beta_w_dry = 0.9
1359        domain.beta_uh = 0.9
1360        domain.beta_uh_dry = 0.9
1361        domain.beta_vh = 0.9
1362        domain.beta_vh_dry = 0.9
1363
1364        # FIXME (Ole): Need tests where these two are commented out
1365        domain.H0 = 0                         # Backwards compatibility (6/2/7)
1366        domain.tight_slope_limiters = False   # Backwards compatibility (14/4/7)
1367        domain.use_centroid_velocities = False # Backwards compatibility (7/5/8)
1368        domain.use_edge_limiter = False       # Backwards compatibility (9/5/8)
1369
1370
1371        # Bed-slope and friction at vertices (and interpolated elsewhere)
1372        def x_slope(x, y):
1373            return -x/3
1374
1375        domain.set_quantity('elevation', x_slope)
1376
1377        # Boundary conditions
1378        Br = Reflective_boundary(domain)
1379        domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
1380
1381        # Initial condition
1382        domain.set_quantity('stage', expression='elevation+0.05')
1383        domain.check_integrity()
1384
1385        # Evolution
1386        for t in domain.evolve(yieldstep=0.05, finaltime=0.1):
1387            pass
1388
1389        assert num.allclose(domain.quantities['stage'].centroid_values[:4],
1390                            [0.00206836, 0.01296714, 0.00363415, 0.01438924])
1391        assert num.allclose(domain.quantities['xmomentum'].centroid_values[:4],
1392                            [0.01360154, 0.00671133, 0.01264578, 0.00648503])
1393        assert num.allclose(domain.quantities['ymomentum'].centroid_values[:4],
1394                            [-1.19201077e-003, -7.23647546e-004,
1395                             -6.39083123e-005, 6.29815168e-005])
1396
1397        os.remove(domain.get_name() + '.sww')
1398
1399    def test_complex_bed(self):
1400        # No friction is tested here
1401
1402        from mesh_factory import rectangular
1403
1404        N = 12
1405        points, vertices, boundary = rectangular(N, N/2, len1=1.2, len2=0.6,
1406                                                 origin=(-0.07, 0))
1407
1408
1409        domain = Domain(points, vertices, boundary)
1410        domain.smooth = False
1411        domain.default_order = 2
1412
1413        inflow_stage = 0.1
1414        Z = Weir(inflow_stage)
1415        domain.set_quantity('elevation', Z)
1416
1417        Br = Reflective_boundary(domain)
1418        Bd = Dirichlet_boundary([inflow_stage, 0.0, 0.0])
1419        domain.set_boundary({'left': Bd, 'right': Br, 'bottom': Br, 'top': Br})
1420
1421        domain.set_quantity('stage', expression='elevation')
1422
1423        for t in domain.evolve(yieldstep=0.02, finaltime=0.2):
1424            pass
1425
1426        #FIXME: These numbers were from version before 25/10
1427        #assert allclose(domain.quantities['stage'].centroid_values,
1428# [3.95822638e-002,  5.61022588e-002,  4.66437868e-002,  5.73081011e-002,
1429#  4.72394613e-002,  5.74684939e-002,  4.74309483e-002,  5.77458084e-002,
1430#  4.80628177e-002,  5.85656225e-002,  4.90498542e-002,  6.02609831e-002,
1431#  1.18470315e-002,  1.75136443e-002,  1.18035266e-002,  2.15565695e-002,
1432#  1.31620268e-002,  2.14351640e-002,  1.32351076e-002,  2.15450687e-002,
1433#  1.36414028e-002,  2.24274619e-002,  1.51689511e-002,  2.21789655e-002,
1434# -7.54337535e-003, -6.86362021e-004, -7.74146760e-003, -1.83756530e-003,
1435# -8.16773628e-003, -4.49916813e-004, -8.08202599e-003, -3.91118720e-004,
1436# -8.10292716e-003, -3.88584984e-004, -7.35226124e-003,  2.73985295e-004,
1437#  1.86166683e-001,  8.74070369e-002,  1.86166712e-001,  8.74035875e-002,
1438#  6.11666935e-002, -3.76173225e-002, -6.38333276e-002, -3.76147365e-002,
1439#  6.11666725e-002,  8.73846774e-002,  1.86166697e-001,  8.74171550e-002,
1440# -4.83333333e-002,  1.18333333e-001, -4.83333333e-002,  1.18333333e-001,
1441# -4.83333333e-002, -6.66666667e-003, -1.73333333e-001, -1.31666667e-001,
1442# -1.73333333e-001, -6.66666667e-003, -4.83333333e-002,  1.18333333e-001,
1443# -2.48333333e-001, -2.31666667e-001, -2.48333333e-001, -2.31666667e-001,
1444# -2.48333333e-001, -2.31666667e-001, -2.48333333e-001, -2.31666667e-001,
1445# -2.48333333e-001, -2.31666667e-001, -2.48333333e-001, -2.31666667e-001,
1446# -4.65000000e-001, -3.65000000e-001, -4.65000000e-001, -3.65000000e-001,
1447# -4.65000000e-001, -3.65000000e-001, -4.65000000e-001, -3.65000000e-001,
1448# -4.65000000e-001, -3.65000000e-001, -4.65000000e-001, -3.65000000e-001,
1449# -5.98333333e-001, -5.81666667e-001, -5.98333333e-001, -5.81666667e-001,
1450# -5.98333333e-001, -5.81666667e-001, -5.98333333e-001, -5.81666667e-001,
1451# -5.98333333e-001, -5.81666667e-001, -5.98333333e-001, -5.81666667e-001,
1452# -6.48333333e-001, -6.31666667e-001, -6.48333333e-001, -6.31666667e-001,
1453# -6.48333333e-001, -6.31666667e-001, -6.48333333e-001, -6.31666667e-001,
1454# -6.48333333e-001, -6.31666667e-001, -6.48333333e-001, -6.31666667e-001,
1455# -5.31666667e-001, -5.98333333e-001, -5.31666667e-001, -5.98333333e-001,
1456# -5.31666667e-001, -5.98333333e-001, -5.31666667e-001, -5.98333333e-001,
1457# -5.31666667e-001, -5.98333333e-001, -5.31666667e-001, -5.98333333e-001,
1458# -4.98333333e-001, -4.81666667e-001, -4.98333333e-001, -4.81666667e-001,
1459# -4.98333333e-001, -4.81666667e-001, -4.98333333e-001, -4.81666667e-001,
1460# -4.98333333e-001, -4.81666667e-001, -4.98333333e-001, -4.81666667e-001,
1461# -5.48333333e-001, -5.31666667e-001, -5.48333333e-001, -5.31666667e-001,
1462# -5.48333333e-001, -5.31666667e-001, -5.48333333e-001, -5.31666667e-001,
1463# -5.48333333e-001, -5.31666667e-001, -5.48333333e-001, -5.31666667e-001])
1464
1465        os.remove(domain.get_name() + '.sww')
1466
1467
1468    def test_tight_slope_limiters(self):
1469        """Test that new slope limiters (Feb 2007) don't induce extremely
1470        small timesteps. This test actually reveals the problem as it
1471        was in March-April 2007
1472        """
1473        import time, os
1474        from Scientific.IO.NetCDF import NetCDFFile
1475        from data_manager import extent_sww
1476        from mesh_factory import rectangular
1477
1478        # Create basic mesh
1479        points, vertices, boundary = rectangular(2, 2)
1480
1481        # Create shallow water domain
1482        domain = Domain(points, vertices, boundary)
1483        domain.default_order = 2
1484
1485        # This will pass
1486        #domain.tight_slope_limiters = 1
1487        #domain.H0 = 0.01
1488
1489        # This will fail
1490        #domain.tight_slope_limiters = 1
1491        #domain.H0 = 0.001
1492
1493        # This will pass provided C extension implements limiting of
1494        # momentum in _compute_speeds
1495        domain.tight_slope_limiters = 1
1496        domain.H0 = 0.001
1497        domain.protect_against_isolated_degenerate_timesteps = True
1498
1499        # Set some field values
1500        domain.set_quantity('elevation', lambda x,y: -x)
1501        domain.set_quantity('friction', 0.03)
1502
1503        # Boundary conditions
1504        B = Transmissive_boundary(domain)
1505        domain.set_boundary({'left': B, 'right': B, 'top': B, 'bottom': B})
1506
1507        # Initial condition - with jumps
1508        bed = domain.quantities['elevation'].vertex_values
1509        stage = num.zeros(bed.shape, num.float)
1510
1511        h = 0.3
1512        for i in range(stage.shape[0]):
1513            if i % 2 == 0:
1514                stage[i,:] = bed[i,:] + h
1515            else:
1516                stage[i,:] = bed[i,:]
1517
1518        domain.set_quantity('stage', stage)
1519
1520        domain.distribute_to_vertices_and_edges()
1521
1522        domain.set_name('tight_limiters')
1523        domain.smooth = True
1524        domain.reduction = mean
1525        domain.set_datadir('.')
1526        domain.smooth = False
1527        domain.store = True
1528
1529        # Evolution
1530        for t in domain.evolve(yieldstep=0.1, finaltime=0.3):
1531            #domain.write_time(track_speeds=True)
1532            stage = domain.quantities['stage'].vertex_values
1533
1534            # Get NetCDF
1535            fid = NetCDFFile(domain.writer.filename, netcdf_mode_r)
1536            stage_file = fid.variables['stage']
1537
1538            fid.close()
1539
1540        os.remove(domain.writer.filename)
1541
1542
1543
1544    def test_pmesh2Domain(self):
1545         import os
1546         import tempfile
1547
1548         fileName = tempfile.mktemp(".tsh")
1549         file = open(fileName, "w")
1550         file.write("4 3 # <vertex #> <x> <y> [attributes]\n \
15510 0.0 0.0 0.0 0.0 0.01 \n \
15521 1.0 0.0 10.0 10.0 0.02  \n \
15532 0.0 1.0 0.0 10.0 0.03  \n \
15543 0.5 0.25 8.0 12.0 0.04  \n \
1555# Vert att title  \n \
1556elevation  \n \
1557stage  \n \
1558friction  \n \
15592 # <triangle #> [<vertex #>] [<neigbouring triangle #>]  \n\
15600 0 3 2 -1  -1  1 dsg\n\
15611 0 1 3 -1  0 -1   ole nielsen\n\
15624 # <segment #> <vertex #>  <vertex #> [boundary tag] \n\
15630 1 0 2 \n\
15641 0 2 3 \n\
15652 2 3 \n\
15663 3 1 1 \n\
15673 0 # <x> <y> [attributes] ...Mesh Vertices... \n \
15680 216.0 -86.0 \n \
15691 160.0 -167.0 \n \
15702 114.0 -91.0 \n \
15713 # <vertex #>  <vertex #> [boundary tag] ...Mesh Segments... \n \
15720 0 1 0 \n \
15731 1 2 0 \n \
15742 2 0 0 \n \
15750 # <x> <y> ...Mesh Holes... \n \
15760 # <x> <y> <attribute>...Mesh Regions... \n \
15770 # <x> <y> <attribute>...Mesh Regions, area... \n\
1578#Geo reference \n \
157956 \n \
1580140 \n \
1581120 \n")
1582         file.close()
1583
1584         tags = {}
1585         b1 =  Dirichlet_boundary(conserved_quantities = num.array([0.0]))
1586         b2 =  Dirichlet_boundary(conserved_quantities = num.array([1.0]))
1587         b3 =  Dirichlet_boundary(conserved_quantities = num.array([2.0]))
1588         tags["1"] = b1
1589         tags["2"] = b2
1590         tags["3"] = b3
1591
1592         domain = Domain(mesh_filename=fileName)
1593                         # verbose=True, use_cache=True)
1594
1595         ## check the quantities
1596         answer = [[0., 8., 0.],
1597                   [0., 10., 8.]]
1598         assert num.allclose(domain.quantities['elevation'].vertex_values,
1599                             answer)
1600
1601         answer = [[0., 12., 10.],
1602                   [0., 10., 12.]]
1603         assert num.allclose(domain.quantities['stage'].vertex_values,
1604                             answer)
1605
1606         answer = [[0.01, 0.04, 0.03],
1607                   [0.01, 0.02, 0.04]]
1608         assert num.allclose(domain.quantities['friction'].vertex_values,
1609                             answer)
1610
1611         tagged_elements = domain.get_tagged_elements()
1612         assert num.allclose(tagged_elements['dsg'][0], 0)
1613         assert num.allclose(tagged_elements['ole nielsen'][0], 1)
1614
1615         msg = "test_tags_to_boundaries failed. Single boundary wasn't added."
1616         self.failUnless( domain.boundary[(1, 0)]  == '1', msg)
1617         self.failUnless( domain.boundary[(1, 2)]  == '2', msg)
1618         self.failUnless( domain.boundary[(0, 1)]  == '3', msg)
1619         self.failUnless( domain.boundary[(0, 0)]  == 'exterior', msg)
1620         msg = "test_pmesh2Domain Too many boundaries"
1621         self.failUnless( len(domain.boundary)  == 4, msg)
1622
1623         # FIXME change to use get_xllcorner
1624         msg = 'Bad geo-reference'
1625         self.failUnless(domain.geo_reference.xllcorner  == 140.0, msg)
1626
1627         domain = Domain(fileName)
1628
1629         answer = [[0., 8., 0.],
1630                   [0., 10., 8.]]
1631         assert num.allclose(domain.quantities['elevation'].vertex_values,
1632                             answer)
1633
1634         answer = [[0., 12., 10.],
1635                   [0., 10., 12.]]
1636         assert num.allclose(domain.quantities['stage'].vertex_values,
1637                             answer)
1638
1639         answer = [[0.01, 0.04, 0.03],
1640                   [0.01, 0.02, 0.04]]
1641         assert num.allclose(domain.quantities['friction'].vertex_values,
1642                             answer)
1643
1644         tagged_elements = domain.get_tagged_elements()
1645         assert num.allclose(tagged_elements['dsg'][0], 0)
1646         assert num.allclose(tagged_elements['ole nielsen'][0], 1)
1647
1648         msg = "test_tags_to_boundaries failed. Single boundary wasn't added."
1649         self.failUnless(domain.boundary[(1, 0)]  == '1', msg)
1650         self.failUnless(domain.boundary[(1, 2)]  == '2', msg)
1651         self.failUnless(domain.boundary[(0, 1)]  == '3', msg)
1652         self.failUnless(domain.boundary[(0, 0)]  == 'exterior', msg)
1653         msg = "test_pmesh2Domain Too many boundaries"
1654         self.failUnless(len(domain.boundary)  == 4, msg)
1655
1656         # FIXME change to use get_xllcorner
1657         msg = 'Bad geo_reference'
1658         self.failUnless(domain.geo_reference.xllcorner  == 140.0, msg)
1659
1660         os.remove(fileName)
1661
1662    def test_get_lone_vertices(self):
1663        a = [0.0, 0.0]
1664        b = [0.0, 2.0]
1665        c = [2.0, 0.0]
1666        d = [0.0, 4.0]
1667        e = [2.0, 2.0]
1668        f = [4.0, 0.0]
1669
1670        points = [a, b, c, d, e, f]
1671        #             bac,     bce,     ecf,     dbe
1672        vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4] ]
1673        boundary = {(0, 0): 'Third',
1674                    (0, 2): 'First',
1675                    (2, 0): 'Second',
1676                    (2, 1): 'Second',
1677                    (3, 1): 'Second',
1678                    (3, 2): 'Third'}
1679
1680        domain = Domain(points, vertices, boundary)
1681        domain.get_lone_vertices()
1682
1683    def test_fitting_using_shallow_water_domain(self):
1684        #Mesh in zone 56 (absolute coords)
1685
1686        x0 = 314036.58727982
1687        y0 = 6224951.2960092
1688
1689        a = [x0+0.0, y0+0.0]
1690        b = [x0+0.0, y0+2.0]
1691        c = [x0+2.0, y0+0.0]
1692        d = [x0+0.0, y0+4.0]
1693        e = [x0+2.0, y0+2.0]
1694        f = [x0+4.0, y0+0.0]
1695
1696        points = [a, b, c, d, e, f]
1697
1698        #             bac,     bce,     ecf,     dbe
1699        elements = [[1,0,2], [1,2,4], [4,2,5], [3,1,4] ]
1700
1701        # absolute going in ..
1702        mesh4 = Domain(points, elements, geo_reference=Geo_reference(56, 0, 0))
1703        mesh4.check_integrity()
1704        quantity = Quantity(mesh4)
1705
1706        # Get (enough) datapoints (relative to georef)
1707        data_points_rel = [[ 0.66666667, 0.66666667],
1708                           [ 1.33333333, 1.33333333],
1709                           [ 2.66666667, 0.66666667],
1710                           [ 0.66666667, 2.66666667],
1711                           [ 0.0,        1.0],
1712                           [ 0.0,        3.0],
1713                           [ 1.0,        0.0],
1714                           [ 1.0,        1.0],
1715                           [ 1.0,        2.0],
1716                           [ 1.0,        3.0],
1717                           [ 2.0,        1.0],
1718                           [ 3.0,        0.0],
1719                           [ 3.0,        1.0]]
1720
1721        data_geo_spatial = Geospatial_data(data_points_rel,
1722                                           geo_reference=Geo_reference(56,
1723                                                                       x0,
1724                                                                       y0))
1725        data_points_absolute = data_geo_spatial.get_data_points(absolute=True)
1726        attributes = linear_function(data_points_absolute)
1727        att = 'spam_and_eggs'
1728
1729        # Create .txt file
1730        ptsfile = tempfile.mktemp(".txt")
1731        file = open(ptsfile, "w")
1732        file.write(" x,y," + att + " \n")
1733        for data_point, attribute in map(None, data_points_absolute, attributes):
1734            row = (str(data_point[0]) + ',' +
1735                   str(data_point[1]) + ',' +
1736                   str(attribute))
1737            file.write(row + "\n")
1738        file.close()
1739
1740        # Check that values can be set from file
1741        quantity.set_values(filename=ptsfile, attribute_name=att, alpha=0)
1742        answer = linear_function(quantity.domain.get_vertex_coordinates())
1743
1744        assert num.allclose(quantity.vertex_values.flat, answer)
1745
1746        # Check that values can be set from file using default attribute
1747        quantity.set_values(filename = ptsfile, alpha = 0)
1748        assert num.allclose(quantity.vertex_values.flat, answer)
1749
1750        # Cleanup
1751        import os
1752        os.remove(ptsfile)
1753
1754    def test_fitting_example_that_crashed(self):
1755        """This unit test has been derived from a real world example
1756        (the Towradgi '98 rainstorm simulation).
1757
1758        It shows a condition where fitting as called from set_quantity crashes
1759        when ANUGA mesh is reused. The test passes in the case where a new mesh
1760        is created.
1761
1762        See ticket:314
1763        """
1764
1765        verbose = False
1766
1767        from anuga.shallow_water import Domain
1768        from anuga.pmesh.mesh_interface import create_mesh_from_regions
1769        from anuga.geospatial_data.geospatial_data import Geospatial_data
1770
1771       
1772        # Get path where this test is run
1773        path = get_pathname_from_package('anuga.shallow_water')       
1774       
1775       
1776        #----------------------------------------------------------------------
1777        # Create domain
1778        #--------------------------------------------------------------------
1779        W = 303400
1780        N = 6195800
1781        E = 308640
1782        S = 6193120
1783        bounding_polygon = [[W, S], [E, S], [E, N], [W, N]]
1784
1785        offending_regions = []
1786
1787        # From culvert 8
1788        offending_regions.append([[307611.43896231, 6193631.6894806],
1789                                  [307600.11394969, 6193608.2855474],
1790                                  [307597.41349586, 6193609.59227963],
1791                                  [307608.73850848, 6193632.99621282]])
1792        offending_regions.append([[307633.69143231, 6193620.9216536],
1793                                  [307622.36641969, 6193597.5177204],
1794                                  [307625.06687352, 6193596.21098818],
1795                                  [307636.39188614, 6193619.61492137]])
1796
1797        # From culvert 9
1798        offending_regions.append([[306326.69660524, 6194818.62900522],
1799                                  [306324.67939476, 6194804.37099478],
1800                                  [306323.75856492, 6194804.50127295],
1801                                  [306325.7757754, 6194818.7592834]])
1802        offending_regions.append([[306365.57160524, 6194813.12900522],
1803                                  [306363.55439476, 6194798.87099478],
1804                                  [306364.4752246, 6194798.7407166],
1805                                  [306366.49243508, 6194812.99872705]])
1806
1807        # From culvert 10
1808        offending_regions.append([[306955.071019428608, 6194465.704096679576],
1809                                  [306951.616980571358, 6194457.295903320424],
1810                                  [306950.044491164153, 6194457.941873183474],
1811                                  [306953.498530021403, 6194466.350066542625]])
1812        offending_regions.append([[307002.540019428649, 6194446.204096679576],
1813                                  [306999.085980571399, 6194437.795903320424],
1814                                  [307000.658469978604, 6194437.149933457375],
1815                                  [307004.112508835853, 6194445.558126816526]])
1816
1817        interior_regions = []
1818        for polygon in offending_regions:
1819            interior_regions.append( [polygon, 100] ) 
1820
1821        meshname = os.path.join(path, 'offending_mesh.msh')
1822        create_mesh_from_regions(bounding_polygon,
1823                                 boundary_tags={'south': [0], 'east': [1],
1824                                                'north': [2], 'west': [3]},
1825                                 maximum_triangle_area=1000000,
1826                                 interior_regions=interior_regions,
1827                                 filename=meshname,
1828                                 use_cache=False,
1829                                 verbose=verbose)
1830
1831        domain = Domain(meshname, use_cache=False, verbose=verbose)
1832
1833        #----------------------------------------------------------------------
1834        # Fit data point to mesh
1835        #----------------------------------------------------------------------
1836
1837        points_file = os.path.join(path, 'offending_point.pts')
1838
1839        # Offending point
1840        G = Geospatial_data(data_points=[[306953.344, 6194461.5]],
1841                            attributes=[1])
1842        G.export_points_file(points_file)
1843       
1844        try:
1845            domain.set_quantity('elevation', filename=points_file,
1846                                use_cache=False, verbose=verbose, alpha=0.01)
1847        except RuntimeError, e:
1848            msg = 'Test failed: %s' % str(e)
1849            raise Exception, msg
1850            # clean up in case raise fails
1851            os.remove(meshname)
1852            os.remove(points_file)
1853        else:
1854            os.remove(meshname)
1855            os.remove(points_file)           
1856       
1857        #finally:
1858            # Cleanup regardless
1859            #FIXME(Ole): Finally does not work like this in python2.4
1860            #FIXME(Ole): Reinstate this when Python2.4 is out of the way
1861            #FIXME(Ole): Python 2.6 apparently introduces something called 'with'           
1862            #os.remove(meshname)
1863            #os.remove(points_file)
1864
1865
1866    def test_fitting_example_that_crashed_2(self):
1867        """test_fitting_example_that_crashed_2
1868       
1869        This unit test has been derived from a real world example
1870        (the JJKelly study, by Petar Milevski).
1871       
1872        It shows a condition where set_quantity crashes due to AtA
1873        not being built properly
1874       
1875        See ticket:314
1876        """
1877
1878        verbose = False       
1879
1880        from anuga.shallow_water import Domain
1881        from anuga.pmesh.mesh_interface import create_mesh_from_regions
1882        from anuga.geospatial_data import Geospatial_data
1883       
1884        # Get path where this test is run
1885        path = get_pathname_from_package('anuga.shallow_water')       
1886
1887        meshname = os.path.join(path, 'test_mesh.msh')
1888        W = 304180
1889        S = 6185270
1890        E = 307650
1891        N = 6189040
1892        maximum_triangle_area = 1000000
1893
1894        bounding_polygon = [[W, S], [E, S], [E, N], [W, N]]
1895
1896        create_mesh_from_regions(bounding_polygon,
1897                                 boundary_tags={'south': [0], 
1898                                                'east': [1], 
1899                                                'north': [2], 
1900                                                'west': [3]},
1901                                 maximum_triangle_area=maximum_triangle_area,
1902                                 filename=meshname,
1903                                 use_cache=False,
1904                                 verbose=verbose)
1905
1906        domain = Domain(meshname, use_cache=True, verbose=verbose)
1907       
1908        # Large test set revealed one problem
1909        points_file = os.path.join(path, 'test_points_large.csv')
1910        try:
1911            domain.set_quantity('elevation', filename=points_file,
1912                                use_cache=False, verbose=verbose)
1913        except AssertionError, e:
1914            msg = 'Test failed: %s' % str(e)
1915            raise Exception, msg
1916            # Cleanup in case this failed
1917            os.remove(meshname)
1918
1919        # Small test set revealed another problem
1920        points_file = os.path.join(path, 'test_points_small.csv')
1921        try:   
1922            domain.set_quantity('elevation', filename=points_file,
1923                                use_cache=False, verbose=verbose)                           
1924        except AssertionError, e:
1925            msg = 'Test failed: %s' % str(e)
1926            raise Exception, msg
1927            # Cleanup in case this failed
1928            os.remove(meshname)
1929        else:
1930            os.remove(meshname)
1931
1932
1933
1934
1935    def test_variable_elevation(self):           
1936        """test_variable_elevation
1937
1938        This will test that elevagtion van be stored in sww files
1939        as a time dependent quantity.
1940       
1941        It will also chck that storage of other quantities
1942        can be controlled this way.
1943        """
1944
1945        #---------------------------------------------------------------------
1946        # Import necessary modules
1947        #---------------------------------------------------------------------
1948        from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular_cross
1949        from anuga.shallow_water import Domain
1950        from anuga.shallow_water import Reflective_boundary
1951        from anuga.shallow_water import Dirichlet_boundary
1952        from anuga.shallow_water import Time_boundary
1953
1954        #---------------------------------------------------------------------
1955        # Setup computational domain
1956        #---------------------------------------------------------------------
1957        length = 8.
1958        width = 6.
1959        dx = dy = 1    # Resolution: Length of subdivisions on both axes
1960       
1961        inc = 0.05 # Elevation increment
1962
1963        points, vertices, boundary = rectangular_cross(int(length/dx), 
1964                                                       int(width/dy),
1965                                                       len1=length, 
1966                                                       len2=width)
1967        domain = Domain(points, vertices, boundary)
1968        domain.set_name('channel_variable_test')  # Output name
1969        domain.set_quantities_to_be_stored({'elevation': 2,
1970                                            'stage': 2})
1971
1972        #---------------------------------------------------------------------
1973        # Setup initial conditions
1974        #---------------------------------------------------------------------
1975
1976        def pole_increment(x,y):
1977            """This provides a small increment to a pole located mid stream
1978            For use with variable elevation data
1979            """
1980           
1981            z = 0.0*x
1982
1983            N = len(x)
1984            for i in range(N):
1985                # Pole
1986                if (x[i] - 4)**2 + (y[i] - 2)**2 < 1.0**2:
1987                    z[i] += inc
1988            return z
1989           
1990        domain.set_quantity('elevation', 0.0)    # Flat bed initially
1991        domain.set_quantity('friction', 0.01)    # Constant friction
1992        domain.set_quantity('stage', 0.0)        # Dry initial condition
1993
1994        #------------------------------------------------------------------
1995        # Setup boundary conditions
1996        #------------------------------------------------------------------
1997        Bi = Dirichlet_boundary([0.4, 0, 0])          # Inflow
1998        Br = Reflective_boundary(domain)              # Solid reflective wall
1999        Bo = Dirichlet_boundary([-5, 0, 0])           # Outflow
2000
2001        domain.set_boundary({'left': Bi, 'right': Bo, 'top': Br, 'bottom': Br})
2002
2003        #-------------------------------------------------------------------
2004        # Evolve system through time
2005        #-------------------------------------------------------------------
2006
2007        for t in domain.evolve(yieldstep=1, finaltime=3.0):
2008            #print domain.timestepping_statistics()
2009
2010            domain.add_quantity('elevation', pole_increment)
2011       
2012           
2013        # Check that quantities have been stored correctly   
2014        from Scientific.IO.NetCDF import NetCDFFile
2015        sww_file = domain.get_name() + '.sww'
2016        fid = NetCDFFile(sww_file)
2017
2018        x = fid.variables['x'][:]
2019        y = fid.variables['y'][:]
2020        stage = fid.variables['stage'][:]
2021        elevation = fid.variables['elevation'][:]
2022        fid.close()
2023
2024        os.remove(sww_file)
2025       
2026                   
2027        assert len(stage.shape) == 2
2028        assert len(elevation.shape) == 2       
2029       
2030        M, N = stage.shape
2031               
2032        for i in range(M): 
2033            # For each timestep
2034            assert num.allclose(max(elevation[i,:]), i * inc) 
2035
2036
2037
2038#################################################################################
2039
2040if __name__ == "__main__":
2041    suite = unittest.makeSuite(Test_swb_basic, 'test')
2042    runner = unittest.TextTestRunner(verbosity=1)
2043    runner.run(suite)
Note: See TracBrowser for help on using the repository browser.