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

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

Played with culvert class and added update_interval

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