source: anuga_core/source/anuga/shallow_water/test_shallow_water_domain.py @ 6492

Last change on this file since 6492 was 6492, checked in by ole, 15 years ago

Reinstated compatibility code for Python2.4

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