source: branches/numpy/anuga/shallow_water/test_shallow_water_domain.py @ 6517

Last change on this file since 6517 was 6517, checked in by rwilson, 15 years ago

Hand-merged recent changes in main trunk. Still work to be done in shallow_water.

File size: 239.2 KB
Line 
1#!/usr/bin/env python
2
3import unittest, os
4from math import pi, sqrt
5import tempfile
6
7from anuga.config import g, epsilon
8from anuga.config import netcdf_mode_r, netcdf_mode_w, netcdf_mode_a
9from anuga.utilities.numerical_tools import mean
10from anuga.utilities.polygon import is_inside_polygon
11from anuga.coordinate_transforms.geo_reference import Geo_reference
12from anuga.abstract_2d_finite_volumes.quantity import Quantity
13from anuga.geospatial_data.geospatial_data import Geospatial_data
14from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular_cross
15
16from shallow_water_domain import *
17
18import numpy as num
19
20# Get gateway to C implementation of flux function for direct testing
21from shallow_water_ext import flux_function_central as flux_function
22
23
24# For test_fitting_using_shallow_water_domain example
25def linear_function(point):
26    point = num.array(point)
27    return point[:,0] + point[:,1]
28
29
30class Weir:
31    """Set a bathymetry for weir with a hole and a downstream gutter
32
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
121class Weir_simple:
122    """Set a bathymetry for weir with a hole and a downstream gutter
123
124    x,y are assumed to be in the unit square
125    """
126
127    def __init__(self, stage):
128        self.inflow_stage = stage
129
130    def __call__(self, x, y):
131        N = len(x)
132        assert N == len(y)
133
134        z = num.zeros(N, num.float)
135        for i in range(N):
136            z[i] = -x[i]    # General slope
137
138            # Flat bit to the left
139            if x[i] < 0.3:
140                z[i] = -x[i]/10    # General slope
141
142            # Weir
143            if x[i] > 0.3 and x[i] < 0.4:
144                z[i] = -x[i] + 0.9
145
146            # Dip
147            if x[i] > 0.6 and x[i] < 0.9:
148                z[i] = -x[i] - 0.5    # -y[i]/5
149
150            # Hole in weir (slightly higher than inflow condition)
151            if x[i] > 0.3 and x[i] < 0.4 and y[i] > 0.2 and y[i] < 0.4:
152                z[i] = -x[i] + self.inflow_stage + 0.05
153
154        return z/2
155
156
157# Variable windfield implemented using functions
158def speed(t, x, y):
159    """Large speeds halfway between center and edges
160
161    Low speeds at center and edges
162    """
163
164    from math import exp, cos, pi
165
166    x = num.array(x)
167    y = num.array(y)
168
169    N = len(x)
170    s = 0 * x    # New array
171
172    for k in range(N):
173        r = num.sqrt(x[k]**2 + y[k]**2)
174        factor = exp(-(r-0.15)**2)
175        s[k] = 4000 * factor * (cos(t*2*pi/150) + 2)
176
177    return s
178
179def scalar_func(t, x, y):
180    """Function that returns a scalar.
181
182    Used to test error message when numeric array is expected
183    """
184
185    return 17.7
186
187def scalar_func_list(t, x, y):
188    """Function that returns a scalar.
189
190    Used to test error message when numeric array is expected
191    """
192
193    return [17.7]
194
195def angle(t, x, y):
196    """Rotating field
197    """
198    from math import atan, pi
199
200    x = num.array(x)
201    y = num.array(y)
202
203    N = len(x)
204    a = 0 * x    # New array
205
206    for k in range(N):
207        r = num.sqrt(x[k]**2 + y[k]**2)
208
209        angle = atan(y[k]/x[k])
210
211        if x[k] < 0:
212            angle += pi
213
214        # Take normal direction
215        angle -= pi/2
216
217        # Ensure positive radians
218        if angle < 0:
219            angle += 2*pi
220
221        a[k] = angle/pi*180
222
223    return a
224
225
226class Test_Shallow_Water(unittest.TestCase):
227    def setUp(self):
228        pass
229
230    def tearDown(self):
231        pass
232
233    def test_rotate(self):
234        normal = num.array([0.0, -1.0])
235
236        q = num.array([1.0, 2.0, 3.0])
237
238        r = rotate(q, normal, direction = 1)
239        assert r[0] == 1
240        assert r[1] == -3
241        assert r[2] == 2
242
243        w = rotate(r, normal, direction = -1)
244        assert num.allclose(w, q)
245
246        # Check error check
247        try:
248            rotate(r, num.array([1, 1, 1]))
249        except:
250            pass
251        else:
252            raise Exception, 'Should have raised an exception'
253
254    # Individual flux tests
255    def test_flux_zero_case(self):
256        ql = num.zeros(3, num.float)
257        qr = num.zeros(3, num.float)
258        normal = num.zeros(2, num.float)
259        edgeflux = num.zeros(3, num.float)
260        zl = zr = 0.
261        H0 = 0.0
262
263        max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux,
264                                  epsilon, g, H0)
265
266        assert num.allclose(edgeflux, [0, 0, 0])
267        assert max_speed == 0.
268
269    def test_flux_constants(self):
270        w = 2.0
271
272        normal = num.array([1., 0])
273        ql = num.array([w, 0, 0])
274        qr = num.array([w, 0, 0])
275        edgeflux = num.zeros(3, num.float)
276        zl = zr = 0.
277        h = w - (zl+zr)/2
278        H0 = 0.0
279
280        max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux,
281                                  epsilon, g, H0)
282
283        assert num.allclose(edgeflux, [0., 0.5*g*h**2, 0.])
284        assert max_speed == num.sqrt(g*h)
285
286    #def test_flux_slope(self):
287    #    #FIXME: TODO
288    #    w = 2.0
289    #
290    #    normal = array([1., 0])
291    #    ql = array([w, 0, 0])
292    #    qr = array([w, 0, 0])
293    #    zl = zr = 0.
294    #    h = w - (zl+zr)/2
295    #
296    #    flux, max_speed = flux_function(normal, ql, qr, zl, zr)
297    #
298    #    assert allclose(flux, [0., 0.5*g*h**2, 0.])
299    #    assert max_speed == sqrt(g*h)
300
301    def test_flux1(self):
302        # Use data from previous version of abstract_2d_finite_volumes
303        normal = num.array([1., 0])
304        ql = num.array([-0.2, 2, 3])
305        qr = num.array([-0.2, 2, 3])
306        zl = zr = -0.5
307        edgeflux = num.zeros(3, num.float)
308
309        H0 = 0.0
310
311        max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux,
312                                  epsilon, g, H0)
313
314        assert num.allclose(edgeflux, [2., 13.77433333, 20.])
315        assert num.allclose(max_speed, 8.38130948661)
316
317    def test_flux2(self):
318        # Use data from previous version of abstract_2d_finite_volumes
319        normal = num.array([0., -1.])
320        ql = num.array([-0.075, 2, 3])
321        qr = num.array([-0.075, 2, 3])
322        zl = zr = -0.375
323
324        edgeflux = num.zeros(3, num.float)
325        H0 = 0.0
326        max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux,
327                                  epsilon, g, H0)
328
329        assert num.allclose(edgeflux, [-3., -20.0, -30.441])
330        assert num.allclose(max_speed, 11.7146428199)
331
332    def test_flux3(self):
333        # Use data from previous version of abstract_2d_finite_volumes
334        normal = num.array([-sqrt(2)/2, sqrt(2)/2])
335        ql = num.array([-0.075, 2, 3])
336        qr = num.array([-0.075, 2, 3])
337        zl = zr = -0.375
338
339        edgeflux = num.zeros(3, num.float)
340        H0 = 0.0
341        max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux,
342                                  epsilon, g, H0)
343
344        assert num.allclose(edgeflux, [sqrt(2)/2, 4.40221112, 7.3829019])
345        assert num.allclose(max_speed, 4.0716654239)
346
347    def test_flux4(self):
348        # Use data from previous version of abstract_2d_finite_volumes
349        normal = num.array([-sqrt(2)/2, sqrt(2)/2])
350        ql = num.array([-0.34319278, 0.10254161, 0.07273855])
351        qr = num.array([-0.30683287, 0.1071986, 0.05930515])
352        zl = zr = -0.375
353
354        edgeflux = num.zeros(3, num.float)
355        H0 = 0.0
356        max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux,
357                                  epsilon, g, H0)
358
359        assert num.allclose(edgeflux, [-0.04072676, -0.07096636, -0.01604364])
360        assert num.allclose(max_speed, 1.31414103233)
361
362    def test_flux_computation(self):
363        """test flux calculation (actual C implementation)
364
365        This one tests the constant case where only the pressure term
366        contributes to each edge and cancels out once the total flux has
367        been summed up.
368        """
369
370        a = [0.0, 0.0]
371        b = [0.0, 2.0]
372        c = [2.0, 0.0]
373        d = [0.0, 4.0]
374        e = [2.0, 2.0]
375        f = [4.0, 0.0]
376
377        points = [a, b, c, d, e, f]
378        #              bac,     bce,     ecf,     dbe
379        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
380
381        domain = Domain(points, vertices)
382        domain.check_integrity()
383
384        # The constant case
385        domain.set_quantity('elevation', -1)
386        domain.set_quantity('stage', 1)
387
388        domain.compute_fluxes()
389        # Central triangle
390        assert num.allclose(domain.get_quantity('stage').explicit_update[1], 0)
391
392        # The more general case
393        def surface(x, y):
394            return -x/2
395
396        domain.set_quantity('elevation', -10)
397        domain.set_quantity('stage', surface)
398        domain.set_quantity('xmomentum', 1)
399
400        domain.compute_fluxes()
401
402        #print domain.get_quantity('stage').explicit_update
403        # FIXME (Ole): TODO the general case
404        #assert allclose(domain.get_quantity('stage').explicit_update[1], ...??)
405
406    def test_sw_domain_simple(self):
407        a = [0.0, 0.0]
408        b = [0.0, 2.0]
409        c = [2.0, 0.0]
410        d = [0.0, 4.0]
411        e = [2.0, 2.0]
412        f = [4.0, 0.0]
413
414        points = [a, b, c, d, e, f]
415        #             bac,     bce,     ecf,     dbe
416        vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
417
418        #from anuga.abstract_2d_finite_volumes.domain import Domain as Generic_domain
419        #msg = 'The class %s is not a subclass of the generic domain class %s'\
420        #      %(DomainClass, Domain)
421        #assert issubclass(DomainClass, Domain), msg
422
423        domain = Domain(points, vertices)
424        domain.check_integrity()
425
426        for name in ['stage', 'xmomentum', 'ymomentum',
427                     'elevation', 'friction']:
428            assert domain.quantities.has_key(name)
429
430        assert num.alltrue(domain.get_conserved_quantities(0, edge=1) == 0.)
431
432    def test_boundary_conditions(self):
433        a = [0.0, 0.0]
434        b = [0.0, 2.0]
435        c = [2.0, 0.0]
436        d = [0.0, 4.0]
437        e = [2.0, 2.0]
438        f = [4.0, 0.0]
439
440        points = [a, b, c, d, e, f]
441        #             bac,     bce,     ecf,     dbe
442        vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
443        boundary = {(0, 0): 'Third',
444                    (0, 2): 'First',
445                    (2, 0): 'Second',
446                    (2, 1): 'Second',
447                    (3, 1): 'Second',
448                    (3, 2): 'Third'}
449
450        domain = Domain(points, vertices, boundary)
451        domain.check_integrity()
452
453        domain.set_quantity('stage', [[1,2,3], [5,5,5], [0,0,9], [-6,3,3]])
454
455        domain.set_quantity('xmomentum', [[1,1,1], [2,2,2], [3,3,3], [4,4,4]])
456
457        domain.set_quantity('ymomentum', [[10,10,10], [20,20,20],
458                                          [30,30,30], [40,40,40]])
459
460        D = Dirichlet_boundary([5, 2, 1])
461        T = Transmissive_boundary(domain)
462        R = Reflective_boundary(domain)
463        domain.set_boundary({'First': D, 'Second': T, 'Third': R})
464
465        domain.update_boundary()
466
467        # Stage
468        assert domain.quantities['stage'].boundary_values[0] == 2.5
469        # Reflective (2.5)
470        assert (domain.quantities['stage'].boundary_values[0] ==
471                domain.get_conserved_quantities(0, edge=0)[0])
472        # Dirichlet
473        assert domain.quantities['stage'].boundary_values[1] == 5.
474        # Transmissive (4.5)
475        assert (domain.quantities['stage'].boundary_values[2] ==
476                domain.get_conserved_quantities(2, edge=0)[0])
477        # Transmissive (4.5)
478        assert (domain.quantities['stage'].boundary_values[3] ==
479                domain.get_conserved_quantities(2, edge=1)[0])
480        # Transmissive (-1.5)
481        assert (domain.quantities['stage'].boundary_values[4] ==
482                domain.get_conserved_quantities(3, edge=1)[0])
483        # Reflective (-1.5)
484        assert (domain.quantities['stage'].boundary_values[5] ==
485                domain.get_conserved_quantities(3, edge=2)[0])
486
487        # Xmomentum
488        # Reflective
489        assert domain.quantities['xmomentum'].boundary_values[0] == 1.0
490        # Dirichlet
491        assert domain.quantities['xmomentum'].boundary_values[1] == 2.
492        # Transmissive
493        assert (domain.quantities['xmomentum'].boundary_values[2] ==
494                domain.get_conserved_quantities(2, edge=0)[1])
495        # Transmissive
496        assert (domain.quantities['xmomentum'].boundary_values[3] ==
497                domain.get_conserved_quantities(2, edge=1)[1])
498        # Transmissive
499        assert (domain.quantities['xmomentum'].boundary_values[4] ==
500                domain.get_conserved_quantities(3, edge=1)[1])
501        # Reflective
502        assert domain.quantities['xmomentum'].boundary_values[5] == -4.0
503
504        # Ymomentum
505        # Reflective
506        assert domain.quantities['ymomentum'].boundary_values[0] == -10.0
507        # Dirichlet
508        assert domain.quantities['ymomentum'].boundary_values[1] == 1.
509        # Transmissive
510        assert domain.quantities['ymomentum'].boundary_values[2] == 30.
511        # Transmissive
512        assert domain.quantities['ymomentum'].boundary_values[3] == 30.
513        # Transmissive
514        assert domain.quantities['ymomentum'].boundary_values[4] == 40.
515        # Reflective
516        assert domain.quantities['ymomentum'].boundary_values[5] == 40.
517
518    def test_boundary_conditionsII(self):
519        a = [0.0, 0.0]
520        b = [0.0, 2.0]
521        c = [2.0, 0.0]
522        d = [0.0, 4.0]
523        e = [2.0, 2.0]
524        f = [4.0, 0.0]
525
526        points = [a, b, c, d, e, f]
527        #             bac,     bce,     ecf,     dbe
528        vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
529        boundary = {(0, 0): 'Third',
530                    (0, 2): 'First',
531                    (2, 0): 'Second',
532                    (2, 1): 'Second',
533                    (3, 1): 'Second',
534                    (3, 2): 'Third',
535                    (0, 1): 'Internal'}
536
537        domain = Domain(points, vertices, boundary)
538        domain.check_integrity()
539
540        domain.set_quantity('stage', [[1,2,3], [5,5,5], [0,0,9], [-6,3,3]])
541
542        domain.set_quantity('xmomentum', [[1,1,1], [2,2,2], [3,3,3], [4,4,4]])
543
544        domain.set_quantity('ymomentum', [[10,10,10], [20,20,20],
545                                          [30,30,30], [40,40,40]])
546
547        D = Dirichlet_boundary([5, 2, 1])
548        T = Transmissive_boundary(domain)
549        R = Reflective_boundary(domain)
550        domain.set_boundary({'First': D, 'Second': T,
551                             'Third': R, 'Internal': None})
552
553        domain.update_boundary()
554        domain.check_integrity()
555
556    def test_boundary_conditionsIII(self):
557        """test_boundary_conditionsIII
558
559        Test Transmissive_stage_zero_momentum_boundary
560        """
561
562        a = [0.0, 0.0]
563        b = [0.0, 2.0]
564        c = [2.0, 0.0]
565        d = [0.0, 4.0]
566        e = [2.0, 2.0]
567        f = [4.0, 0.0]
568
569        points = [a, b, c, d, e, f]
570        #             bac,     bce,     ecf,     dbe
571        vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
572        boundary = {(0, 0): 'Third',
573                    (0, 2): 'First',
574                    (2, 0): 'Second',
575                    (2, 1): 'Second',
576                    (3, 1): 'Second',
577                    (3, 2): 'Third'}
578
579        domain = Domain(points, vertices, boundary)
580        domain.check_integrity()
581
582        domain.set_quantity('stage', [[1,2,3], [5,5,5], [0,0,9], [-6,3,3]])
583
584        domain.set_quantity('xmomentum', [[1,1,1], [2,2,2], [3,3,3], [4,4,4]])
585
586        domain.set_quantity('ymomentum', [[10,10,10], [20,20,20],
587                                          [30,30,30], [40,40,40]])
588
589        D = Dirichlet_boundary([5, 2, 1])
590        T = Transmissive_stage_zero_momentum_boundary(domain)
591        R = Reflective_boundary(domain)
592        domain.set_boundary({'First': D, 'Second': T, 'Third': R})
593
594        domain.update_boundary()
595
596        # Stage
597        # Reflective (2.5)
598        assert domain.quantities['stage'].boundary_values[0] == 2.5
599        assert (domain.quantities['stage'].boundary_values[0] ==
600                domain.get_conserved_quantities(0, edge=0)[0])
601        # Dirichlet
602        assert domain.quantities['stage'].boundary_values[1] == 5.
603        # Transmissive (4.5)
604        assert (domain.quantities['stage'].boundary_values[2] ==
605                domain.get_conserved_quantities(2, edge=0)[0])
606        # Transmissive (4.5)
607        assert (domain.quantities['stage'].boundary_values[3] ==
608                domain.get_conserved_quantities(2, edge=1)[0])
609        # Transmissive (-1.5)
610        assert (domain.quantities['stage'].boundary_values[4] ==
611                domain.get_conserved_quantities(3, edge=1)[0])
612        # Reflective (-1.5)
613        assert (domain.quantities['stage'].boundary_values[5] ==
614                domain.get_conserved_quantities(3, edge=2)[0])
615
616        # Xmomentum
617        # Reflective
618        assert domain.quantities['xmomentum'].boundary_values[0] == 1.0
619        # Dirichlet
620        assert domain.quantities['xmomentum'].boundary_values[1] == 2.
621        assert domain.quantities['xmomentum'].boundary_values[2] == 0.0
622        assert domain.quantities['xmomentum'].boundary_values[3] == 0.0
623        assert domain.quantities['xmomentum'].boundary_values[4] == 0.0
624        # Reflective
625        assert domain.quantities['xmomentum'].boundary_values[5] == -4.0
626
627        # Ymomentum
628        # Reflective
629        assert domain.quantities['ymomentum'].boundary_values[0] == -10.0
630        # Dirichlet
631        assert domain.quantities['ymomentum'].boundary_values[1] == 1.
632        assert domain.quantities['ymomentum'].boundary_values[2] == 0.0
633        assert domain.quantities['ymomentum'].boundary_values[3] == 0.0
634        assert domain.quantities['ymomentum'].boundary_values[4] == 0.0
635        # Reflective
636        assert domain.quantities['ymomentum'].boundary_values[5] == 40.
637
638    def test_boundary_condition_time(self):
639        """test_boundary_condition_time
640
641        This tests that boundary conditions are evaluated
642        using the right time from domain.
643        """
644
645        # Setup computational domain
646        from anuga.abstract_2d_finite_volumes.mesh_factory \
647                import rectangular_cross
648
649        #--------------------------------------------------------------
650        # Setup computational domain
651        #--------------------------------------------------------------
652        N = 5
653        points, vertices, boundary = rectangular_cross(N, N)
654        domain = Domain(points, vertices, boundary)
655
656        #--------------------------------------------------------------
657        # Setup initial conditions
658        #--------------------------------------------------------------
659        domain.set_quantity('elevation', 0.0)
660        domain.set_quantity('friction', 0.0)
661        domain.set_quantity('stage', 0.0)
662
663        #--------------------------------------------------------------
664        # Setup boundary conditions
665        #--------------------------------------------------------------
666        # Time dependent boundary
667        Bt = Time_boundary(domain=domain, f=lambda t: [t, 0.0, 0.0])
668
669        Br = Reflective_boundary(domain)              # Reflective wall
670
671        domain.set_boundary({'left': Bt, 'right': Br, 'top': Br, 'bottom': Br})
672
673        for t in domain.evolve(yieldstep = 10, finaltime = 20.0):
674            q = Bt.evaluate()
675
676            # FIXME (Ole): This test would not have passed in
677            # changeset:5846.
678            msg = 'Time boundary not evaluated correctly'
679            assert num.allclose(t, q[0]), msg
680
681    def test_compute_fluxes0(self):
682        # Do a full triangle and check that fluxes cancel out for
683        # the constant stage case
684
685        a = [0.0, 0.0]
686        b = [0.0, 2.0]
687        c = [2.0, 0.0]
688        d = [0.0, 4.0]
689        e = [2.0, 2.0]
690        f = [4.0, 0.0]
691
692        points = [a, b, c, d, e, f]
693        #             bac,     bce,     ecf,     dbe
694        vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
695
696        domain = Domain(points, vertices)
697        domain.set_quantity('stage', [[2,2,2], [2,2,2], [2,2,2], [2,2,2]])
698        domain.check_integrity()
699
700        assert num.allclose(domain.neighbours,
701                            [[-1,1,-1], [2,3,0], [-1,-1,1],[1,-1,-1]])
702        assert num.allclose(domain.neighbour_edges,
703                            [[-1,2,-1], [2,0,1], [-1,-1,0],[1,-1,-1]])
704
705        zl = zr = 0.     # Assume flat bed
706
707        edgeflux = num.zeros(3, num.float)
708        edgeflux0 = num.zeros(3, num.float)
709        edgeflux1 = num.zeros(3, num.float)
710        edgeflux2 = num.zeros(3, num.float)
711        H0 = 0.0
712
713        # Flux across right edge of volume 1
714        normal = domain.get_normal(1, 0)
715        ql = domain.get_conserved_quantities(vol_id=1, edge=0)
716        qr = domain.get_conserved_quantities(vol_id=2, edge=2)
717        max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux0,
718                                  epsilon, g, H0)
719
720        # Check that flux seen from other triangles is inverse
721        (ql, qr) = (qr, ql)
722        #tmp = qr
723        #qr = ql
724        #ql = tmp
725        normal = domain.get_normal(2, 2)
726        max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux,
727                                  epsilon, g, H0)
728
729        assert num.allclose(edgeflux0 + edgeflux, 0.)
730
731        # Flux across upper edge of volume 1
732        normal = domain.get_normal(1, 1)
733        ql = domain.get_conserved_quantities(vol_id=1, edge=1)
734        qr = domain.get_conserved_quantities(vol_id=3, edge=0)
735        max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux1,
736                                  epsilon, g, H0)
737
738        # Check that flux seen from other triangles is inverse
739        (ql, qr) = (qr, ql)
740        #tmp = qr
741        #qr = ql
742        #ql = tmp
743        normal = domain.get_normal(3, 0)
744        max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux,
745                                  epsilon, g, H0)
746
747        assert num.allclose(edgeflux1 + edgeflux, 0.)
748
749        # Flux across lower left hypotenuse of volume 1
750        normal = domain.get_normal(1, 2)
751        ql = domain.get_conserved_quantities(vol_id=1, edge=2)
752        qr = domain.get_conserved_quantities(vol_id=0, edge=1)
753        max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux2,
754                                  epsilon, g, H0)
755
756        # Check that flux seen from other triangles is inverse
757        (ql, qr) = (qr, ql)
758        #tmp = qr
759        #qr=ql
760        #ql=tmp
761        normal = domain.get_normal(0, 1)
762        max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux,
763                                  epsilon, g, H0)
764        assert num.allclose(edgeflux2 + edgeflux, 0.)
765
766        # Scale by edgelengths, add up anc check that total flux is zero
767        e0 = domain.edgelengths[1, 0]
768        e1 = domain.edgelengths[1, 1]
769        e2 = domain.edgelengths[1, 2]
770
771        assert num.allclose(e0*edgeflux0 + e1*edgeflux1 + e2*edgeflux2, 0.)
772
773        # Now check that compute_flux yields zeros as well
774        domain.compute_fluxes()
775
776        for name in ['stage', 'xmomentum', 'ymomentum']:
777            assert num.allclose(domain.quantities[name].explicit_update[1], 0)
778
779    def test_compute_fluxes1(self):
780        #Use values from previous version
781        a = [0.0, 0.0]
782        b = [0.0, 2.0]
783        c = [2.0, 0.0]
784        d = [0.0, 4.0]
785        e = [2.0, 2.0]
786        f = [4.0, 0.0]
787
788        points = [a, b, c, d, e, f]
789        #             bac,     bce,     ecf,     dbe
790        vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
791
792        domain = Domain(points, vertices)
793        val0 = 2. + 2.0/3
794        val1 = 4. + 4.0/3
795        val2 = 8. + 2.0/3
796        val3 = 2. + 8.0/3
797
798        domain.set_quantity('stage', [[val0, val0, val0], [val1, val1, val1],
799                                      [val2, val2, val2], [val3, val3, val3]])
800        domain.check_integrity()
801
802        zl = zr = 0.    # Assume flat bed
803
804        edgeflux = num.zeros(3, num.float)
805        edgeflux0 = num.zeros(3, num.float)
806        edgeflux1 = num.zeros(3, num.float)
807        edgeflux2 = num.zeros(3, num.float)
808        H0 = 0.0
809
810        # Flux across right edge of volume 1
811        normal = domain.get_normal(1, 0)    # Get normal 0 of triangle 1
812        assert num.allclose(normal, [1, 0])
813
814        ql = domain.get_conserved_quantities(vol_id=1, edge=0)
815        assert num.allclose(ql, [val1, 0, 0])
816
817        qr = domain.get_conserved_quantities(vol_id=2, edge=2)
818        assert num.allclose(qr, [val2, 0, 0])
819
820        max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux0,
821                                  epsilon, g, H0)
822
823        # Flux across edge in the east direction (as per normal vector)
824        assert num.allclose(edgeflux0, [-15.3598804, 253.71111111, 0.])
825        assert num.allclose(max_speed, 9.21592824046)
826
827        #Flux across edge in the west direction (opposite sign for xmomentum)
828        normal_opposite = domain.get_normal(2, 2)   # Get normal 2 of triangle 2
829        assert num.allclose(normal_opposite, [-1, 0])
830
831        max_speed = flux_function(normal_opposite, ql, qr, zl, zr, edgeflux,
832                                  epsilon, g, H0)
833        assert num.allclose(edgeflux, [-15.3598804, -253.71111111, 0.])
834
835        #Flux across upper edge of volume 1
836        normal = domain.get_normal(1, 1)
837        ql = domain.get_conserved_quantities(vol_id=1, edge=1)
838        qr = domain.get_conserved_quantities(vol_id=3, edge=0)
839        max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux1,
840                                  epsilon, g, H0)
841
842        assert num.allclose(edgeflux1, [2.4098563, 0., 123.04444444])
843        assert num.allclose(max_speed, 7.22956891292)
844
845        #Flux across lower left hypotenuse of volume 1
846        normal = domain.get_normal(1, 2)
847        ql = domain.get_conserved_quantities(vol_id=1, edge=2)
848        qr = domain.get_conserved_quantities(vol_id=0, edge=1)
849        max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux2,
850                                  epsilon, g, H0)
851
852        assert num.allclose(edgeflux2, [9.63942522, -61.59685738, -61.59685738])
853        assert num.allclose(max_speed, 7.22956891292)
854
855        #Scale, add up and check that compute_fluxes is correct for vol 1
856        e0 = domain.edgelengths[1, 0]
857        e1 = domain.edgelengths[1, 1]
858        e2 = domain.edgelengths[1, 2]
859
860        total_flux = -(e0*edgeflux0 +
861                       e1*edgeflux1 +
862                       e2*edgeflux2) / domain.areas[1]
863
864        assert num.allclose(total_flux, [-0.68218178, -166.6, -35.93333333])
865
866        domain.compute_fluxes()
867
868        for i, name in enumerate(['stage', 'xmomentum', 'ymomentum']):
869            assert num.allclose(total_flux[i],
870                                domain.quantities[name].explicit_update[1])
871
872        assert num.allclose(domain.quantities['stage'].explicit_update,
873                            [0., -0.68218178, -111.77316251, -35.68522449])
874        assert num.allclose(domain.quantities['xmomentum'].explicit_update,
875                            [-69.68888889, -166.6, 69.68888889, 0])
876        assert num.allclose(domain.quantities['ymomentum'].explicit_update,
877                            [-69.68888889, -35.93333333, 0., 69.68888889])
878
879    def test_compute_fluxes2(self):
880        #Random values, incl momentum
881        a = [0.0, 0.0]
882        b = [0.0, 2.0]
883        c = [2.0, 0.0]
884        d = [0.0, 4.0]
885        e = [2.0, 2.0]
886        f = [4.0, 0.0]
887
888        points = [a, b, c, d, e, f]
889        #             bac,     bce,     ecf,     dbe
890        vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
891
892        domain = Domain(points, vertices)
893        val0 = 2. + 2.0/3
894        val1 = 4. + 4.0/3
895        val2 = 8. + 2.0/3
896        val3 = 2. + 8.0/3
897
898        zl = zr = 0    # Assume flat zero bed
899        edgeflux = num.zeros(3, num.float)
900        edgeflux0 = num.zeros(3, num.float)
901        edgeflux1 = num.zeros(3, num.float)
902        edgeflux2 = num.zeros(3, num.float)
903        H0 = 0.0
904
905        domain.set_quantity('elevation', zl*num.ones((4, 3), num.int)) #array default#
906
907        domain.set_quantity('stage', [[val0, val0-1, val0-2],
908                                      [val1, val1+1, val1],
909                                      [val2, val2-2, val2],
910                                      [val3-0.5, val3, val3]])
911
912        domain.set_quantity('xmomentum',
913                            [[1,2,3], [3,4,5], [1,-1,0], [0,-2,2]])
914
915        domain.set_quantity('ymomentum',
916                            [[1,-1,0], [0,-3,2], [0,1,0], [-1,2,2]])
917
918        domain.check_integrity()
919
920        # Flux across right edge of volume 1
921        normal = domain.get_normal(1, 0)
922        ql = domain.get_conserved_quantities(vol_id=1, edge=0)
923        qr = domain.get_conserved_quantities(vol_id=2, edge=2)
924        max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux0,
925                                  epsilon, g, H0)
926
927        # Flux across upper edge of volume 1
928        normal = domain.get_normal(1, 1)
929        ql = domain.get_conserved_quantities(vol_id=1, edge=1)
930        qr = domain.get_conserved_quantities(vol_id=3, edge=0)
931        max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux1,
932                                  epsilon, g, H0)
933
934        # Flux across lower left hypotenuse of volume 1
935        normal = domain.get_normal(1, 2)
936        ql = domain.get_conserved_quantities(vol_id=1, edge=2)
937        qr = domain.get_conserved_quantities(vol_id=0, edge=1)
938        max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux2,
939                                  epsilon, g, H0)
940
941        # Scale, add up and check that compute_fluxes is correct for vol 1
942        e0 = domain.edgelengths[1, 0]
943        e1 = domain.edgelengths[1, 1]
944        e2 = domain.edgelengths[1, 2]
945
946        total_flux = -(e0*edgeflux0 +
947                       e1*edgeflux1 +
948                       e2*edgeflux2) / domain.areas[1]
949
950        domain.compute_fluxes()
951
952        for i, name in enumerate(['stage', 'xmomentum', 'ymomentum']):
953            assert num.allclose(total_flux[i],
954                                domain.quantities[name].explicit_update[1])
955
956    # FIXME (Ole): Need test like this for fluxes in very shallow water.
957    def test_compute_fluxes3(self):
958        #Random values, incl momentum
959        a = [0.0, 0.0]
960        b = [0.0, 2.0]
961        c = [2.0, 0.0]
962        d = [0.0, 4.0]
963        e = [2.0, 2.0]
964        f = [4.0, 0.0]
965
966        points = [a, b, c, d, e, f]
967        #             bac,     bce,     ecf,     dbe
968        vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
969
970        domain = Domain(points, vertices)
971
972        val0 = 2.+2.0/3
973        val1 = 4.+4.0/3
974        val2 = 8.+2.0/3
975        val3 = 2.+8.0/3
976
977        zl = zr = -3.75    # Assume constant bed (must be less than stage)
978        domain.set_quantity('elevation', zl*num.ones((4, 3), num.int)) #array default#
979
980        edgeflux = num.zeros(3, num.float)
981        edgeflux0 = num.zeros(3, num.float)
982        edgeflux1 = num.zeros(3, num.float)
983        edgeflux2 = num.zeros(3, num.float)
984        H0 = 0.0
985
986        domain.set_quantity('stage', [[val0, val0-1, val0-2],
987                                      [val1, val1+1, val1],
988                                      [val2, val2-2, val2],
989                                      [val3-0.5, val3, val3]])
990
991        domain.set_quantity('xmomentum',
992                            [[1,2,3], [3,4,5], [1,-1,0], [0,-2,2]])
993
994        domain.set_quantity('ymomentum',
995                            [[1,-1,0], [0,-3,2], [0,1,0], [-1,2,2]])
996
997        domain.check_integrity()
998
999        # Flux across right edge of volume 1
1000        normal = domain.get_normal(1, 0)
1001        ql = domain.get_conserved_quantities(vol_id=1, edge=0)
1002        qr = domain.get_conserved_quantities(vol_id=2, edge=2)
1003        max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux0,
1004                                  epsilon, g, H0)
1005
1006        # Flux across upper edge of volume 1
1007        normal = domain.get_normal(1, 1)
1008        ql = domain.get_conserved_quantities(vol_id=1, edge=1)
1009        qr = domain.get_conserved_quantities(vol_id=3, edge=0)
1010        max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux1,
1011                                  epsilon, g, H0)
1012
1013        # Flux across lower left hypotenuse of volume 1
1014        normal = domain.get_normal(1, 2)
1015        ql = domain.get_conserved_quantities(vol_id=1, edge=2)
1016        qr = domain.get_conserved_quantities(vol_id=0, edge=1)
1017        max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux2,
1018                                  epsilon, g, H0)
1019
1020        # Scale, add up and check that compute_fluxes is correct for vol 1
1021        e0 = domain.edgelengths[1, 0]
1022        e1 = domain.edgelengths[1, 1]
1023        e2 = domain.edgelengths[1, 2]
1024
1025        total_flux = -(e0*edgeflux0 +
1026                       e1*edgeflux1 +
1027                       e2*edgeflux2) / domain.areas[1]
1028
1029        domain.compute_fluxes()
1030
1031        for i, name in enumerate(['stage', 'xmomentum', 'ymomentum']):
1032            assert num.allclose(total_flux[i],
1033                                domain.quantities[name].explicit_update[1])
1034
1035    def xtest_catching_negative_heights(self):
1036        #OBSOLETE
1037
1038        a = [0.0, 0.0]
1039        b = [0.0, 2.0]
1040        c = [2.0, 0.0]
1041        d = [0.0, 4.0]
1042        e = [2.0, 2.0]
1043        f = [4.0, 0.0]
1044
1045        points = [a, b, c, d, e, f]
1046        #             bac,     bce,     ecf,     dbe
1047        vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
1048
1049        domain = Domain(points, vertices)
1050        val0 = 2. + 2.0/3
1051        val1 = 4. + 4.0/3
1052        val2 = 8. + 2.0/3
1053        val3 = 2. + 8.0/3
1054
1055        zl = zr = 4    # Too large
1056        domain.set_quantity('elevation', zl*num.ones((4, 3), num.int)) #array default#
1057        domain.set_quantity('stage', [[val0, val0-1, val0-2],
1058                                      [val1, val1+1, val1],
1059                                      [val2, val2-2, val2],
1060                                      [val3-0.5, val3, val3]])
1061
1062        #Should fail
1063        try:
1064            domain.check_integrity()
1065        except:
1066            pass
1067
1068    def test_get_wet_elements(self):
1069        a = [0.0, 0.0]
1070        b = [0.0, 2.0]
1071        c = [2.0, 0.0]
1072        d = [0.0, 4.0]
1073        e = [2.0, 2.0]
1074        f = [4.0, 0.0]
1075
1076        points = [a, b, c, d, e, f]
1077        #             bac,     bce,     ecf,     dbe
1078        vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
1079
1080        domain = Domain(points, vertices)
1081
1082        val0 = 2. + 2.0/3
1083        val1 = 4. + 4.0/3
1084        val2 = 8. + 2.0/3
1085        val3 = 2. + 8.0/3
1086
1087        zl = zr = 5
1088        domain.set_quantity('elevation', zl*num.ones((4, 3), num.int)) #array default#
1089        domain.set_quantity('stage', [[val0, val0-1, val0-2],
1090                                      [val1, val1+1, val1],
1091                                      [val2, val2-2, val2],
1092                                      [val3-0.5, val3, val3]])
1093
1094        domain.check_integrity()
1095
1096        indices = domain.get_wet_elements()
1097        assert num.allclose(indices, [1, 2])
1098
1099        indices = domain.get_wet_elements(indices=[0, 1, 3])
1100        assert num.allclose(indices, [1])
1101
1102    def test_get_maximum_inundation_1(self):
1103        a = [0.0, 0.0]
1104        b = [0.0, 2.0]
1105        c = [2.0, 0.0]
1106        d = [0.0, 4.0]
1107        e = [2.0, 2.0]
1108        f = [4.0, 0.0]
1109
1110        points = [a, b, c, d, e, f]
1111        #             bac,     bce,     ecf,     dbe
1112        vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
1113
1114        domain = Domain(points, vertices)
1115
1116        domain.set_quantity('elevation', lambda x, y: x+2*y)    # 2 4 4 6
1117        domain.set_quantity('stage', 3)
1118
1119        domain.check_integrity()
1120
1121        indices = domain.get_wet_elements()
1122        assert num.allclose(indices, [0])
1123
1124        q = domain.get_maximum_inundation_elevation()
1125        assert num.allclose(q, domain.get_quantity('elevation').\
1126                                   get_values(location='centroids')[0])
1127
1128        x, y = domain.get_maximum_inundation_location()
1129        assert num.allclose([x, y], domain.get_centroid_coordinates()[0])
1130
1131    def test_get_maximum_inundation_2(self):
1132        """test_get_maximum_inundation_2(self)
1133
1134        Test multiple wet cells with same elevation
1135        """
1136
1137        a = [0.0, 0.0]
1138        b = [0.0, 2.0]
1139        c = [2.0, 0.0]
1140        d = [0.0, 4.0]
1141        e = [2.0, 2.0]
1142        f = [4.0, 0.0]
1143
1144        points = [a, b, c, d, e, f]
1145        #             bac,     bce,     ecf,     dbe
1146        vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
1147
1148        domain = Domain(points, vertices)
1149
1150        domain.set_quantity('elevation', lambda x, y: x+2*y)    # 2 4 4 6
1151        domain.set_quantity('stage', 4.1)
1152
1153        domain.check_integrity()
1154
1155        indices = domain.get_wet_elements()
1156        assert num.allclose(indices, [0, 1, 2])
1157
1158        q = domain.get_maximum_inundation_elevation()
1159        assert num.allclose(q, 4)
1160
1161        x, y = domain.get_maximum_inundation_location()
1162        assert num.allclose([x, y], domain.get_centroid_coordinates()[1])
1163
1164    def test_get_maximum_inundation_3(self):
1165        """test_get_maximum_inundation_3(self)
1166
1167        Test of real runup example:
1168        """
1169
1170        from anuga.abstract_2d_finite_volumes.mesh_factory \
1171                import rectangular_cross
1172
1173        initial_runup_height = -0.4
1174        final_runup_height = -0.3
1175
1176        #--------------------------------------------------------------
1177        # Setup computational domain
1178        #--------------------------------------------------------------
1179        N = 5
1180        points, vertices, boundary = rectangular_cross(N, N)
1181        domain = Domain(points, vertices, boundary)
1182        domain.set_maximum_allowed_speed(1.0)
1183
1184        #--------------------------------------------------------------
1185        # Setup initial conditions
1186        #--------------------------------------------------------------
1187        def topography(x, y):
1188            return -x/2                             # linear bed slope
1189
1190        # Use function for elevation
1191        domain.set_quantity('elevation', topography)
1192        domain.set_quantity('friction', 0.)                # Zero friction
1193        # Constant negative initial stage
1194        domain.set_quantity('stage', initial_runup_height)
1195
1196        #--------------------------------------------------------------
1197        # Setup boundary conditions
1198        #--------------------------------------------------------------
1199        Br = Reflective_boundary(domain)                       # Reflective wall
1200        Bd = Dirichlet_boundary([final_runup_height, 0, 0])    # Constant inflow
1201
1202        # All reflective to begin with (still water)
1203        domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
1204
1205        #--------------------------------------------------------------
1206        # Test initial inundation height
1207        #--------------------------------------------------------------
1208
1209        indices = domain.get_wet_elements()
1210        z = domain.get_quantity('elevation').get_values(location='centroids',
1211                                                        indices=indices)
1212        assert num.alltrue(z < initial_runup_height)
1213
1214        q = domain.get_maximum_inundation_elevation()
1215        # First order accuracy
1216        assert num.allclose(q, initial_runup_height, rtol=1.0/N)
1217
1218        x, y = domain.get_maximum_inundation_location()
1219
1220        qref = domain.get_quantity('elevation').\
1221                     get_values(interpolation_points=[[x, y]])
1222        assert num.allclose(q, qref)
1223
1224        wet_elements = domain.get_wet_elements()
1225        wet_elevations = domain.get_quantity('elevation').\
1226                                    get_values(location='centroids',
1227                                               indices=wet_elements)
1228        assert num.alltrue(wet_elevations < initial_runup_height)
1229        assert num.allclose(wet_elevations, z)
1230
1231        #--------------------------------------------------------------
1232        # Let triangles adjust
1233        #--------------------------------------------------------------
1234        for t in domain.evolve(yieldstep = 0.1, finaltime = 1.0):
1235            pass
1236
1237        #--------------------------------------------------------------
1238        # Test inundation height again
1239        #--------------------------------------------------------------
1240        indices = domain.get_wet_elements()
1241        z = domain.get_quantity('elevation').get_values(location='centroids',
1242                                                        indices=indices)
1243
1244        assert num.alltrue(z < initial_runup_height)
1245
1246        q = domain.get_maximum_inundation_elevation()
1247        # First order accuracy
1248        assert num.allclose(q, initial_runup_height, rtol=1.0/N)
1249
1250        x, y = domain.get_maximum_inundation_location()
1251        qref = domain.get_quantity('elevation').\
1252                        get_values(interpolation_points=[[x, y]])
1253        assert num.allclose(q, qref)
1254
1255        #--------------------------------------------------------------
1256        # Update boundary to allow inflow
1257        #--------------------------------------------------------------
1258        domain.set_boundary({'right': Bd})
1259
1260        #--------------------------------------------------------------
1261        # Evolve system through time
1262        #--------------------------------------------------------------
1263        for t in domain.evolve(yieldstep = 0.1, finaltime = 3.0):
1264            pass
1265
1266        #--------------------------------------------------------------
1267        # Test inundation height again
1268        #--------------------------------------------------------------
1269        indices = domain.get_wet_elements()
1270        z = domain.get_quantity('elevation').\
1271                    get_values(location='centroids', indices=indices)
1272
1273        assert num.alltrue(z < final_runup_height)
1274
1275        q = domain.get_maximum_inundation_elevation()
1276        # First order accuracy
1277        assert num.allclose(q, final_runup_height, rtol=1.0/N)
1278
1279        x, y = domain.get_maximum_inundation_location()
1280        qref = domain.get_quantity('elevation').\
1281                        get_values(interpolation_points=[[x, y]])
1282        assert num.allclose(q, qref)
1283
1284        wet_elements = domain.get_wet_elements()
1285        wet_elevations = domain.get_quantity('elevation').\
1286                             get_values(location='centroids',
1287                                        indices=wet_elements)
1288        assert num.alltrue(wet_elevations < final_runup_height)
1289        assert num.allclose(wet_elevations, z)
1290
1291    def test_get_maximum_inundation_from_sww(self):
1292        """test_get_maximum_inundation_from_sww(self)
1293
1294        Test of get_maximum_inundation_elevation()
1295        and get_maximum_inundation_location() from data_manager.py
1296
1297        This is based on test_get_maximum_inundation_3(self) but works with the
1298        stored results instead of with the internal data structure.
1299
1300        This test uses the underlying get_maximum_inundation_data for tests
1301        """
1302
1303        from anuga.abstract_2d_finite_volumes.mesh_factory \
1304                import rectangular_cross
1305        from data_manager import get_maximum_inundation_elevation
1306        from data_manager import get_maximum_inundation_location
1307        from data_manager import get_maximum_inundation_data
1308
1309        initial_runup_height = -0.4
1310        final_runup_height = -0.3
1311
1312        #--------------------------------------------------------------
1313        # Setup computational domain
1314        #--------------------------------------------------------------
1315        N = 10
1316        points, vertices, boundary = rectangular_cross(N, N)
1317        domain = Domain(points, vertices, boundary)
1318        domain.set_name('runup_test')
1319        domain.set_maximum_allowed_speed(1.0)
1320
1321        # FIXME: This works better with old limiters so far
1322        domain.tight_slope_limiters = 0
1323
1324        #--------------------------------------------------------------
1325        # Setup initial conditions
1326        #--------------------------------------------------------------
1327        def topography(x, y):
1328            return -x/2                             # linear bed slope
1329
1330        # Use function for elevation
1331        domain.set_quantity('elevation', topography)
1332        domain.set_quantity('friction', 0.)                # Zero friction
1333        # Constant negative initial stage
1334        domain.set_quantity('stage', initial_runup_height)
1335
1336        #--------------------------------------------------------------
1337        # Setup boundary conditions
1338        #--------------------------------------------------------------
1339        Br = Reflective_boundary(domain)                       # Reflective wall
1340        Bd = Dirichlet_boundary([final_runup_height, 0, 0])    # Constant inflow
1341
1342        # All reflective to begin with (still water)
1343        domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
1344
1345        #--------------------------------------------------------------
1346        # Test initial inundation height
1347        #--------------------------------------------------------------
1348        indices = domain.get_wet_elements()
1349        z = domain.get_quantity('elevation').\
1350                get_values(location='centroids', indices=indices)
1351        assert num.alltrue(z < initial_runup_height)
1352
1353        q_ref = domain.get_maximum_inundation_elevation()
1354        # First order accuracy
1355        assert num.allclose(q_ref, initial_runup_height, rtol=1.0/N)
1356
1357        #--------------------------------------------------------------
1358        # Let triangles adjust
1359        #--------------------------------------------------------------
1360        for t in domain.evolve(yieldstep = 0.1, finaltime = 1.0):
1361            pass
1362
1363        #--------------------------------------------------------------
1364        # Test inundation height again
1365        #--------------------------------------------------------------
1366        q_ref = domain.get_maximum_inundation_elevation()
1367        q = get_maximum_inundation_elevation('runup_test.sww')
1368        msg = 'We got %f, should have been %f' % (q, q_ref)
1369        assert num.allclose(q, q_ref, rtol=1.0/N), msg
1370
1371        q = get_maximum_inundation_elevation('runup_test.sww')
1372        msg = 'We got %f, should have been %f' % (q, initial_runup_height)
1373        assert num.allclose(q, initial_runup_height, rtol = 1.0/N), msg
1374
1375        # Test error condition if time interval is out
1376        try:
1377            q = get_maximum_inundation_elevation('runup_test.sww',
1378                                                 time_interval=[2.0, 3.0])
1379        except ValueError:
1380            pass
1381        else:
1382            msg = 'should have caught wrong time interval'
1383            raise Exception, msg
1384
1385        # Check correct time interval
1386        q, loc = get_maximum_inundation_data('runup_test.sww',
1387                                             time_interval=[0.0, 3.0])
1388        msg = 'We got %f, should have been %f' % (q, initial_runup_height)
1389        assert num.allclose(q, initial_runup_height, rtol = 1.0/N), msg
1390        assert num.allclose(-loc[0]/2, q)    # From topography formula
1391
1392        #--------------------------------------------------------------
1393        # Update boundary to allow inflow
1394        #--------------------------------------------------------------
1395        domain.set_boundary({'right': Bd})
1396
1397        #--------------------------------------------------------------
1398        # Evolve system through time
1399        #--------------------------------------------------------------
1400        q_max = None
1401        for t in domain.evolve(yieldstep = 0.1, finaltime = 3.0,
1402                               skip_initial_step = True):
1403            q = domain.get_maximum_inundation_elevation()
1404            if q > q_max:
1405                q_max = q
1406
1407        #--------------------------------------------------------------
1408        # Test inundation height again
1409        #--------------------------------------------------------------
1410        indices = domain.get_wet_elements()
1411        z = domain.get_quantity('elevation').\
1412                get_values(location='centroids', indices=indices)
1413
1414        assert num.alltrue(z < final_runup_height)
1415
1416        q = domain.get_maximum_inundation_elevation()
1417        # First order accuracy
1418        assert num.allclose(q, final_runup_height, rtol=1.0/N)
1419
1420        q, loc = get_maximum_inundation_data('runup_test.sww',
1421                                             time_interval=[3.0, 3.0])
1422        msg = 'We got %f, should have been %f' % (q, final_runup_height)
1423        assert num.allclose(q, final_runup_height, rtol=1.0/N), msg
1424        assert num.allclose(-loc[0]/2, q)    # From topography formula
1425
1426        q = get_maximum_inundation_elevation('runup_test.sww')
1427        loc = get_maximum_inundation_location('runup_test.sww')
1428        msg = 'We got %f, should have been %f' % (q, q_max)
1429        assert num.allclose(q, q_max, rtol=1.0/N), msg
1430        assert num.allclose(-loc[0]/2, q)    # From topography formula
1431
1432        q = get_maximum_inundation_elevation('runup_test.sww',
1433                                             time_interval=[0, 3])
1434        msg = 'We got %f, should have been %f' % (q, q_max)
1435        assert num.allclose(q, q_max, rtol=1.0/N), msg
1436
1437        # Check polygon mode
1438        # Runup region
1439        polygon = [[0.3, 0.0], [0.9, 0.0], [0.9, 1.0], [0.3, 1.0]]
1440        q = get_maximum_inundation_elevation('runup_test.sww',
1441                                             polygon = polygon,
1442                                             time_interval=[0, 3])
1443        msg = 'We got %f, should have been %f' % (q, q_max)
1444        assert num.allclose(q, q_max, rtol=1.0/N), msg
1445
1446        # Offshore region
1447        polygon = [[0.9, 0.0], [1.0, 0.0], [1.0, 1.0], [0.9, 1.0]]
1448        q, loc = get_maximum_inundation_data('runup_test.sww',
1449                                             polygon = polygon,
1450                                             time_interval=[0, 3])
1451        msg = 'We got %f, should have been %f' % (q, -0.475)
1452        assert num.allclose(q, -0.475, rtol=1.0/N), msg
1453        assert is_inside_polygon(loc, polygon)
1454        assert num.allclose(-loc[0]/2, q)    # From topography formula
1455
1456        # Dry region
1457        polygon = [[0.0, 0.0], [0.4, 0.0], [0.4, 1.0], [0.0, 1.0]]
1458        q, loc = get_maximum_inundation_data('runup_test.sww',
1459                                             polygon = polygon,
1460                                             time_interval=[0, 3])
1461        msg = 'We got %s, should have been None' % (q)
1462        assert q is None, msg
1463        msg = 'We got %s, should have been None' % (loc)
1464        assert loc is None, msg
1465
1466        # Check what happens if no time point is within interval
1467        try:
1468            q = get_maximum_inundation_elevation('runup_test.sww',
1469                                                 time_interval=[2.75, 2.75])
1470        except AssertionError:
1471            pass
1472        else:
1473            msg = 'Time interval should have raised an exception'
1474            raise Exception, msg
1475
1476        # Cleanup
1477        try:
1478            os.remove(domain.get_name() + '.' + domain.format)
1479        except:
1480            pass
1481            #FIXME(Ole): Windows won't allow removal of this
1482
1483    def test_get_flow_through_cross_section_with_geo(self):
1484        """test_get_flow_through_cross_section(self):
1485
1486        Test that the total flow through a cross section can be
1487        correctly obtained at run-time from the ANUGA domain.
1488
1489        This test creates a flat bed with a known flow through it and tests
1490        that the function correctly returns the expected flow.
1491
1492        The specifics are
1493        e = -1 m
1494        u = 2 m/s
1495        h = 2 m
1496        w = 3 m (width of channel)
1497
1498        q = u*h*w = 12 m^3/s
1499
1500        This run tries it with georeferencing and with elevation = -1
1501        """
1502
1503        import time, os
1504        from Scientific.IO.NetCDF import NetCDFFile
1505        from mesh_factory import rectangular
1506
1507        # Create basic mesh (20m x 3m)
1508        width = 3
1509        length = 20
1510        t_end = 1
1511        points, vertices, boundary = rectangular(length, width, length, width)
1512
1513        # Create shallow water domain
1514        domain = Domain(points, vertices, boundary,
1515                        geo_reference=Geo_reference(56, 308500, 6189000))
1516
1517        domain.default_order = 2
1518        domain.set_quantities_to_be_stored(None)
1519
1520        e = -1.0
1521        w = 1.0
1522        h = w-e
1523        u = 2.0
1524        uh = u*h
1525
1526        Br = Reflective_boundary(domain)       # Side walls
1527        Bd = Dirichlet_boundary([w, uh, 0])    # 2 m/s across the 3 m inlet:
1528
1529        # Initial conditions
1530        domain.set_quantity('elevation', e)
1531        domain.set_quantity('stage', w)
1532        domain.set_quantity('xmomentum', uh)
1533        domain.set_boundary({'left': Bd, 'right': Bd, 'top': Br, 'bottom': Br})
1534
1535        # Interpolation points down the middle
1536        I = [[0, width/2.],
1537             [length/2., width/2.],
1538             [length, width/2.]]
1539        interpolation_points = domain.geo_reference.get_absolute(I)
1540
1541        # Shortcuts to quantites
1542        stage = domain.get_quantity('stage')
1543        xmomentum = domain.get_quantity('xmomentum')
1544        ymomentum = domain.get_quantity('ymomentum')
1545
1546        for t in domain.evolve(yieldstep=0.1, finaltime=t_end):
1547            # Check that quantities are they should be in the interior
1548            w_t = stage.get_values(interpolation_points)
1549            uh_t = xmomentum.get_values(interpolation_points)
1550            vh_t = ymomentum.get_values(interpolation_points)
1551
1552            assert num.allclose(w_t, w)
1553            assert num.allclose(uh_t, uh)
1554            assert num.allclose(vh_t, 0.0)
1555
1556            # Check flows through the middle
1557            for i in range(5):
1558                x = length/2. + i*0.23674563    # Arbitrary
1559                cross_section = [[x, 0], [x, width]]
1560
1561                cross_section = domain.geo_reference.get_absolute(cross_section)
1562                Q = domain.get_flow_through_cross_section(cross_section,
1563                                                          verbose=False)
1564
1565                assert num.allclose(Q, uh*width)
1566
1567    def test_get_energy_through_cross_section_with_geo(self):
1568        """test_get_energy_through_cross_section(self):
1569
1570        Test that the total and specific energy through a cross section can be
1571        correctly obtained at run-time from the ANUGA domain.
1572
1573        This test creates a flat bed with a known flow through it and tests
1574        that the function correctly returns the expected energies.
1575
1576        The specifics are
1577        e = -1 m
1578        u = 2 m/s
1579        h = 2 m
1580        w = 3 m (width of channel)
1581
1582        q = u*h*w = 12 m^3/s
1583
1584        This run tries it with georeferencing and with elevation = -1
1585        """
1586
1587        import time, os
1588        from Scientific.IO.NetCDF import NetCDFFile
1589        from mesh_factory import rectangular
1590
1591        # Create basic mesh (20m x 3m)
1592        width = 3
1593        length = 20
1594        t_end = 1
1595        points, vertices, boundary = rectangular(length, width, length, width)
1596
1597        # Create shallow water domain
1598        domain = Domain(points, vertices, boundary,
1599                        geo_reference=Geo_reference(56, 308500, 6189000))
1600
1601        domain.default_order = 2
1602        domain.set_quantities_to_be_stored(None)
1603
1604        e = -1.0
1605        w = 1.0
1606        h = w-e
1607        u = 2.0
1608        uh = u*h
1609
1610        Br = Reflective_boundary(domain)       # Side walls
1611        Bd = Dirichlet_boundary([w, uh, 0])    # 2 m/s across the 3 m inlet:
1612
1613        # Initial conditions
1614        domain.set_quantity('elevation', e)
1615        domain.set_quantity('stage', w)
1616        domain.set_quantity('xmomentum', uh)
1617        domain.set_boundary({'left': Bd, 'right': Bd, 'top': Br, 'bottom': Br})
1618
1619        # Interpolation points down the middle
1620        I = [[0, width/2.],
1621             [length/2., width/2.],
1622             [length, width/2.]]
1623        interpolation_points = domain.geo_reference.get_absolute(I)
1624
1625        # Shortcuts to quantites
1626        stage = domain.get_quantity('stage')
1627        xmomentum = domain.get_quantity('xmomentum')
1628        ymomentum = domain.get_quantity('ymomentum')
1629
1630        for t in domain.evolve(yieldstep=0.1, finaltime=t_end):
1631            # Check that quantities are they should be in the interior
1632            w_t = stage.get_values(interpolation_points)
1633            uh_t = xmomentum.get_values(interpolation_points)
1634            vh_t = ymomentum.get_values(interpolation_points)
1635
1636            assert num.allclose(w_t, w)
1637            assert num.allclose(uh_t, uh)
1638            assert num.allclose(vh_t, 0.0)
1639
1640            # Check energies through the middle
1641            for i in range(5):
1642                x = length/2. + i*0.23674563    # Arbitrary
1643                cross_section = [[x, 0], [x, width]]
1644
1645                cross_section = domain.geo_reference.get_absolute(cross_section)
1646                Es = domain.get_energy_through_cross_section(cross_section,
1647                                                             kind='specific',
1648                                                             verbose=False)
1649
1650                assert num.allclose(Es, h + 0.5*u*u/g)
1651
1652                Et = domain.get_energy_through_cross_section(cross_section,
1653                                                             kind='total',
1654                                                             verbose=False)
1655                assert num.allclose(Et, w + 0.5*u*u/g)
1656
1657    def test_another_runup_example(self):
1658        """test_another_runup_example
1659
1660        Test runup example where actual timeseries at interpolated
1661        points are tested.
1662        """
1663
1664        from anuga.pmesh.mesh_interface import create_mesh_from_regions
1665        from anuga.abstract_2d_finite_volumes.mesh_factory \
1666                import rectangular_cross
1667        from anuga.shallow_water import Domain
1668        from anuga.shallow_water import Reflective_boundary
1669        from anuga.shallow_water import Dirichlet_boundary
1670
1671        #-----------------------------------------------------------------
1672        # Setup computational domain
1673        #-----------------------------------------------------------------
1674        points, vertices, boundary = rectangular_cross(10, 10) # Basic mesh
1675        domain = Domain(points, vertices, boundary) # Create domain
1676        domain.set_default_order(1)       
1677        domain.set_quantities_to_be_stored(None)
1678        domain.set_maximum_allowed_speed(100) #FIXME (Ole): try to remove this
1679
1680        # FIXME (Ole): Need tests where this is commented out
1681        domain.tight_slope_limiters = 0    # Backwards compatibility (14/4/7)
1682        domain.H0 = 0    # Backwards compatibility (6/2/7)
1683        domain.use_centroid_velocities = 0    # Backwards compatibility (7/5/8)
1684
1685        #-----------------------------------------------------------------
1686        # Setup initial conditions
1687        #-----------------------------------------------------------------
1688        def topography(x, y):
1689            return -x/2                              # linear bed slope
1690
1691        domain.set_quantity('elevation', topography)
1692        domain.set_quantity('friction', 0.0)
1693        domain.set_quantity('stage', expression='elevation')
1694
1695        #----------------------------------------------------------------
1696        # Setup boundary conditions
1697        #----------------------------------------------------------------
1698        Br = Reflective_boundary(domain)           # Solid reflective wall
1699        Bd = Dirichlet_boundary([-0.2, 0., 0.])    # Constant boundary values
1700        domain.set_boundary({'left': Br, 'right': Bd, 'top': Br, 'bottom': Br})
1701
1702        #----------------------------------------------------------------
1703        # Evolve system through time
1704        #----------------------------------------------------------------
1705        interpolation_points = [[0.4,0.5], [0.6,0.5], [0.8,0.5], [0.9,0.5]]
1706        gauge_values = []
1707        for _ in interpolation_points:
1708            gauge_values.append([])
1709
1710        time = []
1711        for t in domain.evolve(yieldstep=0.1, finaltime=5.0):
1712            # Record time series at known points
1713            time.append(domain.get_time())
1714
1715            stage = domain.get_quantity('stage')
1716            w = stage.get_values(interpolation_points=interpolation_points)
1717
1718            for i, _ in enumerate(interpolation_points):
1719                gauge_values[i].append(w[i])
1720
1721        #Reference (nautilus 26/6/2008)
1722        G0 = [-0.20000000000000001, -0.20000000000000001, -0.19920600846161715,
1723              -0.19153647344085376, -0.19127622768281194, -0.1770671909675095,
1724              -0.16739412133181927, -0.16196038919122191, -0.15621633053131384,
1725              -0.15130021599977705, -0.13930978857215484, -0.19349274358263582,
1726              -0.19975307598803765, -0.19999897143103357, -0.1999999995532111,
1727              -0.19999999999949952, -0.19999999999949952, -0.19999999999949952,
1728              -0.19997270012494556, -0.19925805948554556, -0.19934513778450533,
1729              -0.19966484196394893, -0.1997352860102084,  -0.19968260481750394,
1730              -0.19980280797303882, -0.19998804881822749, -0.19999999778075916,
1731              -0.19999999999966167, -0.19999999999966167, -0.19999999999966167,
1732              -0.19999999999966167, -0.19999999999966167, -0.19999999999966167,
1733              -0.19999999999966167, -0.19999999999966167, -0.19999999999966167,
1734              -0.19999999999966167, -0.19999999999966167, -0.19999999999966167,
1735              -0.19999999999966167, -0.19999999999966167, -0.19999999999966167,
1736              -0.19999999999966167, -0.19999999999966167, -0.19999999999966167,
1737              -0.19999999999966167, -0.19999999999966167, -0.19999999999966167,
1738              -0.19999999999966167, -0.19999999999966167, -0.19999999999966167]
1739
1740        G1 = [-0.29999999999999993, -0.29999588068034899, -0.29250047332330331,
1741              -0.28335081844518584, -0.26142206997410805, -0.22656028856329835,
1742              -0.21224087216745585, -0.19934324109114465, -0.1889857939783175,
1743              -0.18146311603911383, -0.17401078727434263, -0.15419361061257214,
1744              -0.16225060576782063, -0.19010941396999181, -0.20901161407004412,
1745              -0.21670683975774699, -0.21771386270738891, -0.21481284465869752,
1746              -0.21063120869004387, -0.20669243364582401, -0.20320707386714859,
1747              -0.19984087691926442, -0.19725417448019505, -0.19633783049069981,
1748              -0.19650494599999785, -0.19708316838336942, -0.19779309449413818,
1749              -0.19853070294429562, -0.19917342167307153, -0.19964814677795845,
1750              -0.19991627610824922, -0.20013162970144974, -0.20029864969405509,
1751              -0.20036259676501131, -0.20030682824965193, -0.20016105135750167,
1752              -0.19997664501985918, -0.19980185871568762, -0.19966836175417696,
1753              -0.19958856744312226, -0.19955954696194517, -0.19956950051110917,
1754              -0.19960377086336181, -0.19964885299433241, -0.19969427478531132,
1755              -0.19973301547655564, -0.19976121574277764, -0.19977765285688653,
1756              -0.19978315117522441, -0.19977994634841739, -0.19977101394878494]
1757
1758        G2 = [-0.40000000000000002, -0.39077401254732241, -0.33350466136630474,
1759              -0.29771023004255281, -0.27605439066140897, -0.25986156218997497,
1760              -0.24502185018573647, -0.231792624329521,   -0.21981564668803993,
1761              -0.20870707082936543, -0.19877739883776599, -0.18980922837977957,
1762              -0.17308011674005838, -0.16306400164013773, -0.17798470933304333,
1763              -0.1929554075869116,  -0.20236705191987037, -0.20695767560655007,
1764              -0.20841025876092567, -0.20792102174869989, -0.20655350005579293,
1765              -0.20492002526259828, -0.20310627026780645, -0.20105983335287836,
1766              -0.19937394565794653, -0.19853917506699659, -0.19836389977624452,
1767              -0.19850305023602796, -0.19877764028836831, -0.19910928131034669,
1768              -0.19943705712418805, -0.19970344172958865, -0.19991076989870474,
1769              -0.20010020127747646, -0.20025937787100062, -0.20035087292905965,
1770              -0.20035829921463297, -0.20029606557316171, -0.20019606915365515,
1771              -0.20009096093399206, -0.20000371608204368, -0.19994495432920584,
1772              -0.19991535665176338, -0.19990981826533513, -0.19992106419898723,
1773              -0.19994189853516578, -0.19996624091229293, -0.19998946016985167,
1774              -0.20000842303470234, -0.20002144460718174, -0.20002815561337187]
1775
1776        G3 = [-0.45000000000000001, -0.37631169657400332, -0.33000044342859486,
1777              -0.30586045469008522, -0.28843572253009941, -0.27215308978603808,
1778              -0.25712951540331219, -0.2431608296216613,  -0.23032023651386374,
1779              -0.2184546873456619,  -0.20735123704254332, -0.19740397194806389,
1780              -0.1859829564064375,  -0.16675980728362105, -0.16951575032846536,
1781              -0.1832860872609344,  -0.19485758939241243, -0.20231368291811427,
1782              -0.20625610376074754, -0.20758116241495619, -0.20721445402086161,
1783              -0.20603406830353785, -0.20450262808396991, -0.2026769581185151,
1784              -0.2007401212066364,  -0.19931160535777592, -0.19863606301128725,
1785              -0.19848511940572691, -0.19860091042948352, -0.19885490669377764,
1786              -0.19916542732701112, -0.19946678238611959, -0.19971209594104697,
1787              -0.19991912886512292, -0.2001058430788881,  -0.20024959409472989,
1788              -0.20032160254609382, -0.20031583165752354, -0.20025051539293123,
1789              -0.2001556115816068,  -0.20005952955420872, -0.1999814429561611,
1790              -0.19992977821558131, -0.19990457708664208, -0.19990104785490476,
1791              -0.19991257153954825, -0.19993258231880562, -0.19995548502882532,
1792              -0.19997700760919687, -0.19999429663503748, -0.20000588800248761]
1793
1794        #FIXME (DSG):This is a hack so the anuga install, not precompiled
1795        # works on DSG's win2000, python 2.3
1796        #The problem is the gauge_values[X] are 52 long, not 51.
1797        #
1798        # This was probably fixed by Stephen in changeset:3804
1799        #if len(gauge_values[0]) == 52: gauge_values[0].pop()
1800        #if len(gauge_values[1]) == 52: gauge_values[1].pop()
1801        #if len(gauge_values[2]) == 52: gauge_values[2].pop()
1802        #if len(gauge_values[3]) == 52: gauge_values[3].pop()
1803
1804        assert num.allclose(gauge_values[0], G0)
1805        assert num.allclose(gauge_values[1], G1)
1806        assert num.allclose(gauge_values[2], G2)
1807        assert num.allclose(gauge_values[3], G3)
1808
1809    #####################################################
1810
1811    def test_flux_optimisation(self):
1812        """test_flux_optimisation
1813
1814        Test that fluxes are correctly computed using
1815        dry and still cell exclusions
1816        """
1817
1818        from anuga.config import g
1819        import copy
1820
1821        a = [0.0, 0.0]
1822        b = [0.0, 2.0]
1823        c = [2.0, 0.0]
1824        d = [0.0, 4.0]
1825        e = [2.0, 2.0]
1826        f = [4.0, 0.0]
1827
1828        points = [a, b, c, d, e, f]
1829        #             bac,     bce,     ecf,     dbe
1830        vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
1831
1832        domain = Domain(points, vertices)
1833
1834        #Set up for a gradient of (3,0) at mid triangle (bce)
1835        def slope(x, y):
1836            return 3*x
1837
1838        h = 0.1
1839        def stage(x, y):
1840            return slope(x, y) + h
1841
1842        domain.set_quantity('elevation', slope)
1843        domain.set_quantity('stage', stage)
1844
1845        # Allow slope limiters to work
1846        # (FIXME (Ole): Shouldn't this be automatic in ANUGA?)
1847        domain.distribute_to_vertices_and_edges()
1848
1849        initial_stage = copy.copy(domain.quantities['stage'].vertex_values)
1850
1851        domain.set_boundary({'exterior': Reflective_boundary(domain)})
1852
1853        #  Check that update arrays are initialised to zero
1854        assert num.allclose(domain.get_quantity('stage').explicit_update, 0)
1855        assert num.allclose(domain.get_quantity('xmomentum').explicit_update, 0)
1856        assert num.allclose(domain.get_quantity('ymomentum').explicit_update, 0)
1857
1858        # Get true values
1859        domain.optimise_dry_cells = False
1860        domain.compute_fluxes()
1861        stage_ref = copy.copy(domain.get_quantity('stage').explicit_update)
1862        xmom_ref = copy.copy(domain.get_quantity('xmomentum').explicit_update)
1863        ymom_ref = copy.copy(domain.get_quantity('ymomentum').explicit_update)
1864
1865        # Try with flux optimisation
1866        domain.optimise_dry_cells = True
1867        domain.compute_fluxes()
1868
1869        assert num.allclose(stage_ref,
1870                            domain.get_quantity('stage').explicit_update)
1871        assert num.allclose(xmom_ref,
1872                            domain.get_quantity('xmomentum').explicit_update)
1873        assert num.allclose(ymom_ref,
1874                            domain.get_quantity('ymomentum').explicit_update)
1875
1876    def test_initial_condition(self):
1877        """test_initial_condition
1878
1879        Test that initial condition is output at time == 0 and that
1880        computed values change as system evolves
1881        """
1882
1883        from anuga.config import g
1884        import copy
1885
1886        a = [0.0, 0.0]
1887        b = [0.0, 2.0]
1888        c = [2.0, 0.0]
1889        d = [0.0, 4.0]
1890        e = [2.0, 2.0]
1891        f = [4.0, 0.0]
1892
1893        points = [a, b, c, d, e, f]
1894        #             bac,     bce,     ecf,     dbe
1895        vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
1896
1897        domain = Domain(points, vertices)
1898
1899        #Set up for a gradient of (3,0) at mid triangle (bce)
1900        def slope(x, y):
1901            return 3*x
1902
1903        h = 0.1
1904        def stage(x, y):
1905            return slope(x, y) + h
1906
1907        domain.set_quantity('elevation', slope)
1908        domain.set_quantity('stage', stage)
1909
1910        # Allow slope limiters to work
1911        # (FIXME (Ole): Shouldn't this be automatic in ANUGA?)
1912        domain.distribute_to_vertices_and_edges()
1913
1914        initial_stage = copy.copy(domain.quantities['stage'].vertex_values)
1915
1916        domain.set_boundary({'exterior': Reflective_boundary(domain)})
1917
1918        domain.optimise_dry_cells = True
1919
1920        #Evolution
1921        for t in domain.evolve(yieldstep=0.5, finaltime=2.0):
1922            stage = domain.quantities['stage'].vertex_values
1923
1924            if t == 0.0:
1925                assert num.allclose(stage, initial_stage)
1926            else:
1927                assert not num.allclose(stage, initial_stage)
1928
1929        os.remove(domain.get_name() + '.sww')
1930
1931    #####################################################
1932
1933    def test_gravity(self):
1934        #Assuming no friction
1935
1936        from anuga.config import g
1937
1938        a = [0.0, 0.0]
1939        b = [0.0, 2.0]
1940        c = [2.0, 0.0]
1941        d = [0.0, 4.0]
1942        e = [2.0, 2.0]
1943        f = [4.0, 0.0]
1944
1945        points = [a, b, c, d, e, f]
1946        #             bac,     bce,     ecf,     dbe
1947        vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
1948
1949        domain = Domain(points, vertices)
1950
1951        #Set up for a gradient of (3,0) at mid triangle (bce)
1952        def slope(x, y):
1953            return 3*x
1954
1955        h = 0.1
1956        def stage(x, y):
1957            return slope(x, y) + h
1958
1959        domain.set_quantity('elevation', slope)
1960        domain.set_quantity('stage', stage)
1961
1962        for name in domain.conserved_quantities:
1963            assert num.allclose(domain.quantities[name].explicit_update, 0)
1964            assert num.allclose(domain.quantities[name].semi_implicit_update, 0)
1965
1966        domain.compute_forcing_terms()
1967
1968        assert num.allclose(domain.quantities['stage'].explicit_update, 0)
1969        assert num.allclose(domain.quantities['xmomentum'].explicit_update,
1970                            -g*h*3)
1971        assert num.allclose(domain.quantities['ymomentum'].explicit_update, 0)
1972
1973    def test_manning_friction(self):
1974        from anuga.config import g
1975
1976        a = [0.0, 0.0]
1977        b = [0.0, 2.0]
1978        c = [2.0, 0.0]
1979        d = [0.0, 4.0]
1980        e = [2.0, 2.0]
1981        f = [4.0, 0.0]
1982
1983        points = [a, b, c, d, e, f]
1984        #             bac,     bce,     ecf,     dbe
1985        vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
1986
1987        domain = Domain(points, vertices)
1988
1989        #Set up for a gradient of (3,0) at mid triangle (bce)
1990        def slope(x, y):
1991            return 3*x
1992
1993        h = 0.1
1994        def stage(x, y):
1995            return slope(x, y) + h
1996
1997        eta = 0.07
1998        domain.set_quantity('elevation', slope)
1999        domain.set_quantity('stage', stage)
2000        domain.set_quantity('friction', eta)
2001
2002        for name in domain.conserved_quantities:
2003            assert num.allclose(domain.quantities[name].explicit_update, 0)
2004            assert num.allclose(domain.quantities[name].semi_implicit_update, 0)
2005
2006        domain.compute_forcing_terms()
2007
2008        assert num.allclose(domain.quantities['stage'].explicit_update, 0)
2009        assert num.allclose(domain.quantities['xmomentum'].explicit_update,
2010                            -g*h*3)
2011        assert num.allclose(domain.quantities['ymomentum'].explicit_update, 0)
2012
2013        assert num.allclose(domain.quantities['stage'].semi_implicit_update, 0)
2014        assert num.allclose(domain.quantities['xmomentum'].semi_implicit_update,
2015                            0)
2016        assert num.allclose(domain.quantities['ymomentum'].semi_implicit_update,
2017                            0)
2018
2019        #Create some momentum for friction to work with
2020        domain.set_quantity('xmomentum', 1)
2021        S = -g*eta**2 / h**(7.0/3)
2022
2023        domain.compute_forcing_terms()
2024        assert num.allclose(domain.quantities['stage'].semi_implicit_update, 0)
2025        assert num.allclose(domain.quantities['xmomentum'].semi_implicit_update,
2026                            S)
2027        assert num.allclose(domain.quantities['ymomentum'].semi_implicit_update,
2028                            0)
2029
2030        #A more complex example
2031        domain.quantities['stage'].semi_implicit_update[:] = 0.0
2032        domain.quantities['xmomentum'].semi_implicit_update[:] = 0.0
2033        domain.quantities['ymomentum'].semi_implicit_update[:] = 0.0
2034
2035        domain.set_quantity('xmomentum', 3)
2036        domain.set_quantity('ymomentum', 4)
2037
2038        S = -g*eta**2*5 / h**(7.0/3)
2039
2040        domain.compute_forcing_terms()
2041
2042        assert num.allclose(domain.quantities['stage'].semi_implicit_update, 0)
2043        assert num.allclose(domain.quantities['xmomentum'].semi_implicit_update,
2044                            3*S)
2045        assert num.allclose(domain.quantities['ymomentum'].semi_implicit_update,
2046                            4*S)
2047
2048    def test_constant_wind_stress(self):
2049        from anuga.config import rho_a, rho_w, eta_w
2050        from math import pi, cos, sin
2051
2052        a = [0.0, 0.0]
2053        b = [0.0, 2.0]
2054        c = [2.0, 0.0]
2055        d = [0.0, 4.0]
2056        e = [2.0, 2.0]
2057        f = [4.0, 0.0]
2058
2059        points = [a, b, c, d, e, f]
2060        #             bac,     bce,     ecf,     dbe
2061        vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
2062
2063        domain = Domain(points, vertices)
2064
2065        #Flat surface with 1m of water
2066        domain.set_quantity('elevation', 0)
2067        domain.set_quantity('stage', 1.0)
2068        domain.set_quantity('friction', 0)
2069
2070        Br = Reflective_boundary(domain)
2071        domain.set_boundary({'exterior': Br})
2072
2073        #Setup only one forcing term, constant wind stress
2074        s = 100
2075        phi = 135
2076        domain.forcing_terms = []
2077        domain.forcing_terms.append(Wind_stress(s, phi))
2078
2079        domain.compute_forcing_terms()
2080
2081        const = eta_w*rho_a / rho_w
2082
2083        #Convert to radians
2084        phi = phi*pi / 180
2085
2086        #Compute velocity vector (u, v)
2087        u = s*cos(phi)
2088        v = s*sin(phi)
2089
2090        #Compute wind stress
2091        S = const * num.sqrt(u**2 + v**2)
2092
2093        assert num.allclose(domain.quantities['stage'].explicit_update, 0)
2094        assert num.allclose(domain.quantities['xmomentum'].explicit_update, S*u)
2095        assert num.allclose(domain.quantities['ymomentum'].explicit_update, S*v)
2096
2097    def test_variable_wind_stress(self):
2098        from anuga.config import rho_a, rho_w, eta_w
2099        from math import pi, cos, sin
2100
2101        a = [0.0, 0.0]
2102        b = [0.0, 2.0]
2103        c = [2.0, 0.0]
2104        d = [0.0, 4.0]
2105        e = [2.0, 2.0]
2106        f = [4.0, 0.0]
2107
2108        points = [a, b, c, d, e, f]
2109        #             bac,     bce,     ecf,     dbe
2110        vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
2111
2112        domain = Domain(points, vertices)
2113
2114        #Flat surface with 1m of water
2115        domain.set_quantity('elevation', 0)
2116        domain.set_quantity('stage', 1.0)
2117        domain.set_quantity('friction', 0)
2118
2119        Br = Reflective_boundary(domain)
2120        domain.set_boundary({'exterior': Br})
2121
2122        domain.time = 5.54    # Take a random time (not zero)
2123
2124        #Setup only one forcing term, constant wind stress
2125        s = 100
2126        phi = 135
2127        domain.forcing_terms = []
2128        domain.forcing_terms.append(Wind_stress(s=speed, phi=angle))
2129
2130        domain.compute_forcing_terms()
2131
2132        #Compute reference solution
2133        const = eta_w*rho_a / rho_w
2134
2135        N = len(domain)    # number_of_triangles
2136
2137        xc = domain.get_centroid_coordinates()
2138        t = domain.time
2139
2140        x = xc[:,0]
2141        y = xc[:,1]
2142        s_vec = speed(t,x,y)
2143        phi_vec = angle(t,x,y)
2144
2145        for k in range(N):
2146            # Convert to radians
2147            phi = phi_vec[k]*pi / 180
2148            s = s_vec[k]
2149
2150            # Compute velocity vector (u, v)
2151            u = s*cos(phi)
2152            v = s*sin(phi)
2153
2154            # Compute wind stress
2155            S = const * num.sqrt(u**2 + v**2)
2156
2157            assert num.allclose(domain.quantities['stage'].explicit_update[k],
2158                                0)
2159            assert num.allclose(domain.quantities['xmomentum'].\
2160                                     explicit_update[k],
2161                                S*u)
2162            assert num.allclose(domain.quantities['ymomentum'].\
2163                                     explicit_update[k],
2164                                S*v)
2165
2166    def test_windfield_from_file(self):
2167        import time
2168        from anuga.config import rho_a, rho_w, eta_w
2169        from math import pi, cos, sin
2170        from anuga.config import time_format
2171        from anuga.abstract_2d_finite_volumes.util import file_function
2172
2173        a = [0.0, 0.0]
2174        b = [0.0, 2.0]
2175        c = [2.0, 0.0]
2176        d = [0.0, 4.0]
2177        e = [2.0, 2.0]
2178        f = [4.0, 0.0]
2179
2180        points = [a, b, c, d, e, f]
2181        #             bac,     bce,     ecf,     dbe
2182        vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
2183
2184        domain = Domain(points, vertices)
2185
2186        # Flat surface with 1m of water
2187        domain.set_quantity('elevation', 0)
2188        domain.set_quantity('stage', 1.0)
2189        domain.set_quantity('friction', 0)
2190
2191        Br = Reflective_boundary(domain)
2192        domain.set_boundary({'exterior': Br})
2193
2194        domain.time = 7    # Take a time that is represented in file (not zero)
2195
2196        # Write wind stress file (ensure that domain.time is covered)
2197        # Take x=1 and y=0
2198        filename = 'test_windstress_from_file'
2199        start = time.mktime(time.strptime('2000', '%Y'))
2200        fid = open(filename + '.txt', 'w')
2201        dt = 1    # One second interval
2202        t = 0.0
2203        while t <= 10.0:
2204            t_string = time.strftime(time_format, time.gmtime(t+start))
2205
2206            fid.write('%s, %f %f\n' %
2207                      (t_string, speed(t,[1],[0])[0], angle(t,[1],[0])[0]))
2208            t += dt
2209
2210        fid.close()
2211
2212        # Convert ASCII file to NetCDF (Which is what we really like!)
2213        from data_manager import timefile2netcdf
2214
2215        timefile2netcdf(filename)
2216        os.remove(filename + '.txt')
2217
2218        # Setup wind stress
2219        F = file_function(filename + '.tms',
2220                          quantities=['Attribute0', 'Attribute1'])
2221        os.remove(filename + '.tms')
2222
2223        W = Wind_stress(F)
2224
2225        domain.forcing_terms = []
2226        domain.forcing_terms.append(W)
2227
2228        domain.compute_forcing_terms()
2229
2230        # Compute reference solution
2231        const = eta_w*rho_a / rho_w
2232
2233        N = len(domain)    # number_of_triangles
2234
2235        t = domain.time
2236
2237        s = speed(t, [1], [0])[0]
2238        phi = angle(t, [1], [0])[0]
2239
2240        # Convert to radians
2241        phi = phi*pi / 180
2242
2243        # Compute velocity vector (u, v)
2244        u = s*cos(phi)
2245        v = s*sin(phi)
2246
2247        # Compute wind stress
2248        S = const * num.sqrt(u**2 + v**2)
2249
2250        for k in range(N):
2251            assert num.allclose(domain.quantities['stage'].explicit_update[k],
2252                                0)
2253            assert num.allclose(domain.quantities['xmomentum'].\
2254                                    explicit_update[k],
2255                                S*u)
2256            assert num.allclose(domain.quantities['ymomentum'].\
2257                                    explicit_update[k],
2258                                S*v)
2259
2260    def test_windfield_from_file_seconds(self):
2261        import time
2262        from anuga.config import rho_a, rho_w, eta_w
2263        from math import pi, cos, sin
2264        from anuga.config import time_format
2265        from anuga.abstract_2d_finite_volumes.util import file_function
2266
2267        a = [0.0, 0.0]
2268        b = [0.0, 2.0]
2269        c = [2.0, 0.0]
2270        d = [0.0, 4.0]
2271        e = [2.0, 2.0]
2272        f = [4.0, 0.0]
2273
2274        points = [a, b, c, d, e, f]
2275        #             bac,     bce,     ecf,     dbe
2276        vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
2277
2278        domain = Domain(points, vertices)
2279
2280        # Flat surface with 1m of water
2281        domain.set_quantity('elevation', 0)
2282        domain.set_quantity('stage', 1.0)
2283        domain.set_quantity('friction', 0)
2284
2285        Br = Reflective_boundary(domain)
2286        domain.set_boundary({'exterior': Br})
2287
2288        domain.time = 7    # Take a time that is represented in file (not zero)
2289
2290        # Write wind stress file (ensure that domain.time is covered)
2291        # Take x=1 and y=0
2292        filename = 'test_windstress_from_file'
2293        start = time.mktime(time.strptime('2000', '%Y'))
2294        fid = open(filename + '.txt', 'w')
2295        dt = 0.5    # Half second interval
2296        t = 0.0
2297        while t <= 10.0:
2298            fid.write('%s, %f %f\n'
2299                      % (str(t), speed(t, [1], [0])[0], angle(t, [1], [0])[0]))
2300            t += dt
2301
2302        fid.close()
2303
2304        # Convert ASCII file to NetCDF (Which is what we really like!)
2305        from data_manager import timefile2netcdf
2306
2307        timefile2netcdf(filename, time_as_seconds=True)
2308        os.remove(filename + '.txt')
2309
2310        # Setup wind stress
2311        F = file_function(filename + '.tms',
2312                          quantities=['Attribute0', 'Attribute1'])
2313        os.remove(filename + '.tms')
2314
2315        W = Wind_stress(F)
2316
2317        domain.forcing_terms = []
2318        domain.forcing_terms.append(W)
2319
2320        domain.compute_forcing_terms()
2321
2322        # Compute reference solution
2323        const = eta_w*rho_a / rho_w
2324
2325        N = len(domain)    # number_of_triangles
2326
2327        t = domain.time
2328
2329        s = speed(t, [1], [0])[0]
2330        phi = angle(t, [1], [0])[0]
2331
2332        # Convert to radians
2333        phi = phi*pi / 180
2334
2335        # Compute velocity vector (u, v)
2336        u = s*cos(phi)
2337        v = s*sin(phi)
2338
2339        # Compute wind stress
2340        S = const * num.sqrt(u**2 + v**2)
2341
2342        for k in range(N):
2343            assert num.allclose(domain.quantities['stage'].explicit_update[k],
2344                                0)
2345            assert num.allclose(domain.quantities['xmomentum'].\
2346                                    explicit_update[k],
2347                                S*u)
2348            assert num.allclose(domain.quantities['ymomentum'].\
2349                                    explicit_update[k],
2350                                S*v)
2351
2352    def test_wind_stress_error_condition(self):
2353        """Test that windstress reacts properly when forcing functions
2354        are wrong - e.g. returns a scalar
2355        """
2356
2357        from math import pi, cos, sin
2358        from anuga.config import rho_a, rho_w, eta_w
2359
2360        a = [0.0, 0.0]
2361        b = [0.0, 2.0]
2362        c = [2.0, 0.0]
2363        d = [0.0, 4.0]
2364        e = [2.0, 2.0]
2365        f = [4.0, 0.0]
2366
2367        points = [a, b, c, d, e, f]
2368        #             bac,     bce,     ecf,     dbe
2369        vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
2370
2371        domain = Domain(points, vertices)
2372
2373        # Flat surface with 1m of water
2374        domain.set_quantity('elevation', 0)
2375        domain.set_quantity('stage', 1.0)
2376        domain.set_quantity('friction', 0)
2377
2378        Br = Reflective_boundary(domain)
2379        domain.set_boundary({'exterior': Br})
2380
2381        domain.time = 5.54    # Take a random time (not zero)
2382
2383        # Setup only one forcing term, bad func
2384        domain.forcing_terms = []
2385
2386        try:
2387            domain.forcing_terms.append(Wind_stress(s=scalar_func_list,
2388                                                    phi=angle))
2389        except AssertionError:
2390            pass
2391        else:
2392            msg = 'Should have raised exception'
2393            raise Exception, msg
2394
2395        try:
2396            domain.forcing_terms.append(Wind_stress(s=speed, phi=scalar_func))
2397        except Exception:
2398            pass
2399        else:
2400            msg = 'Should have raised exception'
2401            raise Exception, msg
2402
2403        try:
2404            domain.forcing_terms.append(Wind_stress(s=speed, phi='xx'))
2405        except:
2406            pass
2407        else:
2408            msg = 'Should have raised exception'
2409            raise Exception, msg
2410
2411    def test_rainfall(self):
2412        from math import pi, cos, sin
2413
2414        a = [0.0, 0.0]
2415        b = [0.0, 2.0]
2416        c = [2.0, 0.0]
2417        d = [0.0, 4.0]
2418        e = [2.0, 2.0]
2419        f = [4.0, 0.0]
2420
2421        points = [a, b, c, d, e, f]
2422        #             bac,     bce,     ecf,     dbe
2423        vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
2424
2425        domain = Domain(points, vertices)
2426
2427        # Flat surface with 1m of water
2428        domain.set_quantity('elevation', 0)
2429        domain.set_quantity('stage', 1.0)
2430        domain.set_quantity('friction', 0)
2431
2432        Br = Reflective_boundary(domain)
2433        domain.set_boundary({'exterior': Br})
2434
2435        # Setup only one forcing term, constant rainfall
2436        domain.forcing_terms = []
2437        domain.forcing_terms.append(Rainfall(domain, rate=2.0))
2438
2439        domain.compute_forcing_terms()
2440        assert num.allclose(domain.quantities['stage'].explicit_update,
2441                            2.0/1000)
2442
2443    def test_rainfall_restricted_by_polygon(self):
2444        from math import pi, cos, sin
2445
2446        a = [0.0, 0.0]
2447        b = [0.0, 2.0]
2448        c = [2.0, 0.0]
2449        d = [0.0, 4.0]
2450        e = [2.0, 2.0]
2451        f = [4.0, 0.0]
2452
2453        points = [a, b, c, d, e, f]
2454        #             bac,     bce,     ecf,     dbe
2455        vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
2456
2457        domain = Domain(points, vertices)
2458
2459        # Flat surface with 1m of water
2460        domain.set_quantity('elevation', 0)
2461        domain.set_quantity('stage', 1.0)
2462        domain.set_quantity('friction', 0)
2463
2464        Br = Reflective_boundary(domain)
2465        domain.set_boundary({'exterior': Br})
2466
2467        # Setup only one forcing term, constant rainfall
2468        # restricted to a polygon enclosing triangle #1 (bce)
2469        domain.forcing_terms = []
2470        R = Rainfall(domain, rate=2.0, polygon=[[1,1], [2,1], [2,2], [1,2]])
2471
2472        assert num.allclose(R.exchange_area, 1)
2473
2474        domain.forcing_terms.append(R)
2475
2476        domain.compute_forcing_terms()
2477
2478        assert num.allclose(domain.quantities['stage'].explicit_update[1],
2479                            2.0/1000)
2480        assert num.allclose(domain.quantities['stage'].explicit_update[0], 0)
2481        assert num.allclose(domain.quantities['stage'].explicit_update[2:], 0)
2482
2483    def test_time_dependent_rainfall_restricted_by_polygon(self):
2484        a = [0.0, 0.0]
2485        b = [0.0, 2.0]
2486        c = [2.0, 0.0]
2487        d = [0.0, 4.0]
2488        e = [2.0, 2.0]
2489        f = [4.0, 0.0]
2490
2491        points = [a, b, c, d, e, f]
2492        #             bac,     bce,     ecf,     dbe
2493        vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
2494
2495        domain = Domain(points, vertices)
2496
2497        # Flat surface with 1m of water
2498        domain.set_quantity('elevation', 0)
2499        domain.set_quantity('stage', 1.0)
2500        domain.set_quantity('friction', 0)
2501
2502        Br = Reflective_boundary(domain)
2503        domain.set_boundary({'exterior': Br})
2504
2505        # Setup only one forcing term, time dependent rainfall
2506        # restricted to a polygon enclosing triangle #1 (bce)
2507        domain.forcing_terms = []
2508        R = Rainfall(domain, rate=lambda t: 3*t + 7,
2509                     polygon=[[1,1], [2,1], [2,2], [1,2]])
2510
2511        assert num.allclose(R.exchange_area, 1)
2512
2513        domain.forcing_terms.append(R)
2514
2515        domain.time = 10.
2516
2517        domain.compute_forcing_terms()
2518
2519        assert num.allclose(domain.quantities['stage'].explicit_update[1],
2520                            (3*domain.time + 7)/1000)
2521        assert num.allclose(domain.quantities['stage'].explicit_update[0], 0)
2522        assert num.allclose(domain.quantities['stage'].explicit_update[2:], 0)
2523
2524    def test_time_dependent_rainfall_using_starttime(self):
2525        rainfall_poly = ensure_numeric([[1,1], [2,1], [2,2], [1,2]], num.float)
2526
2527        a = [0.0, 0.0]
2528        b = [0.0, 2.0]
2529        c = [2.0, 0.0]
2530        d = [0.0, 4.0]
2531        e = [2.0, 2.0]
2532        f = [4.0, 0.0]
2533
2534        points = [a, b, c, d, e, f]
2535        #             bac,     bce,     ecf,     dbe
2536        vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
2537
2538        domain = Domain(points, vertices)
2539
2540        # Flat surface with 1m of water
2541        domain.set_quantity('elevation', 0)
2542        domain.set_quantity('stage', 1.0)
2543        domain.set_quantity('friction', 0)
2544
2545        Br = Reflective_boundary(domain)
2546        domain.set_boundary({'exterior': Br})
2547
2548        # Setup only one forcing term, time dependent rainfall
2549        # restricted to a polygon enclosing triangle #1 (bce)
2550        domain.forcing_terms = []
2551        R = Rainfall(domain, rate=lambda t: 3*t + 7,
2552                     polygon=rainfall_poly)
2553
2554        assert num.allclose(R.exchange_area, 1)
2555
2556        domain.forcing_terms.append(R)
2557
2558        # This will test that time used in the forcing function takes
2559        # startime into account.
2560        domain.starttime = 5.0
2561
2562        domain.time = 7.
2563
2564        domain.compute_forcing_terms()
2565
2566        assert num.allclose(domain.quantities['stage'].explicit_update[1],
2567                            (3*domain.get_time() + 7)/1000)
2568        assert num.allclose(domain.quantities['stage'].explicit_update[1],
2569                            (3*(domain.time + domain.starttime) + 7)/1000)
2570
2571        # Using internal time her should fail
2572        assert not num.allclose(domain.quantities['stage'].explicit_update[1],
2573                                (3*domain.time + 7)/1000)
2574
2575        assert num.allclose(domain.quantities['stage'].explicit_update[0], 0)
2576        assert num.allclose(domain.quantities['stage'].explicit_update[2:], 0)
2577
2578    def test_time_dependent_rainfall_using_georef(self):
2579        """test_time_dependent_rainfall_using_georef
2580
2581        This will also test the General forcing term using georef
2582        """
2583
2584        # Mesh in zone 56 (absolute coords)
2585        x0 = 314036.58727982
2586        y0 = 6224951.2960092
2587
2588        rainfall_poly = ensure_numeric([[1,1], [2,1], [2,2], [1,2]], num.float)
2589        rainfall_poly += [x0, y0]
2590
2591        a = [0.0, 0.0]
2592        b = [0.0, 2.0]
2593        c = [2.0, 0.0]
2594        d = [0.0, 4.0]
2595        e = [2.0, 2.0]
2596        f = [4.0, 0.0]
2597
2598        points = [a, b, c, d, e, f]
2599        #             bac,     bce,     ecf,     dbe
2600        vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
2601
2602        domain = Domain(points, vertices,
2603                        geo_reference=Geo_reference(56, x0, y0))
2604
2605        # Flat surface with 1m of water
2606        domain.set_quantity('elevation', 0)
2607        domain.set_quantity('stage', 1.0)
2608        domain.set_quantity('friction', 0)
2609
2610        Br = Reflective_boundary(domain)
2611        domain.set_boundary({'exterior': Br})
2612
2613        # Setup only one forcing term, time dependent rainfall
2614        # restricted to a polygon enclosing triangle #1 (bce)
2615        domain.forcing_terms = []
2616        R = Rainfall(domain, rate=lambda t: 3*t + 7, polygon=rainfall_poly)
2617
2618        assert num.allclose(R.exchange_area, 1)
2619
2620        domain.forcing_terms.append(R)
2621
2622        # This will test that time used in the forcing function takes
2623        # startime into account.
2624        domain.starttime = 5.0
2625
2626        domain.time = 7.
2627
2628        domain.compute_forcing_terms()
2629
2630        assert num.allclose(domain.quantities['stage'].explicit_update[1],
2631                            (3*domain.get_time() + 7)/1000)
2632        assert num.allclose(domain.quantities['stage'].explicit_update[1],
2633                            (3*(domain.time + domain.starttime) + 7)/1000)
2634
2635        # Using internal time her should fail
2636        assert not num.allclose(domain.quantities['stage'].explicit_update[1],
2637                                (3*domain.time + 7)/1000)
2638
2639        assert num.allclose(domain.quantities['stage'].explicit_update[0], 0)
2640        assert num.allclose(domain.quantities['stage'].explicit_update[2:], 0)
2641
2642    def test_time_dependent_rainfall_restricted_by_polygon_with_default(self):
2643        """
2644        Test that default rainfall can be used when given rate runs out of data.
2645        """
2646
2647        a = [0.0, 0.0]
2648        b = [0.0, 2.0]
2649        c = [2.0, 0.0]
2650        d = [0.0, 4.0]
2651        e = [2.0, 2.0]
2652        f = [4.0, 0.0]
2653
2654        points = [a, b, c, d, e, f]
2655        #             bac,     bce,     ecf,     dbe
2656        vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
2657
2658        domain = Domain(points, vertices)
2659
2660        # Flat surface with 1m of water
2661        domain.set_quantity('elevation', 0)
2662        domain.set_quantity('stage', 1.0)
2663        domain.set_quantity('friction', 0)
2664
2665        Br = Reflective_boundary(domain)
2666        domain.set_boundary({'exterior': Br})
2667
2668        # Setup only one forcing term, time dependent rainfall
2669        # that expires at t==20
2670        from anuga.fit_interpolate.interpolate import Modeltime_too_late
2671
2672        def main_rate(t):
2673            if t > 20:
2674                msg = 'Model time exceeded.'
2675                raise Modeltime_too_late, msg
2676            else:
2677                return 3*t + 7
2678
2679        domain.forcing_terms = []
2680        R = Rainfall(domain, rate=main_rate,
2681                     polygon=[[1,1], [2,1], [2,2], [1,2]], default_rate=5.0)
2682
2683        assert num.allclose(R.exchange_area, 1)
2684
2685        domain.forcing_terms.append(R)
2686
2687        domain.time = 10.
2688
2689        domain.compute_forcing_terms()
2690
2691        assert num.allclose(domain.quantities['stage'].explicit_update[1],
2692                            (3*domain.time+7)/1000)
2693        assert num.allclose(domain.quantities['stage'].explicit_update[0], 0)
2694        assert num.allclose(domain.quantities['stage'].explicit_update[2:], 0)
2695
2696        domain.time = 100.
2697        domain.quantities['stage'].explicit_update[:] = 0.0     # Reset
2698        domain.compute_forcing_terms()
2699
2700        assert num.allclose(domain.quantities['stage'].explicit_update[1],
2701                            5.0/1000) # Default value
2702        assert num.allclose(domain.quantities['stage'].explicit_update[0], 0)
2703        assert num.allclose(domain.quantities['stage'].explicit_update[2:], 0)
2704
2705    def test_rainfall_forcing_with_evolve(self):
2706        """test_rainfall_forcing_with_evolve
2707
2708        Test how forcing terms are called within evolve
2709        """
2710
2711        # FIXME(Ole): This test is just to experiment
2712
2713        a = [0.0, 0.0]
2714        b = [0.0, 2.0]
2715        c = [2.0, 0.0]
2716        d = [0.0, 4.0]
2717        e = [2.0, 2.0]
2718        f = [4.0, 0.0]
2719
2720        points = [a, b, c, d, e, f]
2721        #             bac,     bce,     ecf,     dbe
2722        vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
2723
2724        domain = Domain(points, vertices)
2725
2726        # Flat surface with 1m of water
2727        domain.set_quantity('elevation', 0)
2728        domain.set_quantity('stage', 1.0)
2729        domain.set_quantity('friction', 0)
2730
2731        Br = Reflective_boundary(domain)
2732        domain.set_boundary({'exterior': Br})
2733
2734        # Setup only one forcing term, time dependent rainfall
2735        # that expires at t==20
2736        from anuga.fit_interpolate.interpolate import Modeltime_too_late
2737
2738        def main_rate(t):
2739            if t > 20:
2740                msg = 'Model time exceeded.'
2741                raise Modeltime_too_late, msg
2742            else:
2743                return 3*t + 7
2744
2745        domain.forcing_terms = []
2746        R = Rainfall(domain, rate=main_rate,
2747                     polygon=[[1,1], [2,1], [2,2], [1,2]], default_rate=5.0)
2748
2749        assert num.allclose(R.exchange_area, 1)
2750
2751        domain.forcing_terms.append(R)
2752
2753        for t in domain.evolve(yieldstep=1, finaltime=25):
2754            pass
2755            #FIXME(Ole):  A test here is hard because explicit_update also
2756            # receives updates from the flux calculation.
2757
2758    def test_inflow_using_circle(self):
2759        from math import pi, cos, sin
2760
2761        a = [0.0, 0.0]
2762        b = [0.0, 2.0]
2763        c = [2.0, 0.0]
2764        d = [0.0, 4.0]
2765        e = [2.0, 2.0]
2766        f = [4.0, 0.0]
2767
2768        points = [a, b, c, d, e, f]
2769        #             bac,     bce,     ecf,     dbe
2770        vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
2771
2772        domain = Domain(points, vertices)
2773
2774        # Flat surface with 1m of water
2775        domain.set_quantity('elevation', 0)
2776        domain.set_quantity('stage', 1.0)
2777        domain.set_quantity('friction', 0)
2778
2779        Br = Reflective_boundary(domain)
2780        domain.set_boundary({'exterior': Br})
2781
2782        # Setup only one forcing term, constant inflow of 2 m^3/s
2783        # on a circle affecting triangles #0 and #1 (bac and bce)
2784        domain.forcing_terms = []
2785        domain.forcing_terms.append(Inflow(domain, rate=2.0,
2786                                           center=(1,1), radius=1))
2787
2788        domain.compute_forcing_terms()
2789
2790        assert num.allclose(domain.quantities['stage'].explicit_update[1],
2791                            2.0/pi)
2792        assert num.allclose(domain.quantities['stage'].explicit_update[0],
2793                            2.0/pi)
2794        assert num.allclose(domain.quantities['stage'].explicit_update[2:], 0)
2795
2796    def test_inflow_using_circle_function(self):
2797        from math import pi, cos, sin
2798
2799        a = [0.0, 0.0]
2800        b = [0.0, 2.0]
2801        c = [2.0, 0.0]
2802        d = [0.0, 4.0]
2803        e = [2.0, 2.0]
2804        f = [4.0, 0.0]
2805
2806        points = [a, b, c, d, e, f]
2807        #             bac,     bce,     ecf,     dbe
2808        vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
2809
2810        domain = Domain(points, vertices)
2811
2812        # Flat surface with 1m of water
2813        domain.set_quantity('elevation', 0)
2814        domain.set_quantity('stage', 1.0)
2815        domain.set_quantity('friction', 0)
2816
2817        Br = Reflective_boundary(domain)
2818        domain.set_boundary({'exterior': Br})
2819
2820        # Setup only one forcing term, time dependent inflow of 2 m^3/s
2821        # on a circle affecting triangles #0 and #1 (bac and bce)
2822        domain.forcing_terms = []
2823        domain.forcing_terms.append(Inflow(domain, rate=lambda t: 2.,
2824                                           center=(1,1), radius=1))
2825
2826        domain.compute_forcing_terms()
2827
2828        assert num.allclose(domain.quantities['stage'].explicit_update[1],
2829                            2.0/pi)
2830        assert num.allclose(domain.quantities['stage'].explicit_update[0],
2831                            2.0/pi)
2832        assert num.allclose(domain.quantities['stage'].explicit_update[2:], 0)
2833
2834    def test_inflow_catch_too_few_triangles(self):
2835        """
2836        Test that exception is thrown if no triangles are covered
2837        by the inflow area
2838        """
2839
2840        from math import pi, cos, sin
2841
2842        a = [0.0, 0.0]
2843        b = [0.0, 2.0]
2844        c = [2.0, 0.0]
2845        d = [0.0, 4.0]
2846        e = [2.0, 2.0]
2847        f = [4.0, 0.0]
2848
2849        points = [a, b, c, d, e, f]
2850        #             bac,     bce,     ecf,     dbe
2851        vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
2852
2853        domain = Domain(points, vertices)
2854
2855        # Flat surface with 1m of water
2856        domain.set_quantity('elevation', 0)
2857        domain.set_quantity('stage', 1.0)
2858        domain.set_quantity('friction', 0)
2859
2860        Br = Reflective_boundary(domain)
2861        domain.set_boundary({'exterior': Br})
2862
2863        # Setup only one forcing term, constant inflow of 2 m^3/s
2864        # on a circle affecting triangles #0 and #1 (bac and bce)
2865        try:
2866            Inflow(domain, rate=2.0, center=(1,1.1), radius=0.01)
2867        except:
2868            pass
2869        else:
2870            msg = 'Should have raised exception'
2871            raise Exception, msg
2872
2873    def Xtest_inflow_outflow_conservation(self):
2874        """
2875        Test what happens if water is abstracted from one area and
2876        injected into another - especially if there is not enough
2877        water to match the abstraction.
2878        This tests that the total volume is kept constant under a range of
2879        scenarios.
2880
2881        This test will fail as the problem was only fixed for culverts.
2882        """
2883
2884        from math import pi, cos, sin
2885
2886        length = 20.
2887        width = 10.
2888
2889        dx = dy = 2    # 1 or 2 OK
2890        points, vertices, boundary = rectangular_cross(int(length/dx),
2891                                                       int(width/dy),
2892                                                       len1=length,
2893                                                       len2=width)
2894        domain = Domain(points, vertices, boundary)
2895        domain.set_name('test_inflow_conservation')    # Output name
2896        domain.set_default_order(2)
2897
2898        # Flat surface with 1m of water
2899        stage = 1.0
2900        domain.set_quantity('elevation', 0)
2901        domain.set_quantity('stage', stage)
2902        domain.set_quantity('friction', 0)
2903
2904        Br = Reflective_boundary(domain)
2905        domain.set_boundary({'left': Br, 'right': Br, 'bottom': Br, 'top': Br})
2906
2907        # Setup one forcing term, constant inflow of 2 m^3/s on a circle
2908        domain.forcing_terms = []
2909        domain.forcing_terms.append(Inflow(domain, rate=2.0,
2910                                           center=(5,5), radius=1))
2911
2912        domain.compute_forcing_terms()
2913
2914        # Check that update values are correct
2915        for x in domain.quantities['stage'].explicit_update:
2916            assert num.allclose(x, 2.0/pi) or num.allclose(x, 0.0)
2917
2918        # Check volumes without inflow
2919        domain.forcing_terms = []
2920        initial_volume = domain.quantities['stage'].get_integral()
2921
2922        assert num.allclose(initial_volume, width*length*stage)
2923
2924        for t in domain.evolve(yieldstep = 0.05, finaltime = 5.0):
2925            volume = domain.quantities['stage'].get_integral()
2926            assert num.allclose(volume, initial_volume)
2927
2928        # Now apply the inflow and check volumes for a range of stage values
2929        for stage in [2.0, 1.0, 0.5, 0.25, 0.1, 0.0]:
2930            domain.time = 0.0
2931            domain.set_quantity('stage', stage)
2932            domain.forcing_terms = []
2933            domain.forcing_terms.append(Inflow(domain, rate=2.0,
2934                                               center=(5,5), radius=1))
2935            initial_volume = domain.quantities['stage'].get_integral()
2936            predicted_volume = initial_volume
2937            dt = 0.05
2938            for t in domain.evolve(yieldstep=dt, finaltime=5.0):
2939                volume = domain.quantities['stage'].get_integral()
2940                assert num.allclose (volume, predicted_volume)
2941                predicted_volume = predicted_volume + 2.0/pi/100/dt # Why 100?
2942
2943        # Apply equivalent outflow only and check volumes
2944        # for a range of stage values
2945        for stage in [2.0, 1.0, 0.5, 0.25, 0.1, 0.0]:
2946            print stage
2947
2948            domain.time = 0.0
2949            domain.set_quantity('stage', stage)
2950            domain.forcing_terms = []
2951            domain.forcing_terms.append(Inflow(domain, rate=-2.0,
2952                                               center=(15,5), radius=1))
2953            initial_volume = domain.quantities['stage'].get_integral()
2954            predicted_volume = initial_volume
2955            dt = 0.05
2956            for t in domain.evolve(yieldstep=dt, finaltime=5.0):
2957                volume = domain.quantities['stage'].get_integral()
2958                print t, volume, predicted_volume
2959                assert num.allclose (volume, predicted_volume)
2960                predicted_volume = predicted_volume - 2.0/pi/100/dt # Why 100?
2961
2962        # Apply both inflow and outflow and check volumes being constant for a
2963        # range of stage values
2964        for stage in [2.0, 1.0, 0.5, 0.25, 0.1, 0.0]:
2965            print stage
2966
2967            domain.time = 0.0
2968            domain.set_quantity('stage', stage)
2969            domain.forcing_terms = []
2970            domain.forcing_terms.append(Inflow(domain, rate=2.0,
2971                                               center=(5,5), radius=1))
2972            domain.forcing_terms.append(Inflow(domain, rate=-2.0,
2973                                               center=(15,5), radius=1))
2974            initial_volume = domain.quantities['stage'].get_integral()
2975
2976            dt = 0.05
2977            for t in domain.evolve(yieldstep=dt, finaltime=5.0):
2978                volume = domain.quantities['stage'].get_integral()
2979
2980                print t, volume
2981                assert num.allclose(volume, initial_volume)
2982
2983    #####################################################
2984
2985    def test_first_order_extrapolator_const_z(self):
2986        a = [0.0, 0.0]
2987        b = [0.0, 2.0]
2988        c = [2.0, 0.0]
2989        d = [0.0, 4.0]
2990        e = [2.0, 2.0]
2991        f = [4.0, 0.0]
2992
2993        points = [a, b, c, d, e, f]
2994        #             bac,     bce,     ecf,     dbe
2995        vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
2996
2997        domain = Domain(points, vertices)
2998        val0 = 2. + 2.0/3
2999        val1 = 4. + 4.0/3
3000        val2 = 8. + 2.0/3
3001        val3 = 2. + 8.0/3
3002
3003        zl = zr = -3.75    # Assume constant bed (must be less than stage)
3004        domain.set_quantity('elevation', zl*num.ones((4, 3), num.int)) #array default#
3005        domain.set_quantity('stage', [[val0, val0-1, val0-2],
3006                                      [val1, val1+1, val1],
3007                                      [val2, val2-2, val2],
3008                                      [val3-0.5, val3, val3]])
3009
3010        domain._order_ = 1
3011        domain.distribute_to_vertices_and_edges()
3012
3013        #Check that centroid values were distributed to vertices
3014        C = domain.quantities['stage'].centroid_values
3015        for i in range(3):
3016            assert num.allclose(domain.quantities['stage'].vertex_values[:,i],
3017                                C)
3018
3019    def test_first_order_limiter_variable_z(self):
3020        '''Check that first order limiter follows bed_slope'''
3021
3022        from anuga.config import epsilon
3023
3024        a = [0.0, 0.0]
3025        b = [0.0, 2.0]
3026        c = [2.0, 0.0]
3027        d = [0.0, 4.0]
3028        e = [2.0, 2.0]
3029        f = [4.0, 0.0]
3030
3031        points = [a, b, c, d, e, f]
3032        #             bac,     bce,     ecf,     dbe
3033        vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
3034
3035        domain = Domain(points, vertices)
3036        val0 = 2. + 2.0/3
3037        val1 = 4. + 4.0/3
3038        val2 = 8. + 2.0/3
3039        val3 = 2. + 8.0/3
3040
3041        domain.set_quantity('elevation', [[0,0,0], [6,0,0],
3042                                          [6,6,6], [6,6,6]])
3043        domain.set_quantity('stage', [[val0, val0, val0],
3044                                      [val1, val1, val1],
3045                                      [val2, val2, val2],
3046                                      [val3, val3, val3]])
3047
3048        E = domain.quantities['elevation'].vertex_values
3049        L = domain.quantities['stage'].vertex_values
3050
3051        # Check that some stages are not above elevation (within eps) -
3052        # so that the limiter has something to work with
3053        assert not num.alltrue(num.alltrue(num.greater_equal(L, E-epsilon)))
3054
3055        domain._order_ = 1
3056        domain.distribute_to_vertices_and_edges()
3057
3058        #Check that all stages are above elevation (within eps)
3059        assert num.alltrue(num.alltrue(num.greater_equal(L, E-epsilon)))
3060
3061    #####################################################
3062
3063    def test_distribute_basic(self):
3064        #Using test data generated by abstract_2d_finite_volumes-2
3065        #Assuming no friction and flat bed (0.0)
3066
3067        a = [0.0, 0.0]
3068        b = [0.0, 2.0]
3069        c = [2.0, 0.0]
3070        d = [0.0, 4.0]
3071        e = [2.0, 2.0]
3072        f = [4.0, 0.0]
3073
3074        points = [a, b, c, d, e, f]
3075        #             bac,     bce,     ecf,     dbe
3076        vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
3077
3078        domain = Domain(points, vertices)
3079
3080        val0 = 2.
3081        val1 = 4.
3082        val2 = 8.
3083        val3 = 2.
3084
3085        domain.set_quantity('stage', [val0, val1, val2, val3],
3086                            location='centroids')
3087        L = domain.quantities['stage'].vertex_values
3088
3089        # First order
3090        domain._order_ = 1
3091        domain.distribute_to_vertices_and_edges()
3092        assert num.allclose(L[1], val1)
3093
3094        # Second order
3095        domain._order_ = 2
3096        domain.beta_w = 0.9
3097        domain.beta_w_dry = 0.9
3098        domain.beta_uh = 0.9
3099        domain.beta_uh_dry = 0.9
3100        domain.beta_vh = 0.9
3101        domain.beta_vh_dry = 0.9
3102        domain.distribute_to_vertices_and_edges()
3103        assert num.allclose(L[1], [2.2, 4.9, 4.9])
3104
3105    def test_distribute_away_from_bed(self):
3106        #Using test data generated by abstract_2d_finite_volumes-2
3107        #Assuming no friction and flat bed (0.0)
3108
3109        a = [0.0, 0.0]
3110        b = [0.0, 2.0]
3111        c = [2.0, 0.0]
3112        d = [0.0, 4.0]
3113        e = [2.0, 2.0]
3114        f = [4.0, 0.0]
3115
3116        points = [a, b, c, d, e, f]
3117        #             bac,     bce,     ecf,     dbe
3118        vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
3119
3120        domain = Domain(points, vertices)
3121        L = domain.quantities['stage'].vertex_values
3122
3123        def stage(x, y):
3124            return x**2
3125
3126        domain.set_quantity('stage', stage, location='centroids')
3127
3128        domain.quantities['stage'].compute_gradients()
3129
3130        a, b = domain.quantities['stage'].get_gradients()
3131
3132        assert num.allclose(a[1], 3.33333334)
3133        assert num.allclose(b[1], 0.0)
3134
3135        domain._order_ = 1
3136        domain.distribute_to_vertices_and_edges()
3137        assert num.allclose(L[1], 1.77777778)
3138
3139        domain._order_ = 2
3140        domain.beta_w = 0.9
3141        domain.beta_w_dry = 0.9
3142        domain.beta_uh = 0.9
3143        domain.beta_uh_dry = 0.9
3144        domain.beta_vh = 0.9
3145        domain.beta_vh_dry = 0.9
3146        domain.distribute_to_vertices_and_edges()
3147        assert num.allclose(L[1], [0.57777777, 2.37777778, 2.37777778])
3148
3149    def test_distribute_away_from_bed1(self):
3150        #Using test data generated by abstract_2d_finite_volumes-2
3151        #Assuming no friction and flat bed (0.0)
3152
3153        a = [0.0, 0.0]
3154        b = [0.0, 2.0]
3155        c = [2.0, 0.0]
3156        d = [0.0, 4.0]
3157        e = [2.0, 2.0]
3158        f = [4.0, 0.0]
3159
3160        points = [a, b, c, d, e, f]
3161        #             bac,     bce,     ecf,     dbe
3162        vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
3163
3164        domain = Domain(points, vertices)
3165        L = domain.quantities['stage'].vertex_values
3166
3167        def stage(x, y):
3168            return x**4 + y**2
3169
3170        domain.set_quantity('stage', stage, location='centroids')
3171
3172        domain.quantities['stage'].compute_gradients()
3173        a, b = domain.quantities['stage'].get_gradients()
3174        assert num.allclose(a[1], 25.18518519)
3175        assert num.allclose(b[1], 3.33333333)
3176
3177        domain._order_ = 1
3178        domain.distribute_to_vertices_and_edges()
3179        assert num.allclose(L[1], 4.9382716)
3180
3181        domain._order_ = 2
3182        domain.beta_w = 0.9
3183        domain.beta_w_dry = 0.9
3184        domain.beta_uh = 0.9
3185        domain.beta_uh_dry = 0.9
3186        domain.beta_vh = 0.9
3187        domain.beta_vh_dry = 0.9
3188        domain.distribute_to_vertices_and_edges()
3189        assert num.allclose(L[1], [1.07160494, 6.46058131, 7.28262855])
3190
3191    def test_distribute_near_bed(self):
3192        a = [0.0, 0.0]
3193        b = [0.0, 2.0]
3194        c = [2.0, 0.0]
3195        d = [0.0, 4.0]
3196        e = [2.0, 2.0]
3197        f = [4.0, 0.0]
3198
3199        points = [a, b, c, d, e, f]
3200        #             bac,     bce,     ecf,     dbe
3201        vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
3202
3203        domain = Domain(points, vertices)
3204
3205        # Set up for a gradient of (10,0) at mid triangle (bce)
3206        def slope(x, y):
3207            return 10*x
3208
3209        h = 0.1
3210        def stage(x, y):
3211            return slope(x, y) + h
3212
3213        domain.set_quantity('elevation', slope)
3214        domain.set_quantity('stage', stage, location='centroids')
3215
3216        E = domain.quantities['elevation'].vertex_values
3217        L = domain.quantities['stage'].vertex_values
3218
3219        # Get reference values
3220        volumes = []
3221        for i in range(len(L)):
3222            volumes.append(num.sum(L[i])/3)
3223            assert num.allclose(volumes[i],
3224                                domain.quantities['stage'].centroid_values[i])
3225
3226        domain._order_ = 1
3227
3228        domain.tight_slope_limiters = 0
3229        domain.distribute_to_vertices_and_edges()
3230        assert num.allclose(L[1], [0.1, 20.1, 20.1])
3231        for i in range(len(L)):
3232            assert num.allclose(volumes[i], num.sum(L[i])/3)
3233
3234        # Allow triangle to be flatter (closer to bed)
3235        domain.tight_slope_limiters = 1
3236
3237        domain.distribute_to_vertices_and_edges()
3238        assert num.allclose(L[1], [0.298, 20.001, 20.001])
3239        for i in range(len(L)):
3240            assert num.allclose(volumes[i], num.sum(L[i])/3)
3241
3242        domain._order_ = 2
3243
3244        domain.tight_slope_limiters = 0
3245        domain.distribute_to_vertices_and_edges()
3246        assert num.allclose(L[1], [0.1, 20.1, 20.1])
3247        for i in range(len(L)):
3248            assert num.allclose(volumes[i], num.sum(L[i])/3)
3249
3250        # Allow triangle to be flatter (closer to bed)
3251        domain.tight_slope_limiters = 1
3252
3253        domain.distribute_to_vertices_and_edges()
3254        assert num.allclose(L[1], [0.298, 20.001, 20.001])
3255        for i in range(len(L)):
3256            assert num.allclose(volumes[i], num.sum(L[i])/3)
3257
3258    def test_distribute_near_bed1(self):
3259        a = [0.0, 0.0]
3260        b = [0.0, 2.0]
3261        c = [2.0, 0.0]
3262        d = [0.0, 4.0]
3263        e = [2.0, 2.0]
3264        f = [4.0, 0.0]
3265
3266        points = [a, b, c, d, e, f]
3267        #             bac,     bce,     ecf,     dbe
3268        vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
3269
3270        domain = Domain(points, vertices)
3271
3272        # Set up for a gradient of (8,2) at mid triangle (bce)
3273        def slope(x, y):
3274            return x**4 + y**2
3275
3276        h = 0.1
3277        def stage(x, y):
3278            return slope(x, y) + h
3279
3280        domain.set_quantity('elevation', slope)
3281        domain.set_quantity('stage', stage)
3282
3283        E = domain.quantities['elevation'].vertex_values
3284        L = domain.quantities['stage'].vertex_values
3285
3286        # Get reference values
3287        volumes = []
3288        for i in range(len(L)):
3289            volumes.append(num.sum(L[i])/3)
3290            assert num.allclose(volumes[i],
3291                                domain.quantities['stage'].centroid_values[i])
3292
3293        domain._order_ = 1
3294
3295        domain.tight_slope_limiters = 0
3296        domain.distribute_to_vertices_and_edges()
3297        assert num.allclose(L[1], [4.1, 16.1, 20.1])
3298        for i in range(len(L)):
3299            assert num.allclose(volumes[i], num.sum(L[i])/3)
3300
3301        # Allow triangle to be flatter (closer to bed)
3302        domain.tight_slope_limiters = 1
3303
3304        domain.distribute_to_vertices_and_edges()
3305        assert num.allclose(L[1], [4.2386, 16.0604, 20.001])
3306        for i in range(len(L)):
3307            assert num.allclose(volumes[i], num.sum(L[i])/3)
3308
3309        domain._order_ = 2
3310
3311        domain.tight_slope_limiters = 0
3312        domain.distribute_to_vertices_and_edges()
3313        assert num.allclose(L[1], [4.1, 16.1, 20.1])
3314        for i in range(len(L)):
3315            assert num.allclose(volumes[i], num.sum(L[i])/3)
3316
3317        # Allow triangle to be flatter (closer to bed)
3318        domain.tight_slope_limiters = 1
3319
3320        domain.distribute_to_vertices_and_edges()
3321        assert (num.allclose(L[1], [4.23370103, 16.06529897, 20.001]) or
3322                num.allclose(L[1], [4.18944138, 16.10955862, 20.001]) or
3323                num.allclose(L[1], [4.19351461, 16.10548539, 20.001]))
3324        # old limiters
3325
3326        for i in range(len(L)):
3327            assert num.allclose(volumes[i], num.sum(L[i])/3)
3328
3329    def test_second_order_distribute_real_data(self):
3330        #Using test data generated by abstract_2d_finite_volumes-2
3331        #Assuming no friction and flat bed (0.0)
3332
3333        a = [0.0, 0.0]
3334        b = [0.0, 1.0/5]
3335        c = [0.0, 2.0/5]
3336        d = [1.0/5, 0.0]
3337        e = [1.0/5, 1.0/5]
3338        f = [1.0/5, 2.0/5]
3339        g = [2.0/5, 2.0/5]
3340
3341        points = [a, b, c, d, e, f, g]
3342        #             bae,     efb,     cbf,     feg
3343        vertices = [[1,0,4], [4,5,1], [2,1,5], [5,4,6]]
3344
3345        domain = Domain(points, vertices)
3346
3347        def slope(x, y):
3348            return -x/3
3349
3350        domain.set_quantity('elevation', slope)
3351        domain.set_quantity('stage',
3352                            [0.01298164, 0.00365611,
3353                             0.01440365, -0.0381856437096],
3354                            location='centroids')
3355        domain.set_quantity('xmomentum',
3356                            [0.00670439, 0.01263789,
3357                             0.00647805, 0.0178180740668],
3358                            location='centroids')
3359        domain.set_quantity('ymomentum',
3360                            [-7.23510980e-004, -6.30413883e-005,
3361                             6.30413883e-005, 0.000200907255866],
3362                            location='centroids')
3363
3364        E = domain.quantities['elevation'].vertex_values
3365        L = domain.quantities['stage'].vertex_values
3366        X = domain.quantities['xmomentum'].vertex_values
3367        Y = domain.quantities['ymomentum'].vertex_values
3368
3369        domain._order_ = 2
3370        domain.beta_w = 0.9
3371        domain.beta_w_dry = 0.9
3372        domain.beta_uh = 0.9
3373        domain.beta_uh_dry = 0.9
3374        domain.beta_vh = 0.9
3375        domain.beta_vh_dry = 0.9
3376
3377        # FIXME (Ole): Need tests where this is commented out
3378        domain.tight_slope_limiters = 0       # Backwards compatibility (14/4/7)
3379        domain.use_centroid_velocities = 0    # Backwards compatibility (7/5/8)
3380
3381        domain.distribute_to_vertices_and_edges()
3382
3383        assert num.allclose(L[1,:], [-0.00825735775384,
3384                                     -0.00801881482869,
3385                                     0.0272445025825])
3386        assert num.allclose(X[1,:], [0.0143507718962,
3387                                     0.0142502147066,
3388                                     0.00931268339717])
3389        assert num.allclose(Y[1,:], [-0.000117062180693,
3390                                     7.94434448109e-005,
3391                                     -0.000151505429018])
3392
3393    def test_balance_deep_and_shallow(self):
3394        """Test that balanced limiters preserve conserved quantites.
3395        This test is using old depth based balanced limiters
3396        """
3397
3398        import copy
3399
3400        a = [0.0, 0.0]
3401        b = [0.0, 2.0]
3402        c = [2.0, 0.0]
3403        d = [0.0, 4.0]
3404        e = [2.0, 2.0]
3405        f = [4.0, 0.0]
3406
3407        points = [a, b, c, d, e, f]
3408        #             bac,     bce,     ecf,     dbe
3409        elements = [[1,0,2], [1,2,4], [4,2,5], [3,1,4] ]
3410
3411        domain = Domain(points, elements)
3412        domain.check_integrity()
3413
3414        # Create a deliberate overshoot
3415        domain.set_quantity('stage', [[3,0,3], [2,2,6], [5,3,8], [8,3,5]])
3416        domain.set_quantity('elevation', 0)    # Flat bed
3417        stage = domain.quantities['stage']
3418
3419        ref_centroid_values = copy.copy(stage.centroid_values[:])    # Copy
3420
3421        # Limit
3422        domain.tight_slope_limiters = 0
3423        domain.distribute_to_vertices_and_edges()
3424
3425        # Assert that quantities are conserved
3426        for k in range(len(domain)):
3427            assert num.allclose(ref_centroid_values[k],
3428                                num.sum(stage.vertex_values[k,:])/3)
3429
3430        # Now try with a non-flat bed - closely hugging initial stage in places
3431        # This will create alphas in the range [0, 0.478260, 1]
3432        domain.set_quantity('stage', [[3,0,3], [2,2,6], [5,3,8], [8,3,5]])
3433        domain.set_quantity('elevation', [[0,0,0],
3434                                          [1.8,1.9,5.9],
3435                                          [4.6,0,0],
3436                                          [0,2,4]])
3437        stage = domain.quantities['stage']
3438
3439        ref_centroid_values = copy.copy(stage.centroid_values[:])    # Copy
3440        ref_vertex_values = copy.copy(stage.vertex_values[:])        # Copy
3441
3442        # Limit
3443        domain.tight_slope_limiters = 0
3444        domain.distribute_to_vertices_and_edges()
3445
3446        # Assert that all vertex quantities have changed
3447        for k in range(len(domain)):
3448            assert not num.allclose(ref_vertex_values[k,:],
3449                                    stage.vertex_values[k,:])
3450        # and assert that quantities are still conserved
3451        for k in range(len(domain)):
3452            assert num.allclose(ref_centroid_values[k],
3453                                num.sum(stage.vertex_values[k,:])/3)
3454
3455        # Check actual results
3456        assert num.allclose(stage.vertex_values,
3457                            [[2,2,2],
3458                             [1.93333333, 2.03333333, 6.03333333],
3459                             [6.93333333, 4.53333333, 4.53333333],
3460                             [5.33333333, 5.33333333, 5.33333333]]) 
3461
3462    def test_balance_deep_and_shallow_tight_SL(self):
3463        """Test that balanced limiters preserve conserved quantites.
3464        This test is using Tight Slope Limiters
3465        """
3466
3467        import copy
3468
3469        a = [0.0, 0.0]
3470        b = [0.0, 2.0]
3471        c = [2.0, 0.0]
3472        d = [0.0, 4.0]
3473        e = [2.0, 2.0]
3474        f = [4.0, 0.0]
3475
3476        points = [a, b, c, d, e, f]
3477        #             bac,     bce,     ecf,     dbe
3478        elements = [[1,0,2], [1,2,4], [4,2,5], [3,1,4] ]
3479
3480        domain = Domain(points, elements)
3481        domain.check_integrity()
3482
3483        # Create a deliberate overshoot
3484        domain.set_quantity('stage', [[3,0,3], [2,2,6], [5,3,8], [8,3,5]])
3485        domain.set_quantity('elevation', 0)    # Flat bed
3486        stage = domain.quantities['stage']
3487
3488        ref_centroid_values = copy.copy(stage.centroid_values[:])    # Copy
3489
3490        # Limit
3491        domain.tight_slope_limiters = 1
3492        domain.distribute_to_vertices_and_edges()
3493
3494        # Assert that quantities are conserved
3495        for k in range(len(domain)):
3496            assert num.allclose (ref_centroid_values[k],
3497                                 num.sum(stage.vertex_values[k,:])/3)
3498
3499        # Now try with a non-flat bed - closely hugging initial stage in places
3500        # This will create alphas in the range [0, 0.478260, 1]
3501        domain.set_quantity('stage', [[3,0,3], [2,2,6], [5,3,8], [8,3,5]])
3502        domain.set_quantity('elevation', [[0,0,0],
3503                                          [1.8,1.9,5.9],
3504                                          [4.6,0,0],
3505                                          [0,2,4]])
3506        stage = domain.quantities['stage']
3507
3508        ref_centroid_values = copy.copy(stage.centroid_values[:])    # Copy
3509        ref_vertex_values = copy.copy(stage.vertex_values[:])        # Copy
3510
3511        # Limit
3512        domain.tight_slope_limiters = 1
3513        domain.distribute_to_vertices_and_edges()
3514
3515        # Assert that all vertex quantities have changed
3516        for k in range(len(domain)):
3517            assert not num.allclose(ref_vertex_values[k,:],
3518                                    stage.vertex_values[k,:])
3519        # and assert that quantities are still conserved
3520        for k in range(len(domain)):
3521            assert num.allclose(ref_centroid_values[k],
3522                                num.sum(stage.vertex_values[k,:])/3)
3523
3524    def test_balance_deep_and_shallow_Froude(self):
3525        """Test that balanced limiters preserve conserved quantites -
3526        and also that excessive Froude numbers are dealt with.
3527        This test is using tight slope limiters.
3528        """
3529
3530        import copy
3531
3532        a = [0.0, 0.0]
3533        b = [0.0, 2.0]
3534        c = [2.0, 0.0]
3535        d = [0.0, 4.0]
3536        e = [2.0, 2.0]
3537        f = [4.0, 0.0]
3538
3539        points = [a, b, c, d, e, f]
3540        #             bac,     bce,     ecf,     dbe
3541        elements = [[1,0,2], [1,2,4], [4,2,5], [3,1,4] ]
3542
3543        domain = Domain(points, elements)
3544        domain.check_integrity()
3545        domain.tight_slope_limiters = True
3546        domain.use_centroid_velocities = True
3547
3548        # Create non-flat bed - closely hugging initial stage in places
3549        # This will create alphas in the range [0, 0.478260, 1]
3550        domain.set_quantity('stage', [[3,0,3], [2,2,6], [5,3,8], [8,3,5]])
3551        domain.set_quantity('elevation', [[0,0,0],
3552                                          [1.8,1.999,5.999],
3553                                          [4.6,0,0],
3554                                          [0,2,4]])
3555
3556        # Create small momenta, that nonetheless will generate large speeds
3557        # due to shallow depth at isolated vertices
3558        domain.set_quantity('xmomentum', -0.0058)
3559        domain.set_quantity('ymomentum', 0.0890)
3560
3561        stage = domain.quantities['stage']
3562        elevation = domain.quantities['elevation']
3563        xmomentum = domain.quantities['xmomentum']
3564        ymomentum = domain.quantities['ymomentum']
3565
3566        # Setup triangle #1 to mimick real Froude explosion observed
3567        # in the Onslow example 13 Nov 2007.
3568        stage.vertex_values[1,:] = [1.6385, 1.6361, 1.2953]
3569        elevation.vertex_values[1,:] = [1.6375, 1.6336, 0.4647]
3570        xmomentum.vertex_values[1,:] = [-0.0058, -0.0050, -0.0066]
3571        ymomentum.vertex_values[1,:] = [0.0890, 0.0890, 0.0890]
3572
3573        xmomentum.interpolate()
3574        ymomentum.interpolate()
3575        stage.interpolate()
3576        elevation.interpolate()
3577
3578        # Verify interpolation
3579        assert num.allclose(stage.centroid_values[1], 1.5233)
3580        assert num.allclose(elevation.centroid_values[1], 1.2452667)
3581        assert num.allclose(xmomentum.centroid_values[1], -0.0058)
3582        assert num.allclose(ymomentum.centroid_values[1], 0.089)
3583
3584        # Derived quantities
3585        depth = stage-elevation
3586        u = xmomentum/depth
3587        v = ymomentum/depth
3588
3589        denom = (depth*g)**0.5
3590        Fx = u/denom
3591        Fy = v/denom
3592
3593        # Verify against Onslow example (14 Nov 2007)
3594        assert num.allclose(depth.centroid_values[1], 0.278033)
3595        assert num.allclose(u.centroid_values[1], -0.0208608)
3596        assert num.allclose(v.centroid_values[1], 0.3201055)
3597
3598        assert num.allclose(denom.centroid_values[1],
3599                            num.sqrt(depth.centroid_values[1]*g))
3600
3601        assert num.allclose(u.centroid_values[1]/denom.centroid_values[1],
3602                            -0.012637746977)
3603        assert num.allclose(Fx.centroid_values[1],
3604                            u.centroid_values[1]/denom.centroid_values[1])
3605
3606        # Check that Froude numbers are small at centroids.
3607        assert num.allclose(Fx.centroid_values[1], -0.012637746977)
3608        assert num.allclose(Fy.centroid_values[1], 0.193924048435)
3609
3610        # But Froude numbers are huge at some vertices and edges
3611        assert num.allclose(Fx.vertex_values[1,:], [-5.85888475e+01,
3612                                                    -1.27775313e+01,
3613                                                    -2.78511420e-03])
3614
3615        assert num.allclose(Fx.edge_values[1,:], [-6.89150773e-03,
3616                                                  -7.38672488e-03,
3617                                                  -2.35626238e+01])
3618
3619        assert num.allclose(Fy.vertex_values[1,:], [8.99035764e+02,
3620                                                    2.27440057e+02,
3621                                                    3.75568430e-02])
3622
3623        assert num.allclose(Fy.edge_values[1,:], [1.05748998e-01,
3624                                                  1.06035244e-01,
3625                                                  3.88346947e+02])
3626
3627
3628        # The task is now to arrange the limiters such that Froude numbers
3629        # remain under control whil at the same time obeying the conservation
3630        # laws.
3631        ref_centroid_values = copy.copy(stage.centroid_values[:])    # Copy
3632        ref_vertex_values = copy.copy(stage.vertex_values[:])        # Copy
3633
3634        # Limit (and invoke balance_deep_and_shallow)
3635        domain.tight_slope_limiters = 1
3636        domain.distribute_to_vertices_and_edges()
3637
3638        # Redo derived quantities
3639        depth = stage - elevation
3640        u = xmomentum/depth
3641        v = ymomentum/depth
3642
3643        # Assert that all vertex velocities stay within one
3644        # order of magnitude of centroid velocities.
3645        assert num.alltrue(num.absolute(u.vertex_values[1,:]) <=
3646                           num.absolute(u.centroid_values[1])*10)
3647        assert num.alltrue(num.absolute(v.vertex_values[1,:]) <=
3648                           num.absolute(v.centroid_values[1])*10)
3649
3650        denom = (depth*g)**0.5
3651        Fx = u/denom
3652        Fy = v/denom
3653
3654        # Assert that Froude numbers are less than max value (TBA)
3655        # at vertices, edges and centroids.
3656        from anuga.config import maximum_froude_number
3657
3658        assert num.alltrue(num.absolute(Fx.vertex_values[1,:]) <
3659                           maximum_froude_number)
3660        assert num.alltrue(num.absolute(Fy.vertex_values[1,:]) <
3661                           maximum_froude_number)
3662
3663        # Assert that all vertex quantities have changed
3664        for k in range(len(domain)):
3665            assert not num.allclose(ref_vertex_values[k,:],
3666                                    stage.vertex_values[k,:])
3667
3668        # Assert that quantities are still conserved
3669        for k in range(len(domain)):
3670            assert num.allclose(ref_centroid_values[k],
3671                                num.sum(stage.vertex_values[k,:])/3)
3672
3673        return
3674
3675        qwidth = 12
3676        for k in [1]:    # range(len(domain)):
3677            print 'Triangle %d (C, V, E)' % k
3678
3679            print ('stage'.ljust(qwidth), stage.centroid_values[k],
3680                   stage.vertex_values[k,:], stage.edge_values[k,:])
3681            print ('elevation'.ljust(qwidth), elevation.centroid_values[k],
3682                   elevation.vertex_values[k,:], elevation.edge_values[k,:])
3683            print ('depth'.ljust(qwidth), depth.centroid_values[k],
3684                   depth.vertex_values[k,:], depth.edge_values[k,:])
3685            print ('xmomentum'.ljust(qwidth), xmomentum.centroid_values[k],
3686                   xmomentum.vertex_values[k,:], xmomentum.edge_values[k,:])
3687            print ('ymomentum'.ljust(qwidth), ymomentum.centroid_values[k],
3688                   ymomentum.vertex_values[k,:], ymomentum.edge_values[k,:])
3689            print ('u'.ljust(qwidth),u.centroid_values[k],
3690                   u.vertex_values[k,:], u.edge_values[k,:])
3691            print ('v'.ljust(qwidth), v.centroid_values[k],
3692                   v.vertex_values[k,:], v.edge_values[k,:])
3693            print ('Fx'.ljust(qwidth), Fx.centroid_values[k],
3694                   Fx.vertex_values[k,:], Fx.edge_values[k,:])
3695            print ('Fy'.ljust(qwidth), Fy.centroid_values[k],
3696                   Fy.vertex_values[k,:], Fy.edge_values[k,:])
3697
3698    def test_conservation_1(self):
3699        """Test that stage is conserved globally
3700
3701        This one uses a flat bed, reflective bdries and a suitable
3702        initial condition
3703        """
3704
3705        from mesh_factory import rectangular
3706
3707        # Create basic mesh
3708        points, vertices, boundary = rectangular(6, 6)
3709
3710        # Create shallow water domain
3711        domain = Domain(points, vertices, boundary)
3712        domain.smooth = False
3713        domain.default_order = 2
3714
3715        # IC
3716        def x_slope(x, y):
3717            return x/3
3718
3719        domain.set_quantity('elevation', 0)
3720        domain.set_quantity('friction', 0)
3721        domain.set_quantity('stage', x_slope)
3722
3723        # Boundary conditions (reflective everywhere)
3724        Br = Reflective_boundary(domain)
3725        domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
3726
3727        domain.check_integrity()
3728
3729        initial_volume = domain.quantities['stage'].get_integral()
3730        initial_xmom = domain.quantities['xmomentum'].get_integral()
3731
3732        # Evolution
3733        for t in domain.evolve(yieldstep=0.05, finaltime=5.0):
3734            volume = domain.quantities['stage'].get_integral()
3735            assert num.allclose(volume, initial_volume)
3736
3737            #I don't believe that the total momentum should be the same
3738            #It starts with zero and ends with zero though
3739            #xmom = domain.quantities['xmomentum'].get_integral()
3740            #print xmom
3741            #assert allclose (xmom, initial_xmom)
3742
3743        os.remove(domain.get_name() + '.sww')
3744
3745    def test_conservation_2(self):
3746        """Test that stage is conserved globally
3747
3748        This one uses a slopy bed, reflective bdries and a suitable
3749        initial condition
3750        """
3751
3752        from mesh_factory import rectangular
3753
3754        # Create basic mesh
3755        points, vertices, boundary = rectangular(6, 6)
3756
3757        # Create shallow water domain
3758        domain = Domain(points, vertices, boundary)
3759        domain.smooth = False
3760        domain.default_order = 2
3761
3762        # IC
3763        def x_slope(x, y):
3764            return x/3
3765
3766        domain.set_quantity('elevation', x_slope)
3767        domain.set_quantity('friction', 0)
3768        domain.set_quantity('stage', 0.4)    # Steady
3769
3770        # Boundary conditions (reflective everywhere)
3771        Br = Reflective_boundary(domain)
3772        domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
3773
3774        domain.check_integrity()
3775
3776        initial_volume = domain.quantities['stage'].get_integral()
3777        initial_xmom = domain.quantities['xmomentum'].get_integral()
3778
3779        # Evolution
3780        for t in domain.evolve(yieldstep=0.05, finaltime=5.0):
3781            volume = domain.quantities['stage'].get_integral()
3782            assert num.allclose(volume, initial_volume)
3783
3784            #FIXME: What would we expect from momentum
3785            #xmom = domain.quantities['xmomentum'].get_integral()
3786            #print xmom
3787            #assert allclose (xmom, initial_xmom)
3788
3789        os.remove(domain.get_name() + '.sww')
3790
3791    def test_conservation_3(self):
3792        """Test that stage is conserved globally
3793
3794        This one uses a larger grid, convoluted bed, reflective boundaries
3795        and a suitable initial condition
3796        """
3797
3798        from mesh_factory import rectangular
3799
3800        # Create basic mesh
3801        points, vertices, boundary = rectangular(2, 1)
3802
3803        # Create shallow water domain
3804        domain = Domain(points, vertices, boundary)
3805        domain.smooth = False
3806        domain.default_order = 2
3807        domain.set_quantities_to_be_stored(['stage', 'xmomentum', 'ymomentum'])
3808
3809        # IC
3810        def x_slope(x, y):
3811            z = 0*x
3812            for i in range(len(x)):
3813                if x[i] < 0.3:
3814                    z[i] = x[i]/3
3815                if 0.3 <= x[i] < 0.5:
3816                    z[i] = -0.5
3817                if 0.5 <= x[i] < 0.7:
3818                    z[i] = 0.39
3819                if 0.7 <= x[i]:
3820                    z[i] = x[i]/3
3821            return z
3822
3823        domain.set_quantity('elevation', x_slope)
3824        domain.set_quantity('friction', 0)
3825        domain.set_quantity('stage', 0.4) #Steady
3826
3827        # Boundary conditions (reflective everywhere)
3828        Br = Reflective_boundary(domain)
3829        domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
3830
3831        domain.check_integrity()
3832
3833        initial_volume = domain.quantities['stage'].get_integral()
3834        initial_xmom = domain.quantities['xmomentum'].get_integral()
3835
3836        import copy
3837
3838        ref_centroid_values = copy.copy(domain.quantities['stage'].\
3839                                            centroid_values)
3840
3841        domain.distribute_to_vertices_and_edges()
3842
3843        assert num.allclose(domain.quantities['stage'].centroid_values,
3844                            ref_centroid_values)
3845
3846        # Check that initial limiter doesn't violate cons quan
3847        assert num.allclose(domain.quantities['stage'].get_integral(),
3848                            initial_volume)
3849
3850        # Evolution
3851        for t in domain.evolve(yieldstep=0.05, finaltime=10):
3852            volume =  domain.quantities['stage'].get_integral()
3853            assert num.allclose (volume, initial_volume)
3854
3855        os.remove(domain.get_name() + '.sww')
3856
3857    def test_conservation_4(self):
3858        """Test that stage is conserved globally
3859
3860        This one uses a larger grid, convoluted bed, reflective boundaries
3861        and a suitable initial condition
3862        """
3863
3864        from mesh_factory import rectangular
3865
3866        # Create basic mesh
3867        points, vertices, boundary = rectangular(6, 6)
3868
3869        # Create shallow water domain
3870        domain = Domain(points, vertices, boundary)
3871        domain.smooth = False
3872        domain.default_order = 2
3873
3874        # IC
3875        def x_slope(x, y):
3876            z = 0*x
3877            for i in range(len(x)):
3878                if x[i] < 0.3:
3879                    z[i] = x[i]/3
3880                if 0.3 <= x[i] < 0.5:
3881                    z[i] = -0.5
3882                if 0.5 <= x[i] < 0.7:
3883                    #z[i] = 0.3     # OK with beta == 0.2
3884                    z[i] = 0.34     # OK with beta == 0.0
3885                    #z[i] = 0.35    # Fails after 80 timesteps with an error
3886                                    # of the order 1.0e-5
3887                if 0.7 <= x[i]:
3888                    z[i] = x[i]/3
3889            return z
3890
3891        domain.set_quantity('elevation', x_slope)
3892        domain.set_quantity('friction', 0)
3893        domain.set_quantity('stage', 0.4) #Steady
3894
3895        # Boundary conditions (reflective everywhere)
3896        Br = Reflective_boundary(domain)
3897        domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
3898
3899        domain.check_integrity()
3900
3901        initial_volume = domain.quantities['stage'].get_integral()
3902        initial_xmom = domain.quantities['xmomentum'].get_integral()
3903
3904        import copy
3905
3906        ref_centroid_values = copy.copy(domain.quantities['stage'].\
3907                                            centroid_values)
3908
3909        # Test limiter by itself
3910        domain.distribute_to_vertices_and_edges()
3911
3912        # Check that initial limiter doesn't violate cons quan
3913        assert num.allclose(domain.quantities['stage'].get_integral(),
3914                            initial_volume)
3915        # NOTE: This would fail if any initial stage was less than the
3916        # corresponding bed elevation - but that is reasonable.
3917
3918        #Evolution
3919        for t in domain.evolve(yieldstep=0.05, finaltime=10.0):
3920            volume =  domain.quantities['stage'].get_integral()
3921            assert num.allclose (volume, initial_volume)
3922
3923        os.remove(domain.get_name() + '.sww')
3924
3925    def test_conservation_5(self):
3926        """Test that momentum is conserved globally in steady state scenario
3927
3928        This one uses a slopy bed, dirichlet and reflective bdries
3929        """
3930
3931        from mesh_factory import rectangular
3932
3933        # Create basic mesh
3934        points, vertices, boundary = rectangular(6, 6)
3935
3936        # Create shallow water domain
3937        domain = Domain(points, vertices, boundary)
3938        domain.smooth = False
3939        domain.default_order = 2
3940
3941        # IC
3942        def x_slope(x, y):
3943            return x/3
3944
3945        domain.set_quantity('elevation', x_slope)
3946        domain.set_quantity('friction', 0)
3947        domain.set_quantity('stage', 0.4) # Steady
3948
3949        # Boundary conditions (reflective everywhere)
3950        Br = Reflective_boundary(domain)
3951        Bleft = Dirichlet_boundary([0.5, 0, 0])
3952        Bright = Dirichlet_boundary([0.1, 0, 0])
3953        domain.set_boundary({'left': Bleft, 'right': Bright,
3954                             'top': Br, 'bottom': Br})
3955
3956        domain.check_integrity()
3957
3958        initial_volume = domain.quantities['stage'].get_integral()
3959        initial_xmom = domain.quantities['xmomentum'].get_integral()
3960
3961        # Evolution
3962        for t in domain.evolve(yieldstep=0.05, finaltime=15.0):
3963            stage =  domain.quantities['stage'].get_integral()
3964            xmom = domain.quantities['xmomentum'].get_integral()
3965            ymom = domain.quantities['ymomentum'].get_integral()
3966
3967            if num.allclose(t, 6):    # Steady state reached
3968                steady_xmom = domain.quantities['xmomentum'].get_integral()
3969                steady_ymom = domain.quantities['ymomentum'].get_integral()
3970                steady_stage = domain.quantities['stage'].get_integral()
3971
3972            if t > 6:
3973                msg = 'xmom=%.2f, steady_xmom=%.2f' % (xmom, steady_xmom)
3974                assert num.allclose(xmom, steady_xmom), msg
3975                assert num.allclose(ymom, steady_ymom)
3976                assert num.allclose(stage, steady_stage)
3977
3978        os.remove(domain.get_name() + '.sww')
3979
3980    def test_conservation_real(self):
3981        """Test that momentum is conserved globally
3982
3983        Stephen finally made a test that revealed the problem.
3984        This test failed with code prior to 25 July 2005
3985        """
3986
3987        import sys
3988        import os.path
3989        sys.path.append(os.path.join('..', 'abstract_2d_finite_volumes'))
3990        from mesh_factory import rectangular
3991
3992        yieldstep = 0.01
3993        finaltime = 0.05
3994        min_depth = 1.0e-2
3995
3996        #Create shallow water domain
3997        points, vertices, boundary = rectangular(10, 10, len1=500, len2=500)
3998        domain = Domain(points, vertices, boundary)
3999        domain.smooth = False
4000        domain.default_order = 1
4001        domain.minimum_allowed_height = min_depth
4002
4003        # Set initial condition
4004        class Set_IC:
4005            """Set an initial condition with a constant value, for x0<x<x1"""
4006
4007            def __init__(self, x0=0.25, x1=0.5, h=1.0):
4008                self.x0 = x0
4009                self.x1 = x1
4010                self.= h
4011
4012            def __call__(self, x, y):
4013                return self.h*((x > self.x0) & (x < self.x1))
4014
4015        domain.set_quantity('stage', Set_IC(200.0, 300.0, 5.0))
4016
4017        # Boundaries
4018        R = Reflective_boundary(domain)
4019        domain.set_boundary({'left': R, 'right': R, 'top':R, 'bottom': R})
4020
4021        ref = domain.quantities['stage'].get_integral()
4022
4023        # Evolution
4024        for t in domain.evolve(yieldstep=yieldstep, finaltime=finaltime):
4025            pass
4026
4027        now = domain.quantities['stage'].get_integral()
4028
4029        msg = 'Stage not conserved: was %f, now %f' % (ref, now)
4030        assert num.allclose(ref, now), msg
4031
4032        os.remove(domain.get_name() + '.sww')
4033
4034    def test_second_order_flat_bed_onestep(self):
4035        from mesh_factory import rectangular
4036
4037        #Create basic mesh
4038        points, vertices, boundary = rectangular(6, 6)
4039
4040        #Create shallow water domain
4041        domain = Domain(points, vertices, boundary)
4042        domain.smooth = False
4043        domain.default_order = 2
4044        domain.beta_w = 0.9
4045        domain.beta_w_dry = 0.9
4046        domain.beta_uh = 0.9
4047        domain.beta_uh_dry = 0.9
4048        domain.beta_vh = 0.9
4049        domain.beta_vh_dry = 0.9
4050        domain.H0 = 0    # Backwards compatibility (6/2/7)
4051
4052        # Boundary conditions
4053        Br = Reflective_boundary(domain)
4054        Bd = Dirichlet_boundary([0.1, 0., 0.])
4055        domain.set_boundary({'left': Bd, 'right': Br, 'top': Br, 'bottom': Br})
4056
4057        domain.check_integrity()
4058
4059        # Evolution
4060        for t in domain.evolve(yieldstep=0.05, finaltime=0.05):
4061            pass
4062
4063        # Data from earlier version of abstract_2d_finite_volumes
4064        assert num.allclose(domain.min_timestep, 0.0396825396825)
4065        assert num.allclose(domain.max_timestep, 0.0396825396825)
4066
4067        assert num.allclose(domain.quantities['stage'].centroid_values[:12],
4068                            [0.00171396, 0.02656103, 0.00241523, 0.02656103,
4069                             0.00241523, 0.02656103, 0.00241523, 0.02656103,
4070                             0.00241523, 0.02656103, 0.00241523, 0.0272623])
4071
4072        domain.distribute_to_vertices_and_edges()
4073
4074        assert num.allclose(domain.quantities['stage'].vertex_values[:12,0],
4075                            [0.0001714, 0.02656103, 0.00024152,
4076                             0.02656103, 0.00024152, 0.02656103,
4077                             0.00024152, 0.02656103, 0.00024152,
4078                             0.02656103, 0.00024152, 0.0272623])
4079
4080        assert num.allclose(domain.quantities['stage'].vertex_values[:12,1],
4081                            [0.00315012, 0.02656103, 0.00024152, 0.02656103,
4082                             0.00024152, 0.02656103, 0.00024152, 0.02656103,
4083                             0.00024152, 0.02656103, 0.00040506, 0.0272623])
4084
4085        assert num.allclose(domain.quantities['stage'].vertex_values[:12,2],
4086                            [0.00182037, 0.02656103, 0.00676264,
4087                             0.02656103, 0.00676264, 0.02656103,
4088                             0.00676264, 0.02656103, 0.00676264,
4089                             0.02656103, 0.0065991,  0.0272623])
4090
4091        assert num.allclose(domain.quantities['xmomentum'].centroid_values[:12],
4092                            [0.00113961, 0.01302432, 0.00148672,
4093                             0.01302432, 0.00148672, 0.01302432,
4094                             0.00148672, 0.01302432, 0.00148672 ,
4095                             0.01302432, 0.00148672, 0.01337143])
4096
4097        assert num.allclose(domain.quantities['ymomentum'].centroid_values[:12],
4098                            [-2.91240050e-004 , 1.22721531e-004,
4099                             -1.22721531e-004,  1.22721531e-004 ,
4100                             -1.22721531e-004,  1.22721531e-004,
4101                             -1.22721531e-004 , 1.22721531e-004,
4102                             -1.22721531e-004,  1.22721531e-004,
4103                             -1.22721531e-004,  -4.57969873e-005])
4104
4105        os.remove(domain.get_name() + '.sww')
4106
4107    def test_second_order_flat_bed_moresteps(self):
4108        from mesh_factory import rectangular
4109
4110        # Create basic mesh
4111        points, vertices, boundary = rectangular(6, 6)
4112
4113        # Create shallow water domain
4114        domain = Domain(points, vertices, boundary)
4115        domain.smooth = False
4116        domain.default_order = 2
4117
4118        # Boundary conditions
4119        Br = Reflective_boundary(domain)
4120        Bd = Dirichlet_boundary([0.1, 0., 0.])
4121        domain.set_boundary({'left': Bd, 'right': Br, 'top': Br, 'bottom': Br})
4122
4123        domain.check_integrity()
4124
4125        # Evolution
4126        for t in domain.evolve(yieldstep=0.05, finaltime=0.1):
4127            pass
4128
4129        # Data from earlier version of abstract_2d_finite_volumes
4130        #assert allclose(domain.min_timestep, 0.0396825396825)
4131        #assert allclose(domain.max_timestep, 0.0396825396825)
4132        #print domain.quantities['stage'].centroid_values
4133
4134        os.remove(domain.get_name() + '.sww')
4135
4136    def test_flatbed_first_order(self):
4137        from mesh_factory import rectangular
4138
4139        # Create basic mesh
4140        N = 8
4141        points, vertices, boundary = rectangular(N, N)
4142
4143        # Create shallow water domain
4144        domain = Domain(points, vertices, boundary)
4145        domain.smooth = False
4146        domain.default_order = 1
4147        domain.H0 = 0    # Backwards compatibility (6/2/7)
4148
4149        # Boundary conditions
4150        Br = Reflective_boundary(domain)
4151        Bd = Dirichlet_boundary([0.2, 0., 0.])
4152
4153        domain.set_boundary({'left': Bd, 'right': Br, 'top': Br, 'bottom': Br})
4154        domain.check_integrity()
4155
4156        # Evolution
4157        for t in domain.evolve(yieldstep=0.02, finaltime=0.5):
4158            pass
4159
4160        # FIXME: These numbers were from version before 25/10
4161        #assert allclose(domain.min_timestep, 0.0140413643926)
4162        #assert allclose(domain.max_timestep, 0.0140947355753)
4163
4164        for i in range(3):
4165            #assert allclose(domain.quantities['stage'].edge_values[:4,i],
4166            #                [0.10730244,0.12337617,0.11200126,0.12605666])
4167            assert num.allclose(domain.quantities['xmomentum'].\
4168                                    edge_values[:4,i],
4169                                [0.07610894,0.06901572,0.07284461,0.06819712])
4170            #assert allclose(domain.quantities['ymomentum'].edge_values[:4,i],
4171            #                [-0.0060238, -0.00157404, -0.00309633, -0.0001637])
4172
4173        os.remove(domain.get_name() + '.sww')
4174
4175    def test_flatbed_second_order(self):
4176        from mesh_factory import rectangular
4177
4178        # Create basic mesh
4179        N = 8
4180        points, vertices, boundary = rectangular(N, N)
4181
4182        # Create shallow water domain
4183        domain = Domain(points, vertices, boundary)
4184        domain.smooth = False
4185        domain.default_order = 2
4186        domain.beta_w = 0.9
4187        domain.beta_w_dry = 0.9
4188        domain.beta_uh = 0.9
4189        domain.beta_uh_dry = 0.9
4190        domain.beta_vh = 0.9
4191        domain.beta_vh_dry = 0.9
4192        #Makes it like the 'oldstyle' balance
4193        #domain.minimum_allowed_height = 0.0
4194        domain.H0 = 0    # Backwards compatibility (6/2/7)
4195        domain.use_centroid_velocities = False # Backwards compatibility (8/5/8)
4196        domain.set_maximum_allowed_speed(1.0)
4197
4198        # Boundary conditions
4199        Br = Reflective_boundary(domain)
4200        Bd = Dirichlet_boundary([0.2, 0., 0.])
4201
4202        domain.set_boundary({'left': Bd, 'right': Br, 'top': Br, 'bottom': Br})
4203        domain.check_integrity()
4204
4205        # Evolution
4206        for t in domain.evolve(yieldstep=0.01, finaltime=0.03):
4207            pass
4208
4209        msg = 'min step was %f instead of %f' % (domain.min_timestep,
4210                                                 0.0210448446782)
4211
4212        assert num.allclose(domain.min_timestep, 0.0210448446782), msg
4213        assert num.allclose(domain.max_timestep, 0.0210448446782)
4214
4215        #FIXME: These numbers were from version before 25/10
4216        #assert allclose(domain.quantities['stage'].vertex_values[:4,0],
4217        #                [0.00101913,0.05352143,0.00104852,0.05354394])
4218
4219        #FIXME: These numbers were from version before 21/3/6 -
4220        #could be recreated by setting maximum_allowed_speed to 0 maybe
4221        #assert allclose(domain.quantities['xmomentum'].vertex_values[:4,0],
4222        #                [ 0.00064835, 0.03685719, 0.00085073, 0.03687313])
4223
4224        assert num.allclose(domain.quantities['xmomentum'].vertex_values[:4,0],
4225                            [ 0.00090581, 0.03685719, 0.00088303, 0.03687313])
4226
4227        #assert allclose(domain.quantities['xmomentum'].vertex_values[:4,0],
4228        #                [0.00090581,0.03685719,0.00088303,0.03687313])
4229
4230        assert num.allclose(domain.quantities['ymomentum'].vertex_values[:4,0],
4231                            [-0.00139463,0.0006156,-0.00060364,0.00061827])
4232
4233        os.remove(domain.get_name() + '.sww')
4234
4235    def test_flatbed_second_order_vmax_0(self):
4236        from mesh_factory import rectangular
4237
4238        # Create basic mesh
4239        N = 8
4240        points, vertices, boundary = rectangular(N, N)
4241
4242        # Create shallow water domain
4243        domain = Domain(points, vertices, boundary)
4244        domain.smooth = False
4245        domain.default_order = 2
4246        domain.beta_w = 0.9
4247        domain.beta_w_dry = 0.9
4248        domain.beta_uh = 0.9
4249        domain.beta_uh_dry = 0.9
4250        domain.beta_vh = 0.9
4251        domain.beta_vh_dry = 0.9
4252        domain.maximum_allowed_speed = 0.0    # Makes it like the 'oldstyle'
4253        domain.H0 = 0    # Backwards compatibility (6/2/7)
4254        domain.use_centroid_velocities = False # Backwards compatibility (8/5/8)
4255
4256        # Boundary conditions
4257        Br = Reflective_boundary(domain)
4258        Bd = Dirichlet_boundary([0.2, 0., 0.])
4259
4260        domain.set_boundary({'left': Bd, 'right': Br, 'top': Br, 'bottom': Br})
4261        domain.check_integrity()
4262
4263        # Evolution
4264        for t in domain.evolve(yieldstep=0.01, finaltime=0.03):
4265            pass
4266
4267        assert num.allclose(domain.min_timestep, 0.0210448446782)
4268        assert num.allclose(domain.max_timestep, 0.0210448446782)
4269
4270        #FIXME: These numbers were from version before 21/3/6 -
4271        #could be recreated by setting maximum_allowed_speed to 0 maybe
4272        assert num.allclose(domain.quantities['xmomentum'].vertex_values[:4,0],
4273                            [ 0.00064835, 0.03685719, 0.00085073, 0.03687313])
4274
4275
4276        assert num.allclose(domain.quantities['ymomentum'].vertex_values[:4,0],
4277                            [-0.00139463,0.0006156,-0.00060364,0.00061827])
4278
4279        os.remove(domain.get_name() + '.sww')
4280
4281    def test_flatbed_second_order_distribute(self):
4282        #Use real data from anuga.abstract_2d_finite_volumes 2
4283        #painfully setup and extracted.
4284
4285        from mesh_factory import rectangular
4286
4287        # Create basic mesh
4288        N = 8
4289        points, vertices, boundary = rectangular(N, N)
4290
4291        # Create shallow water domain
4292        domain = Domain(points, vertices, boundary)
4293        domain.smooth = False
4294        domain.default_order=domain._order_ = 2
4295        domain.beta_w = 0.9
4296        domain.beta_w_dry = 0.9
4297        domain.beta_uh = 0.9
4298        domain.beta_uh_dry = 0.9
4299        domain.beta_vh = 0.9
4300        domain.beta_vh_dry = 0.9
4301        domain.H0 = 0    # Backwards compatibility (6/2/7)
4302        domain.use_centroid_velocities = False # Backwards compatibility (8/5/8)
4303        domain.set_maximum_allowed_speed(1.0)
4304
4305        # Boundary conditions
4306        Br = Reflective_boundary(domain)
4307        Bd = Dirichlet_boundary([0.2, 0., 0.])
4308
4309        domain.set_boundary({'left': Bd, 'right': Br, 'top': Br, 'bottom': Br})
4310        domain.check_integrity()
4311
4312        for V in [False, True]:
4313            if V:
4314                # Set centroids as if system had been evolved
4315                L = num.zeros(2*N*N, num.float)
4316                L[:32] = [7.21205592e-003, 5.35214298e-002, 1.00910824e-002,
4317                          5.35439433e-002, 1.00910824e-002, 5.35439433e-002,
4318                          1.00910824e-002, 5.35439433e-002, 1.00910824e-002,
4319                          5.35439433e-002, 1.00910824e-002, 5.35439433e-002,
4320                          1.00910824e-002, 5.35393928e-002, 1.02344264e-002,
4321                          5.59605058e-002, 0.00000000e+000, 3.31027800e-004,
4322                          0.00000000e+000, 4.37962142e-005, 0.00000000e+000,
4323                          4.37962142e-005, 0.00000000e+000, 4.37962142e-005,
4324                          0.00000000e+000, 4.37962142e-005, 0.00000000e+000,
4325                          4.37962142e-005, 0.00000000e+000, 4.37962142e-005,
4326                          0.00000000e+000, 5.57305948e-005]
4327
4328                X = num.zeros(2*N*N, num.float)
4329                X[:32] = [6.48351607e-003, 3.68571894e-002, 8.50733285e-003,
4330                          3.68731327e-002, 8.50733285e-003, 3.68731327e-002,
4331                          8.50733285e-003, 3.68731327e-002, 8.50733285e-003,
4332                          3.68731327e-002, 8.50733285e-003, 3.68731327e-002,
4333                          8.50733285e-003, 3.68693861e-002, 8.65220973e-003,
4334                          3.85055387e-002, 0.00000000e+000, 2.86060840e-004,
4335                          0.00000000e+000, 3.58905503e-005, 0.00000000e+000,
4336                          3.58905503e-005, 0.00000000e+000, 3.58905503e-005,
4337                          0.00000000e+000, 3.58905503e-005, 0.00000000e+000,
4338                          3.58905503e-005, 0.00000000e+000, 3.58905503e-005,
4339                          0.00000000e+000, 4.57662812e-005]
4340
4341                Y = num.zeros(2*N*N, num.float)
4342                Y[:32] = [-1.39463104e-003, 6.15600298e-004, -6.03637382e-004,
4343                          6.18272251e-004, -6.03637382e-004, 6.18272251e-004,
4344                          -6.03637382e-004, 6.18272251e-004, -6.03637382e-004,
4345                          6.18272251e-004, -6.03637382e-004, 6.18272251e-004,
4346                          -6.03637382e-004, 6.18599320e-004, -6.74622797e-004,
4347                          -1.48934756e-004, 0.00000000e+000, -5.35079969e-005,
4348                          0.00000000e+000, -2.57264987e-005, 0.00000000e+000,
4349                          -2.57264987e-005, 0.00000000e+000, -2.57264987e-005,
4350                          0.00000000e+000, -2.57264987e-005, 0.00000000e+000,
4351                          -2.57264987e-005, 0.00000000e+000, -2.57264987e-005,
4352                          0.00000000e+000, -2.57635178e-005]
4353
4354                domain.set_quantity('stage', L, location='centroids')
4355                domain.set_quantity('xmomentum', X, location='centroids')
4356                domain.set_quantity('ymomentum', Y, location='centroids')
4357
4358                domain.check_integrity()
4359            else:
4360                # Evolution
4361                for t in domain.evolve(yieldstep=0.01, finaltime=0.03):
4362                    pass
4363                assert num.allclose(domain.min_timestep, 0.0210448446782)
4364                assert num.allclose(domain.max_timestep, 0.0210448446782)
4365
4366            # Centroids were correct but not vertices.
4367            # Hence the check of distribute below.
4368            assert num.allclose(domain.quantities['stage'].centroid_values[:4],
4369                                [0.00721206, 0.05352143,
4370                                 0.01009108, 0.05354394])
4371
4372            assert num.allclose(domain.quantities['xmomentum'].\
4373                                    centroid_values[:4],
4374                                [0.00648352, 0.03685719,
4375                                 0.00850733, 0.03687313])
4376
4377            assert num.allclose(domain.quantities['ymomentum'].\
4378                                    centroid_values[:4],
4379                                [-0.00139463, 0.0006156,
4380                                 -0.00060364, 0.00061827])
4381
4382            #assert allclose(domain.quantities['xmomentum'].\
4383            #                    centroid_values[17],0.00028606084)
4384            if not V:
4385                #FIXME: These numbers were from version before 21/3/6 -
4386                #could be recreated by setting maximum_allowed_speed to 0 maybe
4387                #assert allclose(domain.quantities['xmomentum'].\
4388                #                    centroid_values[17], 0.0)
4389                assert num.allclose(domain.quantities['xmomentum'].\
4390                                        centroid_values[17],
4391                                    0.000286060839592)
4392
4393            else:
4394                assert num.allclose(domain.quantities['xmomentum'].\
4395                                        centroid_values[17],
4396                                    0.00028606084)
4397
4398            import copy
4399
4400            XX = copy.copy(domain.quantities['xmomentum'].centroid_values)
4401            assert num.allclose(domain.quantities['xmomentum'].centroid_values,
4402                                XX)
4403
4404            domain.distribute_to_vertices_and_edges()
4405
4406            #assert allclose(domain.quantities['xmomentum'].centroid_values, XX)
4407            #assert allclose(domain.quantities['xmomentum'].centroid_values[17],
4408            #                0.0)
4409
4410            assert num.allclose(domain.quantities['xmomentum'].\
4411                                    centroid_values[17],
4412                                0.000286060839592)
4413
4414            #FIXME: These numbers were from version before 25/10
4415            #assert allclose(domain.quantities['stage'].vertex_values[:4,0],
4416            #                [0.00101913,0.05352143,0.00104852,0.05354394])
4417
4418            assert num.allclose(domain.quantities['ymomentum'].\
4419                                    vertex_values[:4,0],
4420                                [-0.00139463, 0.0006156,
4421                                 -0.00060364, 0.00061827])
4422
4423            assert num.allclose(domain.quantities['xmomentum'].\
4424                                    vertex_values[:4,0],
4425                                [0.00090581, 0.03685719,
4426                                 0.00088303, 0.03687313])
4427
4428        os.remove(domain.get_name() + '.sww')
4429
4430    def test_bedslope_problem_first_order(self):
4431        from mesh_factory import rectangular
4432
4433        # Create basic mesh
4434        points, vertices, boundary = rectangular(6, 6)
4435
4436        # Create shallow water domain
4437        domain = Domain(points, vertices, boundary)
4438        domain.smooth = False
4439        domain.default_order = 1
4440
4441        # Bed-slope and friction
4442        def x_slope(x, y):
4443            return -x/3
4444
4445        domain.set_quantity('elevation', x_slope)
4446
4447        # Boundary conditions
4448        Br = Reflective_boundary(domain)
4449        domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
4450
4451        # Initial condition
4452        domain.set_quantity('stage', expression='elevation+0.05')
4453        domain.check_integrity()
4454
4455        # Evolution
4456        for t in domain.evolve(yieldstep=0.05, finaltime=0.05):
4457            pass
4458
4459        # FIXME (Ole): Need some other assertion here!
4460        #print domain.min_timestep, domain.max_timestep
4461        #assert allclose(domain.min_timestep, 0.050010003001)
4462        #assert allclose(domain.max_timestep, 0.050010003001)
4463
4464        os.remove(domain.get_name() + '.sww')
4465
4466    def test_bedslope_problem_first_order_moresteps(self):
4467        from mesh_factory import rectangular
4468
4469        # Create basic mesh
4470        points, vertices, boundary = rectangular(6, 6)
4471
4472        # Create shallow water domain
4473        domain = Domain(points, vertices, boundary)
4474        domain.smooth = False
4475        domain.default_order = 1
4476
4477        # FIXME (Ole): Need tests where these two are commented out
4478        domain.H0 = 0                         # Backwards compatibility (6/2/7)
4479        domain.tight_slope_limiters = 0       # Backwards compatibility (14/4/7)
4480        domain.use_centroid_velocities = 0    # Backwards compatibility (7/5/8)
4481
4482        # Bed-slope and friction
4483        def x_slope(x, y):
4484            return -x/3
4485
4486        domain.set_quantity('elevation', x_slope)
4487
4488        # Boundary conditions
4489        Br = Reflective_boundary(domain)
4490        domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
4491
4492        # Initial condition
4493        domain.set_quantity('stage', expression='elevation+0.05')
4494        domain.check_integrity()
4495
4496        # Evolution
4497        for t in domain.evolve(yieldstep=0.05, finaltime=0.5):
4498            pass
4499
4500        #Data from earlier version of abstract_2d_finite_volumes
4501        assert num.allclose(domain.quantities['stage'].centroid_values,
4502                            [-0.02998628, -0.01520652, -0.03043492,
4503                             -0.0149132,  -0.03004706, -0.01476251,
4504                             -0.0298215,  -0.01467976, -0.02988158,
4505                             -0.01474662, -0.03036161, -0.01442995,
4506                             -0.07624583, -0.06297061, -0.07733792,
4507                             -0.06342237, -0.07695439, -0.06289595,
4508                             -0.07635559, -0.0626065,  -0.07633628,
4509                             -0.06280072, -0.07739632, -0.06386738,
4510                             -0.12161738, -0.11028239, -0.1223796,
4511                             -0.11095953, -0.12189744, -0.11048616,
4512                             -0.12074535, -0.10987605, -0.12014311,
4513                             -0.10976691, -0.12096859, -0.11087692,
4514                             -0.16868259, -0.15868061, -0.16801135,
4515                             -0.1588003,  -0.16674343, -0.15813323,
4516                             -0.16457595, -0.15693826, -0.16281096,
4517                             -0.15585154, -0.16283873, -0.15540068,
4518                             -0.17450362, -0.19919913, -0.18062882,
4519                             -0.19764131, -0.17783111, -0.19407213,
4520                             -0.1736915,  -0.19053624, -0.17228678,
4521                             -0.19105634, -0.17920133, -0.1968828,
4522                             -0.14244395, -0.14604641, -0.14473537,
4523                             -0.1506107,  -0.14510055, -0.14919522,
4524                             -0.14175896, -0.14560798, -0.13911658,
4525                             -0.14439383, -0.13924047, -0.14829043])
4526
4527        os.remove(domain.get_name() + '.sww')
4528
4529    def test_bedslope_problem_second_order_one_step(self):
4530        from mesh_factory import rectangular
4531
4532        # Create basic mesh
4533        points, vertices, boundary = rectangular(6, 6)
4534
4535        # Create shallow water domain
4536        domain = Domain(points, vertices, boundary)
4537        domain.smooth = False
4538        domain.default_order = 2
4539        domain.beta_w = 0.9
4540        domain.beta_w_dry = 0.9
4541        domain.beta_uh = 0.9
4542        domain.beta_uh_dry = 0.9
4543        domain.beta_vh = 0.9
4544        domain.beta_vh_dry = 0.9
4545
4546
4547        # FIXME (Ole): Need tests where this is commented out
4548        domain.tight_slope_limiters = 0       # Backwards compatibility (14/4/7)
4549        domain.use_centroid_velocities = 0    # Backwards compatibility (7/5/8)
4550
4551        # Bed-slope and friction at vertices (and interpolated elsewhere)
4552        def x_slope(x, y):
4553            return -x/3
4554
4555        domain.set_quantity('elevation', x_slope)
4556
4557        # Boundary conditions
4558        Br = Reflective_boundary(domain)
4559        domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
4560
4561        #Initial condition
4562        domain.set_quantity('stage', expression='elevation+0.05')
4563        domain.check_integrity()
4564
4565        assert num.allclose(domain.quantities['stage'].centroid_values,
4566                            [ 0.01296296,  0.03148148,  0.01296296,
4567                              0.03148148,  0.01296296,  0.03148148,
4568                              0.01296296,  0.03148148,  0.01296296,
4569                              0.03148148,  0.01296296,  0.03148148,
4570                             -0.04259259, -0.02407407, -0.04259259,
4571                             -0.02407407, -0.04259259, -0.02407407,
4572                             -0.04259259, -0.02407407, -0.04259259,
4573                             -0.02407407, -0.04259259, -0.02407407,
4574                             -0.09814815, -0.07962963, -0.09814815,
4575                             -0.07962963, -0.09814815, -0.07962963,
4576                             -0.09814815, -0.07962963, -0.09814815,
4577                             -0.07962963, -0.09814815, -0.07962963,
4578                             -0.1537037,  -0.13518519, -0.1537037,
4579                             -0.13518519, -0.1537037,  -0.13518519,
4580                             -0.1537037 , -0.13518519, -0.1537037,
4581                             -0.13518519, -0.1537037,  -0.13518519,
4582                             -0.20925926, -0.19074074, -0.20925926,
4583                             -0.19074074, -0.20925926, -0.19074074,
4584                             -0.20925926, -0.19074074, -0.20925926,
4585                             -0.19074074, -0.20925926, -0.19074074,
4586                             -0.26481481, -0.2462963,  -0.26481481,
4587                             -0.2462963,  -0.26481481, -0.2462963,
4588                             -0.26481481, -0.2462963,  -0.26481481,
4589                             -0.2462963,  -0.26481481, -0.2462963])
4590
4591
4592        # Evolution
4593        for t in domain.evolve(yieldstep=0.05, finaltime=0.05):
4594            pass
4595
4596        assert num.allclose(domain.quantities['stage'].centroid_values,
4597                            [ 0.01290985,  0.02356019,  0.01619096,
4598                              0.02356019,  0.01619096,  0.02356019,
4599                              0.01619096,  0.02356019,  0.01619096,
4600                              0.02356019,  0.01619096,  0.0268413,
4601                             -0.04411074, -0.0248011,  -0.04186556,
4602                             -0.0248011,  -0.04186556, -0.0248011,
4603                             -0.04186556, -0.0248011,  -0.04186556,
4604                             -0.0248011,  -0.04186556, -0.02255593,
4605                             -0.09966629, -0.08035666, -0.09742112,
4606                             -0.08035666, -0.09742112, -0.08035666,
4607                             -0.09742112, -0.08035666, -0.09742112,
4608                             -0.08035666, -0.09742112, -0.07811149,
4609                             -0.15522185, -0.13591222, -0.15297667,
4610                             -0.13591222, -0.15297667, -0.13591222,
4611                             -0.15297667, -0.13591222, -0.15297667,
4612                             -0.13591222, -0.15297667, -0.13366704,
4613                             -0.2107774,  -0.19146777, -0.20853223,
4614                             -0.19146777, -0.20853223, -0.19146777,
4615                             -0.20853223, -0.19146777, -0.20853223,
4616                             -0.19146777, -0.20853223, -0.1892226,
4617                             -0.26120669, -0.24776246, -0.25865535,
4618                             -0.24776246, -0.25865535, -0.24776246,
4619                             -0.25865535, -0.24776246, -0.25865535,
4620                             -0.24776246, -0.25865535, -0.24521113])
4621
4622        os.remove(domain.get_name() + '.sww')
4623
4624    def test_bedslope_problem_second_order_two_steps(self):
4625        from mesh_factory import rectangular
4626
4627        # Create basic mesh
4628        points, vertices, boundary = rectangular(6, 6)
4629
4630        # Create shallow water domain
4631        domain = Domain(points, vertices, boundary)
4632        domain.smooth = False
4633        domain.default_order = 2
4634        domain.beta_w = 0.9
4635        domain.beta_w_dry = 0.9
4636        domain.beta_uh = 0.9
4637        domain.beta_uh_dry = 0.9
4638        domain.beta_vh = 0.9
4639        domain.beta_vh_dry = 0.9
4640
4641        # FIXME (Ole): Need tests where this is commented out
4642        domain.tight_slope_limiters = 0    # Backwards compatibility (14/4/7)
4643        domain.H0 = 0    # Backwards compatibility (6/2/7)
4644        domain.use_centroid_velocities = 0    # Backwards compatibility (7/5/8)
4645
4646        # Bed-slope and friction at vertices (and interpolated elsewhere)
4647        def x_slope(x, y):
4648            return -x/3
4649
4650        domain.set_quantity('elevation', x_slope)
4651
4652        # Boundary conditions
4653        Br = Reflective_boundary(domain)
4654        domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
4655
4656        # Initial condition
4657        domain.set_quantity('stage', expression='elevation+0.05')
4658        domain.check_integrity()
4659
4660        assert num.allclose(domain.quantities['stage'].centroid_values,
4661                            [ 0.01296296,  0.03148148,  0.01296296,
4662                              0.03148148,  0.01296296,  0.03148148,
4663                              0.01296296,  0.03148148,  0.01296296,
4664                              0.03148148,  0.01296296,  0.03148148,
4665                             -0.04259259, -0.02407407, -0.04259259,
4666                             -0.02407407, -0.04259259, -0.02407407,
4667                             -0.04259259, -0.02407407, -0.04259259,
4668                             -0.02407407, -0.04259259, -0.02407407,
4669                             -0.09814815, -0.07962963, -0.09814815,
4670                             -0.07962963, -0.09814815, -0.07962963,
4671                             -0.09814815, -0.07962963, -0.09814815,
4672                             -0.07962963, -0.09814815, -0.07962963,
4673                             -0.1537037 , -0.13518519, -0.1537037,
4674                             -0.13518519, -0.1537037,  -0.13518519,
4675                             -0.1537037 , -0.13518519, -0.1537037,
4676                             -0.13518519, -0.1537037,  -0.13518519,
4677                             -0.20925926, -0.19074074, -0.20925926,
4678                             -0.19074074, -0.20925926, -0.19074074,
4679                             -0.20925926, -0.19074074, -0.20925926,
4680                             -0.19074074, -0.20925926, -0.19074074,
4681                             -0.26481481, -0.2462963,  -0.26481481,
4682                             -0.2462963,  -0.26481481, -0.2462963,
4683                             -0.26481481, -0.2462963,  -0.26481481,
4684                             -0.2462963,  -0.26481481, -0.2462963])
4685
4686        # Evolution
4687        for t in domain.evolve(yieldstep=0.05, finaltime=0.1):
4688            pass
4689
4690        # Data from earlier version of abstract_2d_finite_volumes ft=0.1
4691        assert num.allclose(domain.min_timestep, 0.0376895634803)
4692        assert num.allclose(domain.max_timestep, 0.0415635655309)
4693
4694        assert num.allclose(domain.quantities['stage'].centroid_values,
4695                            [ 0.00855788,  0.01575204,  0.00994606,
4696                              0.01717072,  0.01005985,  0.01716362,
4697                              0.01005985,  0.01716299,  0.01007098,
4698                              0.01736248,  0.01216452,  0.02026776,
4699                             -0.04462374, -0.02479045, -0.04199789,
4700                             -0.0229465,  -0.04184033, -0.02295693,
4701                             -0.04184013, -0.02295675, -0.04184486,
4702                             -0.0228168,  -0.04028876, -0.02036486,
4703                             -0.10029444, -0.08170809, -0.09772846,
4704                             -0.08021704, -0.09760006, -0.08022143,
4705                             -0.09759984, -0.08022124, -0.09760261,
4706                             -0.08008893, -0.09603914, -0.07758209,
4707                             -0.15584152, -0.13723138, -0.15327266,
4708                             -0.13572906, -0.15314427, -0.13573349,
4709                             -0.15314405, -0.13573331, -0.15314679,
4710                             -0.13560104, -0.15158523, -0.13310701,
4711                             -0.21208605, -0.19283913, -0.20955631,
4712                             -0.19134189, -0.20942821, -0.19134598,
4713                             -0.20942799, -0.1913458,  -0.20943005,
4714                             -0.19120952, -0.20781177, -0.18869401,
4715                             -0.25384082, -0.2463294,  -0.25047649,
4716                             -0.24464654, -0.25031159, -0.24464253,
4717                             -0.25031112, -0.24464253, -0.25031463,
4718                             -0.24454764, -0.24885323, -0.24286438])
4719
4720        os.remove(domain.get_name() + '.sww')
4721
4722    def test_bedslope_problem_second_order_two_yieldsteps(self):
4723        from mesh_factory import rectangular
4724
4725        #Create basic mesh
4726        points, vertices, boundary = rectangular(6, 6)
4727
4728        #Create shallow water domain
4729        domain = Domain(points, vertices, boundary)
4730        domain.smooth = False
4731        domain.default_order = 2
4732        domain.beta_w = 0.9
4733        domain.beta_w_dry = 0.9
4734        domain.beta_uh = 0.9
4735        domain.beta_uh_dry = 0.9
4736        domain.beta_vh = 0.9
4737        domain.beta_vh_dry = 0.9
4738
4739        # FIXME (Ole): Need tests where this is commented out
4740        domain.tight_slope_limiters = 0    # Backwards compatibility (14/4/7)
4741        domain.H0 = 0    # Backwards compatibility (6/2/7)
4742        domain.use_centroid_velocities = 0    # Backwards compatibility (7/5/8)
4743
4744
4745        # Bed-slope and friction at vertices (and interpolated elsewhere)
4746        def x_slope(x, y):
4747            return -x/3
4748
4749        domain.set_quantity('elevation', x_slope)
4750
4751        # Boundary conditions
4752        Br = Reflective_boundary(domain)
4753        domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
4754
4755        # Initial condition
4756        domain.set_quantity('stage', expression='elevation+0.05')
4757        domain.check_integrity()
4758
4759        assert num.allclose(domain.quantities['stage'].centroid_values,
4760                            [ 0.01296296,  0.03148148,  0.01296296,
4761                              0.03148148,  0.01296296,  0.03148148,
4762                              0.01296296,  0.03148148,  0.01296296,
4763                              0.03148148,  0.01296296,  0.03148148,
4764                             -0.04259259, -0.02407407, -0.04259259,
4765                             -0.02407407, -0.04259259, -0.02407407,
4766                             -0.04259259, -0.02407407, -0.04259259,
4767                             -0.02407407, -0.04259259, -0.02407407,
4768                             -0.09814815, -0.07962963, -0.09814815,
4769                             -0.07962963, -0.09814815, -0.07962963,
4770                             -0.09814815, -0.07962963, -0.09814815,
4771                             -0.07962963, -0.09814815, -0.07962963,
4772                             -0.1537037 , -0.13518519, -0.1537037,
4773                             -0.13518519, -0.1537037,  -0.13518519,
4774                             -0.1537037 , -0.13518519, -0.1537037,
4775                             -0.13518519, -0.1537037,  -0.13518519,
4776                             -0.20925926, -0.19074074, -0.20925926,
4777                             -0.19074074, -0.20925926, -0.19074074,
4778                             -0.20925926, -0.19074074, -0.20925926,
4779                             -0.19074074, -0.20925926, -0.19074074,
4780                             -0.26481481, -0.2462963,  -0.26481481,
4781                             -0.2462963,  -0.26481481, -0.2462963,
4782                             -0.26481481, -0.2462963,  -0.26481481,
4783                             -0.2462963,  -0.26481481, -0.2462963])
4784
4785
4786        # Evolution
4787        for t in domain.evolve(yieldstep=0.05, finaltime=0.1):   #0.05??
4788            pass
4789
4790        assert num.allclose(domain.quantities['stage'].centroid_values,
4791                            [ 0.00855788,  0.01575204,  0.00994606,
4792                              0.01717072,  0.01005985,  0.01716362,
4793                              0.01005985,  0.01716299,  0.01007098,
4794                              0.01736248,  0.01216452,  0.02026776,
4795                             -0.04462374, -0.02479045, -0.04199789,
4796                             -0.0229465,  -0.04184033, -0.02295693,
4797                             -0.04184013, -0.02295675, -0.04184486,
4798                             -0.0228168,  -0.04028876, -0.02036486,
4799                             -0.10029444, -0.08170809, -0.09772846,
4800                             -0.08021704, -0.09760006, -0.08022143,
4801                             -0.09759984, -0.08022124, -0.09760261,
4802                             -0.08008893, -0.09603914, -0.07758209,
4803                             -0.15584152, -0.13723138, -0.15327266,
4804                             -0.13572906, -0.15314427, -0.13573349,
4805                             -0.15314405, -0.13573331, -0.15314679,
4806                             -0.13560104, -0.15158523, -0.13310701,
4807                             -0.21208605, -0.19283913, -0.20955631,
4808                             -0.19134189, -0.20942821, -0.19134598,
4809                             -0.20942799, -0.1913458,  -0.20943005,
4810                             -0.19120952, -0.20781177, -0.18869401,
4811                             -0.25384082, -0.2463294,  -0.25047649,
4812                             -0.24464654, -0.25031159, -0.24464253,
4813                             -0.25031112, -0.24464253, -0.25031463,
4814                             -0.24454764, -0.24885323, -0.24286438])
4815
4816        os.remove(domain.get_name() + '.sww')
4817
4818    def test_bedslope_problem_second_order_more_steps(self):
4819        from mesh_factory import rectangular
4820
4821        # Create basic mesh
4822        points, vertices, boundary = rectangular(6, 6)
4823
4824        # Create shallow water domain
4825        domain = Domain(points, vertices, boundary)
4826        domain.smooth = False
4827        domain.default_order = 2
4828        domain.beta_w = 0.9
4829        domain.beta_w_dry = 0.9
4830        domain.beta_uh = 0.9
4831        domain.beta_uh_dry = 0.9
4832        domain.beta_vh = 0.9
4833        domain.beta_vh_dry = 0.9
4834
4835        # FIXME (Ole): Need tests where these two are commented out
4836        domain.H0 = 0                      # Backwards compatibility (6/2/7)
4837        domain.tight_slope_limiters = 0    # Backwards compatibility (14/4/7)
4838        domain.use_centroid_velocities = 0    # Backwards compatibility (7/5/8)
4839
4840        # Bed-slope and friction at vertices (and interpolated elsewhere)
4841        def x_slope(x, y):
4842            return -x/3
4843
4844        domain.set_quantity('elevation', x_slope)
4845
4846        # Boundary conditions
4847        Br = Reflective_boundary(domain)
4848        domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
4849
4850        # Initial condition
4851        domain.set_quantity('stage', expression='elevation+0.05')
4852        domain.check_integrity()
4853
4854        assert num.allclose(domain.quantities['stage'].centroid_values,
4855                            [ 0.01296296,  0.03148148,  0.01296296,
4856                              0.03148148,  0.01296296,  0.03148148,
4857                              0.01296296,  0.03148148,  0.01296296,
4858                              0.03148148,  0.01296296,  0.03148148,
4859                             -0.04259259, -0.02407407, -0.04259259,
4860                             -0.02407407, -0.04259259, -0.02407407,
4861                             -0.04259259, -0.02407407, -0.04259259,
4862                             -0.02407407, -0.04259259, -0.02407407,
4863                             -0.09814815, -0.07962963, -0.09814815,
4864                             -0.07962963, -0.09814815, -0.07962963,
4865                             -0.09814815, -0.07962963, -0.09814815,
4866                             -0.07962963, -0.09814815, -0.07962963,
4867                             -0.1537037 , -0.13518519, -0.1537037,
4868                             -0.13518519, -0.1537037,  -0.13518519,
4869                             -0.1537037 , -0.13518519, -0.1537037,
4870                             -0.13518519, -0.1537037,  -0.13518519,
4871                             -0.20925926, -0.19074074, -0.20925926,
4872                             -0.19074074, -0.20925926, -0.19074074,
4873                             -0.20925926, -0.19074074, -0.20925926,
4874                             -0.19074074, -0.20925926, -0.19074074,
4875                             -0.26481481, -0.2462963,  -0.26481481,
4876                             -0.2462963,  -0.26481481, -0.2462963,
4877                             -0.26481481, -0.2462963,  -0.26481481,
4878                             -0.2462963,  -0.26481481, -0.2462963])
4879
4880        # Evolution
4881        for t in domain.evolve(yieldstep=0.05, finaltime=0.5):
4882            # Check that diagnostics works
4883            msg = domain.timestepping_statistics(track_speeds=True)
4884            #FIXME(Ole): One might check the contents of msg here.
4885
4886        assert num.allclose(domain.quantities['stage'].centroid_values,
4887                            [-0.02907028, -0.01475478, -0.02973417,
4888                             -0.01447186, -0.02932665, -0.01428285,
4889                             -0.02901975, -0.0141361,  -0.02898816,
4890                             -0.01418135, -0.02961409, -0.01403487,
4891                             -0.07597998, -0.06252591, -0.07664854,
4892                             -0.06312532, -0.07638287, -0.06265139,
4893                             -0.07571145, -0.06235231, -0.0756817,
4894                             -0.06245309, -0.07652292, -0.06289946,
4895                             -0.12367464, -0.11088981, -0.12237277,
4896                             -0.11115338, -0.1218934,  -0.1107174,
4897                             -0.12081485, -0.11000491, -0.12038451,
4898                             -0.11010335, -0.12102113, -0.11012105,
4899                             -0.16909116, -0.15831543, -0.16730214,
4900                             -0.15786249, -0.1665493,  -0.15697919,
4901                             -0.16496618, -0.15559852, -0.16338679,
4902                             -0.15509088, -0.16364092, -0.15424423,
4903                             -0.18771107, -0.19903904, -0.18903759,
4904                             -0.19858437, -0.18701552, -0.19697797,
4905                             -0.1833593,  -0.19505871, -0.1818806,
4906                             -0.19418042, -0.18586159, -0.19576946,
4907                             -0.13986873, -0.14170053, -0.14132188,
4908                             -0.14560674, -0.14095617, -0.14373292,
4909                             -0.13785933, -0.14033364, -0.13592955,
4910                             -0.13936356, -0.13596008, -0.14216296])
4911
4912        assert num.allclose(domain.quantities['xmomentum'].centroid_values,
4913                            [ 0.00831121,  0.00317948,  0.00731797,
4914                              0.00334939,  0.00764717,  0.00348053,
4915                              0.00788729,  0.00356522,  0.00780649,
4916                              0.00341919,  0.00693525,  0.00310375,
4917                              0.02166196,  0.01421475,  0.02017737,
4918                              0.01316839,  0.02037015,  0.01368659,
4919                              0.02106,     0.01399161,  0.02074514,
4920                              0.01354935,  0.01887407,  0.0123113,
4921                              0.03775083,  0.02855197,  0.03689337,
4922                              0.02759782,  0.03732848,  0.02812072,
4923                              0.03872545,  0.02913348,  0.03880939,
4924                              0.02803804,  0.03546499,  0.0260039,
4925                              0.0632131,   0.04730634,  0.0576324,
4926                              0.04592336,  0.05790921,  0.04690514,
4927                              0.05986467,  0.04871165,  0.06170068,
4928                              0.04811572,  0.05657041,  0.04416292,
4929                              0.08489642,  0.07188097,  0.07835261,
4930                              0.06843406,  0.07986412,  0.0698247,
4931                              0.08201071,  0.07216756,  0.08378418,
4932                              0.07273624,  0.080399,    0.06645841,
4933                              0.01631548,  0.04691608,  0.0206632,
4934                              0.044409,    0.02115518,  0.04560305,
4935                              0.02160608,  0.04663725,  0.02174734,
4936                              0.04795559,  0.02281427,  0.05667111])
4937
4938
4939        assert num.allclose(domain.quantities['ymomentum'].centroid_values,
4940                            [ 1.45876601e-004, -3.24627393e-004,
4941                             -1.57572719e-004, -2.92790187e-004,
4942                             -9.90988382e-005, -3.06677335e-004,
4943                             -1.62493106e-004, -3.71310004e-004,
4944                             -1.99445058e-004, -3.28493467e-004,
4945                              6.68217349e-005, -8.42042805e-006,
4946                              5.05093371e-004, -1.42842214e-004,
4947                             -6.81454718e-005, -5.02084057e-004,
4948                             -8.50583861e-005, -4.65443981e-004,
4949                             -1.96406564e-004, -5.88889562e-004,
4950                             -2.70160173e-004, -5.35485454e-004,
4951                              2.60780997e-004,  3.12145471e-005,
4952                              5.16189608e-004,  1.07069062e-004,
4953                              9.29989252e-005, -3.71211119e-004,
4954                              1.16350246e-004, -3.82407830e-004,
4955                             -1.62077969e-004, -6.30906636e-004,
4956                             -4.74025708e-004, -6.94463009e-004,
4957                              6.15092843e-005,  2.22106820e-004,
4958                             -6.29589294e-004,  2.43611937e-004,
4959                             -5.88125094e-004, -6.94293192e-005,
4960                             -4.17914641e-004,  6.64609019e-005,
4961                             -7.68334577e-004, -3.40232101e-004,
4962                             -1.67424308e-003, -7.39485066e-004,
4963                             -1.59966988e-003,  5.68262838e-005,
4964                             -1.48470633e-003, -1.84554882e-003,
4965                             -2.27200099e-003, -1.67506848e-003,
4966                             -1.95610258e-003, -1.47638801e-003,
4967                             -1.73779477e-003, -1.85498791e-003,
4968                             -2.01357843e-003, -2.17675471e-003,
4969                             -1.65783870e-003, -1.15818681e-003,
4970                             -1.18663036e-003, -2.94229849e-003,
4971                             -3.59309018e-003, -5.13496584e-003,
4972                             -6.17359400e-003, -5.98761937e-003,
4973                             -6.00540116e-003, -5.01121966e-003,
4974                             -4.50964850e-003, -3.06319963e-003,
4975                              6.08950810e-004, -4.79537921e-004])
4976
4977        os.remove(domain.get_name() + '.sww')
4978
4979    def NOtest_bedslope_problem_second_order_more_steps_feb_2007(self):
4980        """test_bedslope_problem_second_order_more_steps_feb_2007
4981
4982        Test shallow water finite volumes, using parameters from
4983        feb 2007 rather than backward compatibility ad infinitum
4984        """
4985
4986        from mesh_factory import rectangular
4987
4988        # Create basic mesh
4989        points, vertices, boundary = rectangular(6, 6)
4990
4991        # Create shallow water domain
4992        domain = Domain(points, vertices, boundary)
4993        domain.smooth = False
4994        domain.default_order = 2
4995        domain.beta_w = 0.9
4996        domain.beta_w_dry = 0.9
4997        domain.beta_uh = 0.9
4998        domain.beta_uh_dry = 0.9
4999        domain.beta_vh = 0.9
5000        domain.beta_vh_dry = 0.9
5001        domain.H0 = 0.001
5002        domain.tight_slope_limiters = 1
5003
5004        # Bed-slope and friction at vertices (and interpolated elsewhere)
5005        def x_slope(x, y):
5006            return -x/3
5007
5008        domain.set_quantity('elevation', x_slope)
5009
5010        # Boundary conditions
5011        Br = Reflective_boundary(domain)
5012        domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
5013
5014        # Initial condition
5015        domain.set_quantity('stage', expression='elevation+0.05')
5016        domain.check_integrity()
5017
5018        assert num.allclose(domain.quantities['stage'].centroid_values,
5019                            [ 0.01296296,  0.03148148,  0.01296296,
5020                              0.03148148,  0.01296296,  0.03148148,
5021                              0.01296296,  0.03148148,  0.01296296,
5022                              0.03148148,  0.01296296,  0.03148148,
5023                             -0.04259259, -0.02407407, -0.04259259,
5024                             -0.02407407, -0.04259259, -0.02407407,
5025                             -0.04259259, -0.02407407, -0.04259259,
5026                             -0.02407407, -0.04259259, -0.02407407,
5027                             -0.09814815, -0.07962963, -0.09814815,
5028                             -0.07962963, -0.09814815, -0.07962963,
5029                             -0.09814815, -0.07962963, -0.09814815,
5030                             -0.07962963, -0.09814815, -0.07962963,
5031                             -0.1537037 , -0.13518519, -0.1537037,
5032                             -0.13518519, -0.1537037,  -0.13518519,
5033                             -0.1537037 , -0.13518519, -0.1537037,
5034                             -0.13518519, -0.1537037,  -0.13518519,
5035                             -0.20925926, -0.19074074, -0.20925926,
5036                             -0.19074074, -0.20925926, -0.19074074,
5037                             -0.20925926, -0.19074074, -0.20925926,
5038                             -0.19074074, -0.20925926, -0.19074074,
5039                             -0.26481481, -0.2462963,  -0.26481481,
5040                             -0.2462963,  -0.26481481, -0.2462963,
5041                             -0.26481481, -0.2462963,  -0.26481481,
5042                             -0.2462963,  -0.26481481, -0.2462963])
5043
5044        # Evolution
5045        for t in domain.evolve(yieldstep = 0.05, finaltime = 0.5):
5046            pass
5047
5048        assert num.allclose(domain.quantities['stage'].centroid_values,
5049                            [-0.03348416, -0.01749303, -0.03299091,
5050                             -0.01739241, -0.03246447, -0.01732016,
5051                             -0.03205390, -0.01717833, -0.03146383,
5052                             -0.01699831, -0.03076577, -0.01671795,
5053                             -0.07952656, -0.06684763, -0.07721455,
5054                             -0.06668388, -0.07632976, -0.06600113,
5055                             -0.07523678, -0.06546373, -0.07447040,
5056                             -0.06508861, -0.07438723, -0.06359288,
5057                             -0.12526729, -0.11205668, -0.12179433,
5058                             -0.11068104, -0.12048395, -0.10968948,
5059                             -0.11912023, -0.10862628, -0.11784090,
5060                             -0.10803744, -0.11790629, -0.10742354,
5061                             -0.16859613, -0.15427413, -0.16664444,
5062                             -0.15464452, -0.16570816, -0.15327556,
5063                             -0.16409162, -0.15204092, -0.16264608,
5064                             -0.15102139, -0.16162736, -0.14969205,
5065                             -0.18736511, -0.19874036, -0.18811230,
5066                             -0.19758289, -0.18590182, -0.19580301,
5067                             -0.18234588, -0.19423215, -0.18100376,
5068                             -0.19380116, -0.18509710, -0.19501636,
5069                             -0.13982382, -0.14166819, -0.14132775,
5070                             -0.14528694, -0.14096905, -0.14351126,
5071                             -0.13800356, -0.14027920, -0.13613538,
5072                             -0.13936795, -0.13621902, -0.14204982]) 
5073
5074        assert num.allclose(domain.quantities['xmomentum'].centroid_values,
5075                            [0.00600290,  0.00175780,  0.00591905,
5076                             0.00190903,  0.00644462,  0.00203095,
5077                             0.00684561,  0.00225089,  0.00708208,
5078                             0.00236235,  0.00649095,  0.00222343,
5079                             0.02068693,  0.01164034,  0.01983343,
5080                             0.01159526,  0.02044611,  0.01233252,
5081                             0.02135685,  0.01301289,  0.02161290,
5082                             0.01260280,  0.01867612,  0.01133078,
5083                             0.04091313,  0.02668283,  0.03634781,
5084                             0.02733469,  0.03767692,  0.02836840,
5085                             0.03906338,  0.02958073,  0.04025669,
5086                             0.02953292,  0.03665616,  0.02583565,
5087                             0.06314558,  0.04830935,  0.05663609,
5088                             0.04564362,  0.05756200,  0.04739673,
5089                             0.05967379,  0.04919083,  0.06124330,
5090                             0.04965808,  0.05879240,  0.04629319,
5091                             0.08220739,  0.06924725,  0.07713556,
5092                             0.06782640,  0.07909499,  0.06992544,
5093                             0.08116621,  0.07210181,  0.08281548,
5094                             0.07222669,  0.07941059,  0.06755612,
5095                             0.01581588,  0.04533609,  0.02017939,
5096                             0.04342565,  0.02073232,  0.04476108,
5097                             0.02117439,  0.04573358,  0.02129473,
5098                             0.04694267,  0.02220398,  0.05533458])
5099
5100
5101        assert num.allclose(domain.quantities['ymomentum'].centroid_values,
5102                            [-7.65882069e-005, -1.46087080e-004,
5103                             -1.09630102e-004, -7.80950424e-005,
5104                             -1.15922807e-005, -9.09134899e-005,
5105                             -1.35994542e-004, -1.95673476e-004,
5106                             -4.25779199e-004, -2.95890312e-004,
5107                             -4.00060341e-004, -9.42021290e-005,
5108                             -3.41372596e-004, -1.54560195e-004,
5109                             -2.94810038e-004, -1.08844546e-004,
5110                             -6.97240892e-005,  3.50299623e-005,
5111                             -2.40159184e-004, -2.01805883e-004,
5112                             -7.60732405e-004, -5.10897642e-004,
5113                             -1.00940001e-003, -1.38037759e-004,
5114                             -1.06169131e-003, -3.12307760e-004,
5115                             -9.90602307e-004, -4.21634250e-005,
5116                             -6.02424239e-004,  1.52230578e-004,
5117                             -7.63833035e-004, -1.10273481e-004,
5118                             -1.40187071e-003, -5.57831837e-004,
5119                             -1.63988285e-003, -2.48018092e-004,
5120                             -1.83309840e-003, -6.19360836e-004,
5121                             -1.29955242e-003, -3.76237145e-004,
5122                             -1.00613007e-003, -8.63641918e-005,
5123                             -1.13604124e-003, -3.90589728e-004,
5124                             -1.91457355e-003, -9.43783961e-004,
5125                             -2.28090840e-003, -5.79107025e-004,
5126                             -1.54091533e-003, -2.39785792e-003,
5127                             -2.47947427e-003, -2.02694009e-003,
5128                             -2.10441194e-003, -1.82082650e-003,
5129                             -1.80229336e-003, -2.10418336e-003,
5130                             -1.93104408e-003, -2.23200334e-003,
5131                             -1.57239706e-003, -1.31486358e-003,
5132                             -1.17564993e-003, -2.85846494e-003,
5133                             -3.52956754e-003, -5.12658193e-003,
5134                             -6.24238960e-003, -6.01820113e-003,
5135                             -6.09602201e-003, -5.04787190e-003,
5136                             -4.59373845e-003, -3.01393146e-003,
5137                             5.08550095e-004, -4.35896549e-004])
5138
5139        os.remove(domain.get_name() + '.sww')
5140
5141    def test_temp_play(self):
5142        from mesh_factory import rectangular
5143
5144        # Create basic mesh
5145        points, vertices, boundary = rectangular(5, 5)
5146
5147        # Create shallow water domain
5148        domain = Domain(points, vertices, boundary)
5149        domain.smooth = False
5150        domain.default_order = 2
5151        domain.beta_w = 0.9
5152        domain.beta_w_dry = 0.9
5153        domain.beta_uh = 0.9
5154        domain.beta_uh_dry = 0.9
5155        domain.beta_vh = 0.9
5156        domain.beta_vh_dry = 0.9
5157
5158        # FIXME (Ole): Need tests where these two are commented out
5159        domain.H0 = 0                         # Backwards compatibility (6/2/7)
5160        domain.tight_slope_limiters = False   # Backwards compatibility (14/4/7)
5161        domain.use_centroid_velocities = False # Backwards compatibility (7/5/8)
5162        domain.use_edge_limiter = False       # Backwards compatibility (9/5/8)
5163
5164
5165        # Bed-slope and friction at vertices (and interpolated elsewhere)
5166        def x_slope(x, y):
5167            return -x/3
5168
5169        domain.set_quantity('elevation', x_slope)
5170
5171        # Boundary conditions
5172        Br = Reflective_boundary(domain)
5173        domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
5174
5175        # Initial condition
5176        domain.set_quantity('stage', expression='elevation+0.05')
5177        domain.check_integrity()
5178
5179        # Evolution
5180        for t in domain.evolve(yieldstep=0.05, finaltime=0.1):
5181            pass
5182
5183        assert num.allclose(domain.quantities['stage'].centroid_values[:4],
5184                            [0.00206836, 0.01296714, 0.00363415, 0.01438924])
5185        assert num.allclose(domain.quantities['xmomentum'].centroid_values[:4],
5186                            [0.01360154, 0.00671133, 0.01264578, 0.00648503])
5187        assert num.allclose(domain.quantities['ymomentum'].centroid_values[:4],
5188                            [-1.19201077e-003, -7.23647546e-004,
5189                             -6.39083123e-005, 6.29815168e-005])
5190
5191        os.remove(domain.get_name() + '.sww')
5192
5193    def test_complex_bed(self):
5194        # No friction is tested here
5195
5196        from mesh_factory import rectangular
5197
5198        N = 12
5199        points, vertices, boundary = rectangular(N, N/2, len1=1.2, len2=0.6,
5200                                                 origin=(-0.07, 0))
5201
5202
5203        domain = Domain(points, vertices, boundary)
5204        domain.smooth = False
5205        domain.default_order = 2
5206
5207        inflow_stage = 0.1
5208        Z = Weir(inflow_stage)
5209        domain.set_quantity('elevation', Z)
5210
5211        Br = Reflective_boundary(domain)
5212        Bd = Dirichlet_boundary([inflow_stage, 0.0, 0.0])
5213        domain.set_boundary({'left': Bd, 'right': Br, 'bottom': Br, 'top': Br})
5214
5215        domain.set_quantity('stage', expression='elevation')
5216
5217        for t in domain.evolve(yieldstep=0.02, finaltime=0.2):
5218            pass
5219
5220        #FIXME: These numbers were from version before 25/10
5221        #assert allclose(domain.quantities['stage'].centroid_values,
5222# [3.95822638e-002,  5.61022588e-002,  4.66437868e-002,  5.73081011e-002,
5223#  4.72394613e-002,  5.74684939e-002,  4.74309483e-002,  5.77458084e-002,
5224#  4.80628177e-002,  5.85656225e-002,  4.90498542e-002,  6.02609831e-002,
5225#  1.18470315e-002,  1.75136443e-002,  1.18035266e-002,  2.15565695e-002,
5226#  1.31620268e-002,  2.14351640e-002,  1.32351076e-002,  2.15450687e-002,
5227#  1.36414028e-002,  2.24274619e-002,  1.51689511e-002,  2.21789655e-002,
5228# -7.54337535e-003, -6.86362021e-004, -7.74146760e-003, -1.83756530e-003,
5229# -8.16773628e-003, -4.49916813e-004, -8.08202599e-003, -3.91118720e-004,
5230# -8.10292716e-003, -3.88584984e-004, -7.35226124e-003,  2.73985295e-004,
5231#  1.86166683e-001,  8.74070369e-002,  1.86166712e-001,  8.74035875e-002,
5232#  6.11666935e-002, -3.76173225e-002, -6.38333276e-002, -3.76147365e-002,
5233#  6.11666725e-002,  8.73846774e-002,  1.86166697e-001,  8.74171550e-002,
5234# -4.83333333e-002,  1.18333333e-001, -4.83333333e-002,  1.18333333e-001,
5235# -4.83333333e-002, -6.66666667e-003, -1.73333333e-001, -1.31666667e-001,
5236# -1.73333333e-001, -6.66666667e-003, -4.83333333e-002,  1.18333333e-001,
5237# -2.48333333e-001, -2.31666667e-001, -2.48333333e-001, -2.31666667e-001,
5238# -2.48333333e-001, -2.31666667e-001, -2.48333333e-001, -2.31666667e-001,
5239# -2.48333333e-001, -2.31666667e-001, -2.48333333e-001, -2.31666667e-001,
5240# -4.65000000e-001, -3.65000000e-001, -4.65000000e-001, -3.65000000e-001,
5241# -4.65000000e-001, -3.65000000e-001, -4.65000000e-001, -3.65000000e-001,
5242# -4.65000000e-001, -3.65000000e-001, -4.65000000e-001, -3.65000000e-001,
5243# -5.98333333e-001, -5.81666667e-001, -5.98333333e-001, -5.81666667e-001,
5244# -5.98333333e-001, -5.81666667e-001, -5.98333333e-001, -5.81666667e-001,
5245# -5.98333333e-001, -5.81666667e-001, -5.98333333e-001, -5.81666667e-001,
5246# -6.48333333e-001, -6.31666667e-001, -6.48333333e-001, -6.31666667e-001,
5247# -6.48333333e-001, -6.31666667e-001, -6.48333333e-001, -6.31666667e-001,
5248# -6.48333333e-001, -6.31666667e-001, -6.48333333e-001, -6.31666667e-001,
5249# -5.31666667e-001, -5.98333333e-001, -5.31666667e-001, -5.98333333e-001,
5250# -5.31666667e-001, -5.98333333e-001, -5.31666667e-001, -5.98333333e-001,
5251# -5.31666667e-001, -5.98333333e-001, -5.31666667e-001, -5.98333333e-001,
5252# -4.98333333e-001, -4.81666667e-001, -4.98333333e-001, -4.81666667e-001,
5253# -4.98333333e-001, -4.81666667e-001, -4.98333333e-001, -4.81666667e-001,
5254# -4.98333333e-001, -4.81666667e-001, -4.98333333e-001, -4.81666667e-001,
5255# -5.48333333e-001, -5.31666667e-001, -5.48333333e-001, -5.31666667e-001,
5256# -5.48333333e-001, -5.31666667e-001, -5.48333333e-001, -5.31666667e-001,
5257# -5.48333333e-001, -5.31666667e-001, -5.48333333e-001, -5.31666667e-001])
5258
5259        os.remove(domain.get_name() + '.sww')
5260
5261    def test_spatio_temporal_boundary_1(self):
5262        """Test that boundary values can be read from file and interpolated
5263        in both time and space.
5264
5265        Verify that the same steady state solution is arrived at and that
5266        time interpolation works.
5267
5268        The full solution history is not exactly the same as
5269        file boundary must read and interpolate from *smoothed* version
5270        as stored in sww.
5271        """
5272
5273        # Create sww file of simple propagation from left to right
5274        # through rectangular domain
5275
5276        import time
5277        from mesh_factory import rectangular
5278
5279        # Create basic mesh
5280        points, vertices, boundary = rectangular(3, 3)
5281
5282        # Create shallow water domain
5283        domain1 = Domain(points, vertices, boundary)
5284
5285        domain1.reduction = mean
5286        domain1.smooth = False    # Exact result
5287
5288        domain1.default_order = 2
5289        domain1.store = True
5290        domain1.set_datadir('.')
5291        domain1.set_name('spatio_temporal_boundary_source' + str(time.time()))
5292
5293        # FIXME: This is extremely important!
5294        # How can we test if they weren't stored?
5295        domain1.quantities_to_be_stored = ['stage', 'xmomentum', 'ymomentum']
5296
5297
5298        # Bed-slope and friction at vertices (and interpolated elsewhere)
5299        domain1.set_quantity('elevation', 0)
5300        domain1.set_quantity('friction', 0)
5301
5302        # Boundary conditions
5303        Br = Reflective_boundary(domain1)
5304        Bd = Dirichlet_boundary([0.3,0,0])
5305        domain1.set_boundary({'left': Bd, 'top': Bd, 'right': Br, 'bottom': Br})
5306
5307        # Initial condition
5308        domain1.set_quantity('stage', 0)
5309        domain1.check_integrity()
5310
5311        finaltime = 5
5312
5313        # Evolution  (full domain - large steps)
5314        for t in domain1.evolve(yieldstep=0.671, finaltime=finaltime):
5315            pass
5316
5317        cv1 = domain1.quantities['stage'].centroid_values
5318
5319        # Create a triangle shaped domain (reusing coordinates from domain 1),
5320        # formed from the lower and right hand  boundaries and
5321        # the sw-ne diagonal
5322        # from domain 1. Call it domain2
5323
5324        points = [[0,0], [1.0/3,0], [1.0/3,1.0/3],
5325                  [2.0/3,0], [2.0/3,1.0/3], [2.0/3,2.0/3],
5326                  [1,0], [1,1.0/3], [1,2.0/3], [1,1]]
5327
5328        vertices = [[1,2,0], [3,4,1], [2,1,4], [4,5,2],
5329                    [6,7,3], [4,3,7], [7,8,4], [5,4,8], [8,9,5]]
5330
5331        boundary = {(0,1): 'bottom',   (1,1): 'bottom',   (4,1): 'bottom',
5332                    (4,2): 'right',    (6,2): 'right',    (8,2): 'right',
5333                    (0,0): 'diagonal', (3,0): 'diagonal', (8,0): 'diagonal'}
5334
5335        domain2 = Domain(points, vertices, boundary)
5336
5337        domain2.reduction = domain1.reduction
5338        domain2.smooth = False
5339        domain2.default_order = 2
5340
5341        # Bed-slope and friction at vertices (and interpolated elsewhere)
5342        domain2.set_quantity('elevation', 0)
5343        domain2.set_quantity('friction', 0)
5344        domain2.set_quantity('stage', 0)
5345
5346        # Boundary conditions
5347        Br = Reflective_boundary(domain2)
5348        Bf = Field_boundary(domain1.get_name() + '.' + domain1.format, domain2)
5349        domain2.set_boundary({'right':Br, 'bottom':Br, 'diagonal':Bf})
5350        domain2.check_integrity()
5351
5352        # Evolution (small steps)
5353        for t in domain2.evolve(yieldstep=0.0711, finaltime=finaltime):
5354            pass
5355
5356        # Use output from domain1 as spatio-temporal boundary for domain2
5357        # and verify that results at right hand side are close.
5358        cv2 = domain2.quantities['stage'].centroid_values
5359
5360        assert num.allclose(num.take(cv1, (0,8,16), axis=0),
5361                            num.take(cv2, (0,3,8), axis=0))      # Diag
5362        assert num.allclose(num.take(cv1, (0,6,12), axis=0),
5363                            num.take(cv2, (0,1,4), axis=0))      # Bottom
5364        assert num.allclose(num.take(cv1, (12,14,16), axis=0),
5365                            num.take(cv2, (4,6,8), axis=0))      # RHS
5366
5367        # Cleanup
5368        os.remove(domain1.get_name() + '.' + domain1.format)
5369        os.remove(domain2.get_name() + '.' + domain2.format)
5370
5371    def test_spatio_temporal_boundary_2(self):
5372        """Test that boundary values can be read from file and interpolated
5373        in both time and space.
5374        This is a more basic test, verifying that boundary object
5375        produces the expected results
5376        """
5377
5378        import time
5379        from mesh_factory import rectangular
5380
5381        # Create sww file of simple propagation from left to right
5382        # through rectangular domain
5383
5384        # Create basic mesh
5385        points, vertices, boundary = rectangular(3, 3)
5386
5387        #Create shallow water domain
5388        domain1 = Domain(points, vertices, boundary)
5389
5390        domain1.reduction = mean
5391        domain1.smooth = True    # To mimic MOST output
5392
5393        domain1.default_order = 2
5394        domain1.store = True
5395        domain1.set_datadir('.')
5396        domain1.set_name('spatio_temporal_boundary_source' + str(time.time()))
5397
5398        # FIXME: This is extremely important!
5399        # How can we test if they weren't stored?
5400        domain1.quantities_to_be_stored = ['stage', 'xmomentum', 'ymomentum']
5401
5402        # Bed-slope and friction at vertices (and interpolated elsewhere)
5403        domain1.set_quantity('elevation', 0)
5404        domain1.set_quantity('friction', 0)
5405
5406        # Boundary conditions
5407        Br = Reflective_boundary(domain1)
5408        Bd = Dirichlet_boundary([0.3,0,0])
5409        domain1.set_boundary({'left': Bd, 'top': Bd, 'right': Br, 'bottom': Br})
5410
5411        # Initial condition
5412        domain1.set_quantity('stage', 0)
5413        domain1.check_integrity()
5414
5415        finaltime = 5
5416
5417        # Evolution  (full domain - large steps)
5418        for t in domain1.evolve(yieldstep=1, finaltime=finaltime):
5419            pass
5420
5421        # Create an triangle shaped domain (coinciding with some
5422        # coordinates from domain 1),
5423        # formed from the lower and right hand  boundaries and
5424        # the sw-ne diagonal
5425        # from domain 1. Call it domain2
5426        points = [[0,0],
5427                  [1.0/3,0], [1.0/3,1.0/3],
5428                  [2.0/3,0], [2.0/3,1.0/3], [2.0/3,2.0/3],
5429                  [1,0],     [1,1.0/3],     [1,2.0/3],     [1,1]]
5430
5431        vertices = [[1,2,0], [3,4,1], [2,1,4], [4,5,2],
5432                    [6,7,3], [4,3,7], [7,8,4], [5,4,8], [8,9,5]]
5433
5434        boundary = {(0,1): 'bottom',   (1,1): 'bottom',   (4,1): 'bottom',
5435                    (4,2): 'right',    (6,2): 'right',    (8,2): 'right',
5436                    (0,0): 'diagonal', (3,0): 'diagonal', (8,0): 'diagonal'}
5437
5438        domain2 = Domain(points, vertices, boundary)
5439
5440        domain2.reduction = domain1.reduction
5441        domain2.smooth = False
5442        domain2.default_order = 2
5443
5444        # Bed-slope and friction at vertices (and interpolated elsewhere)
5445        domain2.set_quantity('elevation', 0)
5446        domain2.set_quantity('friction', 0)
5447        domain2.set_quantity('stage', 0)
5448
5449        # Read results for specific timesteps t=1 and t=2
5450        from Scientific.IO.NetCDF import NetCDFFile
5451        fid = NetCDFFile(domain1.get_name() + '.' + domain1.format)
5452
5453        x = fid.variables['x'][:]
5454        y = fid.variables['y'][:]
5455        s1 = fid.variables['stage'][1,:]
5456        s2 = fid.variables['stage'][2,:]
5457        fid.close()
5458
5459        shp = (len(x), 1)
5460        points = num.concatenate((num.reshape(x, shp), num.reshape(y, shp)),
5461                                 axis=1)
5462
5463        # The diagonal points of domain 1 are 0, 5, 10, 15
5464        msg = ('value was\n%s\nshould be\n'
5465               '[[0,0], [1.0/3, 1.0/3],\n'
5466               '[2.0/3, 2.0/3], [1,1]]'
5467               % str(num.take(points, [0,5,10,15], axis=0)))
5468        assert num.allclose(num.take(points, [0,5,10,15], axis=0),
5469                            [[0,0], [1.0/3, 1.0/3], [2.0/3, 2.0/3], [1,1]]), msg
5470
5471        # Boundary conditions
5472        Br = Reflective_boundary(domain2)
5473        Bf = Field_boundary(domain1.get_name() + '.' + domain1.format,
5474                            domain2, verbose=False)
5475        domain2.set_boundary({'right':Br, 'bottom':Br, 'diagonal':Bf})
5476        domain2.check_integrity()
5477
5478        # Test that interpolation points are the mid points of the all boundary
5479        # segments
5480        boundary_midpoints = [[1.0/6, 0], [1.0/2, 0], [5.0/6,0],
5481                              [1.0, 1.0/6], [1.0, 1.0/2], [1.0, 5.0/6],
5482                              [1.0/6, 1.0/6], [0.5, 0.5], [5.0/6, 5.0/6]]
5483
5484        boundary_midpoints.sort()
5485        R = Bf.F.interpolation_points.tolist()
5486        R.sort()
5487        assert num.allclose(boundary_midpoints, R)
5488
5489        # Check spatially interpolated output at time == 1
5490        domain2.time = 1
5491
5492        # First diagonal midpoint
5493        R0 = Bf.evaluate(0, 0)
5494        assert num.allclose(R0[0], (s1[0] + s1[5])/2)
5495
5496        # Second diagonal midpoint
5497        R0 = Bf.evaluate(3, 0)
5498        assert num.allclose(R0[0], (s1[5] + s1[10])/2)
5499
5500        # First diagonal midpoint
5501        R0 = Bf.evaluate(8, 0)
5502        assert num.allclose(R0[0], (s1[10] + s1[15])/2)
5503
5504        # Check spatially interpolated output at time == 2
5505        domain2.time = 2
5506
5507        # First diagonal midpoint
5508        R0 = Bf.evaluate(0, 0)
5509        assert num.allclose(R0[0], (s2[0] + s2[5])/2)
5510
5511        # Second diagonal midpoint
5512        R0 = Bf.evaluate(3, 0)
5513        assert num.allclose(R0[0], (s2[5] + s2[10])/2)
5514
5515        # First diagonal midpoint
5516        R0 = Bf.evaluate(8, 0)
5517        assert num.allclose(R0[0], (s2[10] + s2[15])/2)
5518
5519        # Now check temporal interpolation
5520        domain2.time = 1 + 2.0/3
5521
5522        # First diagonal midpoint
5523        R0 = Bf.evaluate(0,0)
5524        assert num.allclose(R0[0],
5525                            ((s1[0] + s1[5])/2 + 2.0*(s2[0] + s2[5])/2)/3)
5526
5527        # Second diagonal midpoint
5528        R0 = Bf.evaluate(3, 0)
5529        assert num.allclose(R0[0],
5530                            ((s1[5] + s1[10])/2 + 2.0*(s2[5] + s2[10])/2)/3)
5531
5532        # First diagonal midpoint
5533        R0 = Bf.evaluate(8, 0)
5534        assert num.allclose(R0[0],
5535                            ((s1[10] + s1[15])/2 + 2.0*(s2[10] + s2[15])/2)/3)
5536
5537        # Cleanup
5538        os.remove(domain1.get_name() + '.' + domain1.format)
5539
5540    def test_spatio_temporal_boundary_3(self):
5541        """Test that boundary values can be read from file and interpolated
5542        in both time and space.
5543        This is a more basic test, verifying that boundary object
5544        produces the expected results
5545
5546        This tests adjusting using mean_stage
5547        """
5548
5549        #Create sww file of simple propagation from left to right
5550        #through rectangular domain
5551
5552        import time
5553        from mesh_factory import rectangular
5554
5555        mean_stage = 5.2    # Adjust stage by this amount in boundary
5556
5557        # Create basic mesh
5558        points, vertices, boundary = rectangular(3, 3)
5559
5560        # Create shallow water domain
5561        domain1 = Domain(points, vertices, boundary)
5562
5563        domain1.reduction = mean
5564        domain1.smooth = True    # To mimic MOST output
5565
5566        domain1.default_order = 2
5567        domain1.store = True
5568        domain1.set_datadir('.')
5569        domain1.set_name('spatio_temporal_boundary_source' + str(time.time()))
5570
5571        # FIXME: This is extremely important!
5572        # How can we test if they weren't stored?
5573        domain1.quantities_to_be_stored = ['stage', 'xmomentum', 'ymomentum']
5574
5575        # Bed-slope and friction at vertices (and interpolated elsewhere)
5576        domain1.set_quantity('elevation', 0)
5577        domain1.set_quantity('friction', 0)
5578
5579        # Boundary conditions
5580        Br = Reflective_boundary(domain1)
5581        Bd = Dirichlet_boundary([0.3, 0, 0])
5582        domain1.set_boundary({'left': Bd, 'top': Bd, 'right': Br, 'bottom': Br})
5583
5584        # Initial condition
5585        domain1.set_quantity('stage', 0)
5586        domain1.check_integrity()
5587
5588        finaltime = 5
5589
5590        # Evolution  (full domain - large steps)
5591        for t in domain1.evolve(yieldstep=1, finaltime=finaltime):
5592            pass
5593
5594        # Create an triangle shaped domain (coinciding with some
5595        # coordinates from domain 1),
5596        # formed from the lower and right hand  boundaries and
5597        # the sw-ne diagonal
5598        # from domain 1. Call it domain2
5599        points = [[0,0],
5600                  [1.0/3,0], [1.0/3,1.0/3],
5601                  [2.0/3,0], [2.0/3,1.0/3], [2.0/3,2.0/3],
5602                  [1,0],     [1,1.0/3],     [1,2.0/3],     [1,1]]
5603
5604        vertices = [[1,2,0],
5605                    [3,4,1], [2,1,4], [4,5,2],
5606                    [6,7,3], [4,3,7], [7,8,4], [5,4,8], [8,9,5]]
5607
5608        boundary = {(0,1): 'bottom',   (1,1): 'bottom',   (4,1): 'bottom',
5609                    (4,2): 'right',    (6,2): 'right',    (8,2): 'right',
5610                    (0,0): 'diagonal', (3,0): 'diagonal', (8,0): 'diagonal'}
5611
5612        domain2 = Domain(points, vertices, boundary)
5613
5614        domain2.reduction = domain1.reduction
5615        domain2.smooth = False
5616        domain2.default_order = 2
5617
5618        # Bed-slope and friction at vertices (and interpolated elsewhere)
5619        domain2.set_quantity('elevation', 0)
5620        domain2.set_quantity('friction', 0)
5621        domain2.set_quantity('stage', 0)
5622
5623        # Read results for specific timesteps t=1 and t=2
5624        from Scientific.IO.NetCDF import NetCDFFile
5625        fid = NetCDFFile(domain1.get_name() + '.' + domain1.format)
5626
5627        x = fid.variables['x'][:]
5628        y = fid.variables['y'][:]
5629        s1 = fid.variables['stage'][1,:]
5630        s2 = fid.variables['stage'][2,:]
5631        fid.close()
5632
5633        shp = (len(x), 1)
5634        points = num.concatenate((num.reshape(x, shp), num.reshape(y, shp)),
5635                                 axis=1)
5636        #The diagonal points of domain 1 are 0, 5, 10, 15
5637
5638        msg = ('values was\n%s\nshould be\n'
5639               '[[0,0], [1.0/3, 1.0/3],\n'
5640               '[2.0/3, 2.0/3], [1,1]]'
5641               % str(num.take(points, [0,5,10,15], axis=0)))
5642        assert num.allclose(num.take(points, [0,5,10,15], axis=0),
5643                            [[0,0], [1.0/3, 1.0/3], [2.0/3, 2.0/3], [1,1]]), msg
5644
5645        # Boundary conditions
5646        Br = Reflective_boundary(domain2)
5647        Bf = Field_boundary(domain1.get_name() + '.' + domain1.format,
5648                            domain2, mean_stage=mean_stage, verbose=False)
5649
5650        domain2.set_boundary({'right':Br, 'bottom':Br, 'diagonal':Bf})
5651        domain2.check_integrity()
5652
5653        # Test that interpolation points are the mid points of the all boundary
5654        # segments
5655        boundary_midpoints = [[1.0/6, 0], [1.0/2, 0], [5.0/6,0],
5656                              [1.0, 1.0/6], [1.0, 1.0/2], [1.0, 5.0/6],
5657                              [1.0/6, 1.0/6], [0.5, 0.5], [5.0/6, 5.0/6]]
5658
5659        boundary_midpoints.sort()
5660        R = Bf.F.interpolation_points.tolist()
5661        R.sort()
5662        assert num.allclose(boundary_midpoints, R)
5663
5664        # Check spatially interpolated output at time == 1
5665        domain2.time = 1
5666
5667        # First diagonal midpoint
5668        R0 = Bf.evaluate(0, 0)
5669        assert num.allclose(R0[0], (s1[0] + s1[5])/2 + mean_stage)
5670
5671        # Second diagonal midpoint
5672        R0 = Bf.evaluate(3, 0)
5673        assert num.allclose(R0[0], (s1[5] + s1[10])/2 + mean_stage)
5674
5675        # First diagonal midpoint
5676        R0 = Bf.evaluate(8, 0)
5677        assert num.allclose(R0[0], (s1[10] + s1[15])/2 + mean_stage)
5678
5679        # Check spatially interpolated output at time == 2
5680        domain2.time = 2
5681
5682        # First diagonal midpoint
5683        R0 = Bf.evaluate(0, 0)
5684        assert num.allclose(R0[0], (s2[0] + s2[5])/2 + mean_stage)
5685
5686        # Second diagonal midpoint
5687        R0 = Bf.evaluate(3, 0)
5688        assert num.allclose(R0[0], (s2[5] + s2[10])/2 + mean_stage)
5689
5690        # First diagonal midpoint
5691        R0 = Bf.evaluate(8, 0)
5692        assert num.allclose(R0[0], (s2[10] + s2[15])/2 + mean_stage)
5693
5694        #Now check temporal interpolation
5695        domain2.time = 1 + 2.0/3
5696
5697        # First diagonal midpoint
5698        R0 = Bf.evaluate(0, 0)
5699        assert num.allclose(R0[0],
5700                            ((s1[0] + s1[5])/2 + 2.0*(s2[0] + s2[5])/2)/3 +
5701                                mean_stage)
5702
5703        # Second diagonal midpoint
5704        R0 = Bf.evaluate(3, 0)
5705        assert num.allclose(R0[0],
5706                            ((s1[5] + s1[10])/2 + 2.0*(s2[5] + s2[10])/2)/3 +
5707                                mean_stage)
5708
5709        # First diagonal midpoint
5710        R0 = Bf.evaluate(8, 0)
5711        assert num.allclose(R0[0],
5712                            ((s1[10] + s1[15])/2 + 2.0*(s2[10] + s2[15])/2)/3 +
5713                                mean_stage)
5714
5715        # Cleanup
5716        os.remove(domain1.get_name() + '.' + domain1.format)
5717
5718    def test_spatio_temporal_boundary_outside(self):
5719        """Test that field_boundary catches if a point is outside the sww
5720        that defines it
5721        """
5722
5723        import time
5724        from mesh_factory import rectangular
5725
5726        # Create sww file of simple propagation from left to right
5727        # through rectangular domain
5728
5729        # Create basic mesh
5730        points, vertices, boundary = rectangular(3, 3)
5731
5732        # Create shallow water domain
5733        domain1 = Domain(points, vertices, boundary)
5734
5735        domain1.reduction = mean
5736        domain1.smooth = True    # To mimic MOST output
5737
5738        domain1.default_order = 2
5739        domain1.store = True
5740        domain1.set_datadir('.')
5741        domain1.set_name('spatio_temporal_boundary_source' + str(time.time()))
5742
5743        # FIXME: This is extremely important!
5744        # How can we test if they weren't stored?
5745        domain1.quantities_to_be_stored = ['stage', 'xmomentum', 'ymomentum']
5746
5747
5748        # Bed-slope and friction at vertices (and interpolated elsewhere)
5749        domain1.set_quantity('elevation', 0)
5750        domain1.set_quantity('friction', 0)
5751
5752        # Boundary conditions
5753        Br = Reflective_boundary(domain1)
5754        Bd = Dirichlet_boundary([0.3, 0, 0])
5755        domain1.set_boundary({'left': Bd, 'top': Bd, 'right': Br, 'bottom': Br})
5756
5757        # Initial condition
5758        domain1.set_quantity('stage', 0)
5759        domain1.check_integrity()
5760
5761        finaltime = 5
5762
5763        # Evolution  (full domain - large steps)
5764        for t in domain1.evolve(yieldstep=1, finaltime=finaltime):
5765            pass
5766
5767        # Create an triangle shaped domain (coinciding with some
5768        # coordinates from domain 1, but one edge outside!),
5769        # formed from the lower and right hand  boundaries and
5770        # the sw-ne diagonal as in the previous test but scaled
5771        # in the x direction by a factor of 2
5772        points = [[0,0],
5773                  [2.0/3,0], [2.0/3,1.0/3],
5774                  [4.0/3,0], [4.0/3,1.0/3], [4.0/3,2.0/3],
5775                  [2,0],     [2,1.0/3],     [2,2.0/3],     [2,1]]
5776
5777        vertices = [[1,2,0],
5778                    [3,4,1], [2,1,4], [4,5,2],
5779                    [6,7,3], [4,3,7], [7,8,4], [5,4,8], [8,9,5]]
5780
5781        boundary = {(0,1): 'bottom',   (1,1): 'bottom',   (4,1): 'bottom',
5782                    (4,2): 'right',    (6,2): 'right',    (8,2): 'right',
5783                    (0,0): 'diagonal', (3,0): 'diagonal', (8,0): 'diagonal'}
5784
5785        domain2 = Domain(points, vertices, boundary)
5786
5787        domain2.reduction = domain1.reduction
5788        domain2.smooth = False
5789        domain2.default_order = 2
5790
5791        # Bed-slope and friction at vertices (and interpolated elsewhere)
5792        domain2.set_quantity('elevation', 0)
5793        domain2.set_quantity('friction', 0)
5794        domain2.set_quantity('stage', 0)
5795
5796        # Read results for specific timesteps t=1 and t=2
5797        from Scientific.IO.NetCDF import NetCDFFile
5798        fid = NetCDFFile(domain1.get_name() + '.' + domain1.format)
5799
5800        x = fid.variables['x'][:]
5801        y = fid.variables['y'][:]
5802        s1 = fid.variables['stage'][1,:]
5803        s2 = fid.variables['stage'][2,:]
5804        fid.close()
5805
5806        shp = (len(x), 1)
5807        points = num.concatenate((num.reshape(x, shp), num.reshape(y, shp)),
5808                                 axis=1)
5809        #The diagonal points of domain 1 are 0, 5, 10, 15
5810
5811        assert num.allclose(num.take(points, [0,5,10,15], axis=0),
5812                            [[0,0], [1.0/3,1.0/3], [2.0/3,2.0/3], [1,1]])
5813
5814        # Boundary conditions
5815        Br = Reflective_boundary(domain2)
5816        Bf = Field_boundary(domain1.get_name() + '.' + domain1.format,
5817                            domain2, mean_stage=1, verbose=False)
5818
5819        domain2.set_boundary({'right':Br, 'bottom':Br, 'diagonal':Bf})
5820        domain2.check_integrity()
5821
5822        try:
5823            for t in domain2.evolve(yieldstep=1, finaltime=finaltime):
5824                pass
5825        except:
5826            pass
5827        else:
5828            msg = 'This should have caught NAN at boundary'
5829            raise Exception, msg
5830
5831        #Cleanup
5832        os.remove(domain1.get_name() + '.' + domain1.format)
5833
5834    def test_extrema(self):
5835        """Test that extrema of quantities are computed correctly
5836        Extrema are updated at every *internal* timestep
5837        """
5838
5839        from anuga.abstract_2d_finite_volumes.mesh_factory \
5840                import rectangular_cross
5841
5842        initial_runup_height = -0.4
5843        final_runup_height = -0.3
5844
5845        #--------------------------------------------------------------
5846        # Setup computational domain
5847        #--------------------------------------------------------------
5848        N = 5
5849        points, vertices, boundary = rectangular_cross(N, N)
5850        domain = Domain(points, vertices, boundary)
5851        domain.set_name('extrema_test')
5852
5853        #--------------------------------------------------------------
5854        # Setup initial conditions
5855        #--------------------------------------------------------------
5856        def topography(x,y):
5857            return -x/2                             # linear bed slope
5858
5859        domain.set_quantity('elevation', topography)    # function for elevation
5860        domain.set_quantity('friction', 0.)             # Zero friction
5861        # Constant negative initial stage
5862        domain.set_quantity('stage', initial_runup_height)
5863        domain.set_quantities_to_be_monitored(['stage', 'stage-elevation'],
5864                                              time_interval=[0.5, 2.7],
5865                                              polygon=[[0,0], [0,1],
5866                                                       [1,1], [1,0]])
5867
5868        assert len(domain.quantities_to_be_monitored) == 2
5869        assert domain.quantities_to_be_monitored.has_key('stage')
5870        assert domain.quantities_to_be_monitored.has_key('stage-elevation')
5871        for key in domain.quantities_to_be_monitored['stage'].keys():
5872            assert domain.quantities_to_be_monitored['stage'][key] is None
5873
5874        #--------------------------------------------------------------
5875        # Setup boundary conditions
5876        #--------------------------------------------------------------
5877        Br = Reflective_boundary(domain)              # Reflective wall
5878        # Constant inflow
5879        Bd = Dirichlet_boundary([final_runup_height, 0, 0])
5880
5881        # All reflective to begin with (still water)
5882        domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
5883
5884        #--------------------------------------------------------------
5885        # Let triangles adjust and check extrema
5886        #--------------------------------------------------------------
5887        for t in domain.evolve(yieldstep=0.1, finaltime=1.0):
5888            domain.quantity_statistics() # Run it silently
5889
5890        #--------------------------------------------------------------
5891        # Test extrema
5892        #--------------------------------------------------------------
5893        stage = domain.quantities_to_be_monitored['stage']
5894        assert stage['min'] <= stage['max']
5895
5896        assert num.allclose(stage['min'], initial_runup_height,
5897                            rtol=1.0/N)    # First order accuracy
5898
5899        depth = domain.quantities_to_be_monitored['stage-elevation']
5900        assert depth['min'] <= depth['max']
5901        assert depth['min'] >= 0.0
5902        assert depth['max'] >= 0.0
5903
5904        #--------------------------------------------------------------
5905        # Update boundary to allow inflow
5906        #--------------------------------------------------------------
5907        domain.set_boundary({'right': Bd})
5908
5909        #--------------------------------------------------------------
5910        # Evolve system through time
5911        #--------------------------------------------------------------
5912        for t in domain.evolve(yieldstep=0.1, finaltime=3.0):
5913            domain.quantity_statistics()    # Run it silently
5914
5915        #--------------------------------------------------------------
5916        # Test extrema again
5917        #--------------------------------------------------------------
5918        stage = domain.quantities_to_be_monitored['stage']
5919        assert stage['min'] <= stage['max']
5920
5921        assert num.allclose(stage['min'], initial_runup_height,
5922                            rtol = 1.0/N) # First order accuracy
5923
5924        depth = domain.quantities_to_be_monitored['stage-elevation']
5925        assert depth['min'] <= depth['max']
5926        assert depth['min'] >= 0.0
5927        assert depth['max'] >= 0.0
5928
5929        # Cleanup
5930        os.remove(domain.get_name() + '.' + domain.format)
5931
5932    def test_tight_slope_limiters(self):
5933        """Test that new slope limiters (Feb 2007) don't induce extremely
5934        small timesteps. This test actually reveals the problem as it
5935        was in March-April 2007
5936        """
5937        import time, os
5938        from Scientific.IO.NetCDF import NetCDFFile
5939        from data_manager import get_dataobject, extent_sww
5940        from mesh_factory import rectangular
5941
5942        # Create basic mesh
5943        points, vertices, boundary = rectangular(2, 2)
5944
5945        # Create shallow water domain
5946        domain = Domain(points, vertices, boundary)
5947        domain.default_order = 2
5948
5949        # This will pass
5950        #domain.tight_slope_limiters = 1
5951        #domain.H0 = 0.01
5952
5953        # This will fail
5954        #domain.tight_slope_limiters = 1
5955        #domain.H0 = 0.001
5956
5957        # This will pass provided C extension implements limiting of
5958        # momentum in _compute_speeds
5959        domain.tight_slope_limiters = 1
5960        domain.H0 = 0.001
5961        domain.protect_against_isolated_degenerate_timesteps = True
5962
5963        # Set some field values
5964        domain.set_quantity('elevation', lambda x,y: -x)
5965        domain.set_quantity('friction', 0.03)
5966
5967        # Boundary conditions
5968        B = Transmissive_boundary(domain)
5969        domain.set_boundary({'left': B, 'right': B, 'top': B, 'bottom': B})
5970
5971        # Initial condition - with jumps
5972        bed = domain.quantities['elevation'].vertex_values
5973        stage = num.zeros(bed.shape, num.float)
5974
5975        h = 0.3
5976        for i in range(stage.shape[0]):
5977            if i % 2 == 0:
5978                stage[i,:] = bed[i,:] + h
5979            else:
5980                stage[i,:] = bed[i,:]
5981
5982        domain.set_quantity('stage', stage)
5983
5984        domain.distribute_to_vertices_and_edges()
5985
5986        domain.set_name('tight_limiters')
5987        domain.format = 'sww'
5988        domain.smooth = True
5989        domain.reduction = mean
5990        domain.set_datadir('.')
5991        domain.smooth = False
5992        domain.store = True
5993
5994        # Evolution
5995        for t in domain.evolve(yieldstep=0.1, finaltime=0.3):
5996            #domain.write_time(track_speeds=True)
5997            stage = domain.quantities['stage'].vertex_values
5998
5999            # Get NetCDF
6000            fid = NetCDFFile(domain.writer.filename, netcdf_mode_r)
6001            stage_file = fid.variables['stage']
6002
6003            fid.close()
6004
6005        os.remove(domain.writer.filename)
6006
6007    def test_pmesh2Domain(self):
6008         import os
6009         import tempfile
6010
6011         fileName = tempfile.mktemp(".tsh")
6012         file = open(fileName, "w")
6013         file.write("4 3 # <vertex #> <x> <y> [attributes]\n \
60140 0.0 0.0 0.0 0.0 0.01 \n \
60151 1.0 0.0 10.0 10.0 0.02  \n \
60162 0.0 1.0 0.0 10.0 0.03  \n \
60173 0.5 0.25 8.0 12.0 0.04  \n \
6018# Vert att title  \n \
6019elevation  \n \
6020stage  \n \
6021friction  \n \
60222 # <triangle #> [<vertex #>] [<neigbouring triangle #>]  \n\
60230 0 3 2 -1  -1  1 dsg\n\
60241 0 1 3 -1  0 -1   ole nielsen\n\
60254 # <segment #> <vertex #>  <vertex #> [boundary tag] \n\
60260 1 0 2 \n\
60271 0 2 3 \n\
60282 2 3 \n\
60293 3 1 1 \n\
60303 0 # <x> <y> [attributes] ...Mesh Vertices... \n \
60310 216.0 -86.0 \n \
60321 160.0 -167.0 \n \
60332 114.0 -91.0 \n \
60343 # <vertex #>  <vertex #> [boundary tag] ...Mesh Segments... \n \
60350 0 1 0 \n \
60361 1 2 0 \n \
60372 2 0 0 \n \
60380 # <x> <y> ...Mesh Holes... \n \
60390 # <x> <y> <attribute>...Mesh Regions... \n \
60400 # <x> <y> <attribute>...Mesh Regions, area... \n\
6041#Geo reference \n \
604256 \n \
6043140 \n \
6044120 \n")
6045         file.close()
6046
6047         tags = {}
6048         b1 =  Dirichlet_boundary(conserved_quantities = num.array([0.0]))
6049         b2 =  Dirichlet_boundary(conserved_quantities = num.array([1.0]))
6050         b3 =  Dirichlet_boundary(conserved_quantities = num.array([2.0]))
6051         tags["1"] = b1
6052         tags["2"] = b2
6053         tags["3"] = b3
6054
6055         domain = Domain(mesh_filename=fileName)
6056                         # verbose=True, use_cache=True)
6057
6058         ## check the quantities
6059         answer = [[0., 8., 0.],
6060                   [0., 10., 8.]]
6061         assert num.allclose(domain.quantities['elevation'].vertex_values,
6062                             answer)
6063
6064         answer = [[0., 12., 10.],
6065                   [0., 10., 12.]]
6066         assert num.allclose(domain.quantities['stage'].vertex_values,
6067                             answer)
6068
6069         answer = [[0.01, 0.04, 0.03],
6070                   [0.01, 0.02, 0.04]]
6071         assert num.allclose(domain.quantities['friction'].vertex_values,
6072                             answer)
6073
6074         tagged_elements = domain.get_tagged_elements()
6075         assert num.allclose(tagged_elements['dsg'][0], 0)
6076         assert num.allclose(tagged_elements['ole nielsen'][0], 1)
6077
6078         msg = "test_tags_to_boundaries failed. Single boundary wasn't added."
6079         self.failUnless( domain.boundary[(1, 0)]  == '1', msg)
6080         self.failUnless( domain.boundary[(1, 2)]  == '2', msg)
6081         self.failUnless( domain.boundary[(0, 1)]  == '3', msg)
6082         self.failUnless( domain.boundary[(0, 0)]  == 'exterior', msg)
6083         msg = "test_pmesh2Domain Too many boundaries"
6084         self.failUnless( len(domain.boundary)  == 4, msg)
6085
6086         # FIXME change to use get_xllcorner
6087         msg = 'Bad geo-reference'
6088         self.failUnless(domain.geo_reference.xllcorner  == 140.0, msg)
6089
6090         domain = Domain(fileName)
6091
6092         answer = [[0., 8., 0.],
6093                   [0., 10., 8.]]
6094         assert num.allclose(domain.quantities['elevation'].vertex_values,
6095                             answer)
6096
6097         answer = [[0., 12., 10.],
6098                   [0., 10., 12.]]
6099         assert num.allclose(domain.quantities['stage'].vertex_values,
6100                             answer)
6101
6102         answer = [[0.01, 0.04, 0.03],
6103                   [0.01, 0.02, 0.04]]
6104         assert num.allclose(domain.quantities['friction'].vertex_values,
6105                             answer)
6106
6107         tagged_elements = domain.get_tagged_elements()
6108         assert num.allclose(tagged_elements['dsg'][0], 0)
6109         assert num.allclose(tagged_elements['ole nielsen'][0], 1)
6110
6111         msg = "test_tags_to_boundaries failed. Single boundary wasn't added."
6112         self.failUnless( domain.boundary[(1, 0)]  == '1', msg)
6113         self.failUnless( domain.boundary[(1, 2)]  == '2', msg)
6114         self.failUnless( domain.boundary[(0, 1)]  == '3', msg)
6115         self.failUnless( domain.boundary[(0, 0)]  == 'exterior', msg)
6116         msg = "test_pmesh2Domain Too many boundaries"
6117         self.failUnless( len(domain.boundary)  == 4, msg)
6118
6119         # FIXME change to use get_xllcorner
6120         msg = 'Bad geo_reference'
6121         self.failUnless(domain.geo_reference.xllcorner  == 140.0, msg)
6122
6123         os.remove(fileName)
6124
6125    def test_get_lone_vertices(self):
6126        a = [0.0, 0.0]
6127        b = [0.0, 2.0]
6128        c = [2.0, 0.0]
6129        d = [0.0, 4.0]
6130        e = [2.0, 2.0]
6131        f = [4.0, 0.0]
6132
6133        points = [a, b, c, d, e, f]
6134        #             bac,     bce,     ecf,     dbe
6135        vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4] ]
6136        boundary = {(0, 0): 'Third',
6137                    (0, 2): 'First',
6138                    (2, 0): 'Second',
6139                    (2, 1): 'Second',
6140                    (3, 1): 'Second',
6141                    (3, 2): 'Third'}
6142
6143        domain = Domain(points, vertices, boundary)
6144        domain.get_lone_vertices()
6145
6146    def test_fitting_using_shallow_water_domain(self):
6147        #Mesh in zone 56 (absolute coords)
6148
6149        x0 = 314036.58727982
6150        y0 = 6224951.2960092
6151
6152        a = [x0+0.0, y0+0.0]
6153        b = [x0+0.0, y0+2.0]
6154        c = [x0+2.0, y0+0.0]
6155        d = [x0+0.0, y0+4.0]
6156        e = [x0+2.0, y0+2.0]
6157        f = [x0+4.0, y0+0.0]
6158
6159        points = [a, b, c, d, e, f]
6160
6161        #             bac,     bce,     ecf,     dbe
6162        elements = [[1,0,2], [1,2,4], [4,2,5], [3,1,4] ]
6163
6164        # absolute going in ..
6165        mesh4 = Domain(points, elements, geo_reference=Geo_reference(56, 0, 0))
6166        mesh4.check_integrity()
6167        quantity = Quantity(mesh4)
6168
6169        # Get (enough) datapoints (relative to georef)
6170        data_points_rel = [[ 0.66666667, 0.66666667],
6171                           [ 1.33333333, 1.33333333],
6172                           [ 2.66666667, 0.66666667],
6173                           [ 0.66666667, 2.66666667],
6174                           [ 0.0, 1.0],
6175                           [ 0.0, 3.0],
6176                           [ 1.0, 0.0],
6177                           [ 1.0, 1.0],
6178                           [ 1.0, 2.0],
6179                           [ 1.0, 3.0],
6180                           [ 2.0, 1.0],
6181                           [ 3.0, 0.0],
6182                           [ 3.0, 1.0]]
6183
6184        data_geo_spatial = Geospatial_data(data_points_rel,
6185                                           geo_reference=Geo_reference(56,
6186                                                                       x0,
6187                                                                       y0))
6188        data_points_absolute = data_geo_spatial.get_data_points(absolute=True)
6189        attributes = linear_function(data_points_absolute)
6190        att = 'spam_and_eggs'
6191
6192        # Create .txt file
6193        ptsfile = tempfile.mktemp(".txt")
6194        file = open(ptsfile, "w")
6195        file.write(" x,y," + att + " \n")
6196        for data_point,attribute in map(None, data_points_absolute, attributes):
6197            row = (str(data_point[0]) + ',' +
6198                   str(data_point[1]) + ',' +
6199                   str(attribute))
6200            file.write(row + "\n")
6201        file.close()
6202
6203        # Check that values can be set from file
6204        quantity.set_values(filename=ptsfile, attribute_name=att, alpha=0)
6205        answer = linear_function(quantity.domain.get_vertex_coordinates())
6206
6207        assert num.allclose(quantity.vertex_values.flat, answer)
6208
6209        # Check that values can be set from file using default attribute
6210        quantity.set_values(filename = ptsfile, alpha = 0)
6211        assert num.allclose(quantity.vertex_values.flat, answer)
6212
6213        # Cleanup
6214        import os
6215        os.remove(ptsfile)
6216
6217    def test_fitting_example_that_crashed(self):
6218        """This unit test has been derived from a real world example
6219        (the Towradgi '98 rainstorm simulation).
6220
6221        It shows a condition where fitting as called from set_quantity crashes
6222        when ANUGA mesh is reused. The test passes in the case where a new mesh
6223        is created.
6224
6225        See ticket:314
6226        """
6227
6228        verbose = False
6229
6230        from anuga.shallow_water import Domain
6231        from anuga.pmesh.mesh_interface import create_mesh_from_regions
6232        from anuga.geospatial_data.geospatial_data import Geospatial_data
6233
6234        #--------------------------------------------------------------------
6235        # Create domain
6236        #--------------------------------------------------------------------
6237        W = 303400
6238        N = 6195800
6239        E = 308640
6240        S = 6193120
6241        bounding_polygon = [[W, S], [E, S], [E, N], [W, N]]
6242
6243        offending_regions = []
6244
6245        # From culvert 8
6246        offending_regions.append([[307611.43896231, 6193631.6894806],
6247                                  [307600.11394969, 6193608.2855474],
6248                                  [307597.41349586, 6193609.59227963],
6249                                  [307608.73850848, 6193632.99621282]])
6250        offending_regions.append([[307633.69143231, 6193620.9216536],
6251                                  [307622.36641969, 6193597.5177204],
6252                                  [307625.06687352, 6193596.21098818],
6253                                  [307636.39188614, 6193619.61492137]])
6254
6255        # From culvert 9
6256        offending_regions.append([[306326.69660524, 6194818.62900522],
6257                                  [306324.67939476, 6194804.37099478],
6258                                  [306323.75856492, 6194804.50127295],
6259                                  [306325.7757754, 6194818.7592834]])
6260        offending_regions.append([[306365.57160524, 6194813.12900522],
6261                                  [306363.55439476, 6194798.87099478],
6262                                  [306364.4752246, 6194798.7407166],
6263                                  [306366.49243508, 6194812.99872705]])
6264
6265        # From culvert 10
6266        offending_regions.append([[306955.071019428608, 6194465.704096679576],
6267                                  [306951.616980571358, 6194457.295903320424],
6268                                  [306950.044491164153, 6194457.941873183474],
6269                                  [306953.498530021403, 6194466.350066542625]])
6270        offending_regions.append([[307002.540019428649, 6194446.204096679576],
6271                                  [306999.085980571399, 6194437.795903320424],
6272                                  [307000.658469978604, 6194437.149933457375],
6273                                  [307004.112508835853, 6194445.558126816526]])
6274
6275        interior_regions = []
6276        for polygon in offending_regions:
6277            interior_regions.append( [polygon, 100] )
6278
6279        meshname = 'offending_mesh.msh'
6280        create_mesh_from_regions(bounding_polygon,
6281                                 boundary_tags={'south': [0], 'east': [1],
6282                                                'north': [2], 'west': [3]},
6283                                 maximum_triangle_area=1000000,
6284                                 interior_regions=interior_regions,
6285                                 filename=meshname,
6286                                 use_cache=False,
6287                                 verbose=verbose)
6288
6289        domain = Domain(meshname, use_cache=False, verbose=verbose)
6290
6291        #--------------------------------------------------------------------
6292        # Fit data point to mesh
6293        #--------------------------------------------------------------------
6294        points_file = 'offending_point.pts'
6295
6296        # This is the offending point
6297        G = Geospatial_data(data_points=[[306953.344, 6194461.5]],
6298                            attributes=[1])
6299        G.export_points_file(points_file)
6300
6301        domain.set_quantity('elevation', filename=points_file, use_cache=False,
6302                            verbose=verbose, alpha=0.01)
6303
6304if __name__ == "__main__":
6305    suite = unittest.makeSuite(Test_Shallow_Water, 'test')
6306    runner = unittest.TextTestRunner(verbosity=1)
6307    runner.run(suite)
Note: See TracBrowser for help on using the repository browser.