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

Last change on this file since 5527 was 5450, checked in by ole, 16 years ago

Added exception is inlet_region for General flow does not contain any triangles.
An input check was also added.

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