source: inundation/ga/storm_surge/pyvolution/test_shallow_water.py @ 474

Last change on this file since 474 was 443, checked in by ole, 21 years ago

Work towards ensuring that this version (pyvolution mark 3) works as well as previous incarnation (pyvolution mark 2).
This is not true yet as some examples e.g. netherlands display 'blow-up'.
File size: 71.6 KB
Line 
1#!/usr/bin/env python
2
3import unittest
4from math import sqrt, pi
5
6from config import g, epsilon
7from Numeric import allclose, array, zeros, ones, Float
8from shallow_water import *
9
10       
11class TestCase(unittest.TestCase):
12    def setUp(self):
13        pass
14       
15    def tearDown(self):
16        pass
17
18    def test_rotate(self):
19        normal = array([0.0,-1.0])       
20
21        q = array([1.0,2.0,3.0])
22     
23        r = rotate(q, normal, direction = 1)
24        assert r[0] == 1
25        assert r[1] == -3
26        assert r[2] == 2
27
28        w = rotate(r, normal, direction = -1)       
29        assert allclose(w, q)
30
31
32
33    def test_flux_zero_case(self):
34        ql = zeros( 3, Float )
35        qr = zeros( 3, Float )
36        normal = zeros( 2, Float )
37        zl = zr = 0.
38        flux, max_speed = flux_function(normal, ql, qr, zl, zr)
39
40        assert allclose(flux, [0,0,0])
41        assert max_speed == 0.
42
43    def test_flux_constants(self):
44        w = 2.0
45       
46        normal = array([1.,0])       
47        ql = array([w, 0, 0])
48        qr = array([w, 0, 0])
49        zl = zr = 0.
50        h = w - (zl+zr)/2
51       
52        flux, max_speed = flux_function(normal, ql, qr, zl, zr)
53
54        assert allclose(flux, [0., 0.5*g*h**2, 0.])
55        assert max_speed == sqrt(g*h)
56       
57
58    def test_flux1(self):
59        #Use data from previous version of pyvolution
60        normal = array([1.,0])       
61        ql = array([-0.2, 2, 3])
62        qr = array([-0.2, 2, 3])
63        zl = zr = -0.5
64        flux, max_speed = flux_function(normal, ql, qr, zl, zr)
65
66        assert allclose(flux, [2.,13.77433333, 20.])
67        assert allclose(max_speed, 8.38130948661)
68
69
70    def test_flux2(self):
71        #Use data from previous version of pyvolution
72        normal = array([0., -1.])       
73        ql = array([-0.075, 2, 3])
74        qr = array([-0.075, 2, 3])
75        zl = zr = -0.375
76        flux, max_speed = flux_function(normal, ql, qr, zl, zr)
77
78        assert allclose(flux, [-3.,-20.0, -30.441])
79        assert allclose(max_speed, 11.7146428199)
80
81    def test_flux3(self):
82        #Use data from previous version of pyvolution       
83        normal = array([-sqrt(2)/2, sqrt(2)/2])       
84        ql = array([-0.075, 2, 3])
85        qr = array([-0.075, 2, 3])
86        zl = zr = -0.375
87        flux, max_speed = flux_function(normal, ql, qr, zl, zr)
88
89        assert allclose(flux, [sqrt(2)/2, 4.40221112, 7.3829019])
90        assert allclose(max_speed, 4.0716654239)
91
92    def test_flux4(self):
93        #Use data from previous version of pyvolution               
94        normal = array([-sqrt(2)/2, sqrt(2)/2])       
95        ql = array([-0.34319278, 0.10254161, 0.07273855])
96        qr = array([-0.30683287, 0.1071986, 0.05930515])
97        zl = zr = -0.375
98        flux, max_speed = flux_function(normal, ql, qr, zl, zr)
99
100        assert allclose(flux, [-0.04072676, -0.07096636, -0.01604364])
101        assert allclose(max_speed, 1.31414103233)                               
102
103    def test_sw_domain_simple(self):
104        a = [0.0, 0.0]
105        b = [0.0, 2.0]
106        c = [2.0,0.0]
107        d = [0.0, 4.0]
108        e = [2.0, 2.0]
109        f = [4.0,0.0]
110
111        points = [a, b, c, d, e, f]
112        #bac, bce, ecf, dbe, daf, dae
113        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4], [3,0,5], [3,0,4]]
114       
115        domain = Domain(points, vertices)       
116        domain.check_integrity()
117
118        for name in ['level', 'xmomentum', 'ymomentum',
119                     'elevation', 'friction']:
120            assert domain.quantities.has_key(name)
121
122
123        assert domain.get_conserved_quantities(0, edge=1) == 0.   
124
125
126    def test_boundary_conditions(self):
127       
128        a = [0.0, 0.0]
129        b = [0.0, 2.0]
130        c = [2.0,0.0]
131        d = [0.0, 4.0]
132        e = [2.0, 2.0]
133        f = [4.0,0.0]
134
135        points = [a, b, c, d, e, f]
136        #bac, bce, ecf, dbe
137        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4] ]
138        boundary = { (0, 0): 'Third',
139                     (0, 2): 'First',
140                     (2, 0): 'Second',
141                     (2, 1): 'Second',
142                     (3, 1): 'Second',
143                     (3, 2): 'Third'}                                         
144                     
145
146        domain = Domain(points, vertices, boundary)
147        domain.check_integrity()
148
149
150        domain.set_quantity('level', [[1,2,3], [5,5,5],
151                                      [0,0,9], [-6, 3, 3]])
152
153        domain.set_quantity('xmomentum', [[1,1,1], [2,2,2],
154                                          [3,3,3], [4, 4, 4]])
155
156        domain.set_quantity('ymomentum', [[10,10,10], [20,20,20],
157                                          [30,30,30], [40, 40, 40]])       
158
159
160        D = Dirichlet_boundary([5,2,1])
161        T = Transmissive_boundary(domain)
162        R = Reflective_boundary(domain)
163        domain.set_boundary( {'First': D, 'Second': T, 'Third': R})
164       
165        domain.update_boundary()
166
167        #Level
168        assert domain.quantities['level'].boundary_values[0] == 2.5
169        assert domain.quantities['level'].boundary_values[0] ==\
170               domain.get_conserved_quantities(0, edge=0)[0] #Reflective (2.5)
171        assert domain.quantities['level'].boundary_values[1] == 5. #Dirichlet
172        assert domain.quantities['level'].boundary_values[2] ==\
173               domain.get_conserved_quantities(2, edge=0)[0] #Transmissive (4.5)
174        assert domain.quantities['level'].boundary_values[3] ==\
175               domain.get_conserved_quantities(2, edge=1)[0] #Transmissive (4.5)
176        assert domain.quantities['level'].boundary_values[4] ==\
177               domain.get_conserved_quantities(3, edge=1)[0] #Transmissive (-1.5)
178        assert domain.quantities['level'].boundary_values[5] ==\
179               domain.get_conserved_quantities(3, edge=2)[0] #Reflective (-1.5)
180
181        #Xmomentum
182        assert domain.quantities['xmomentum'].boundary_values[0] == 1.0 #Reflective
183        assert domain.quantities['xmomentum'].boundary_values[1] == 2. #Dirichlet
184        assert domain.quantities['xmomentum'].boundary_values[2] ==\
185               domain.get_conserved_quantities(2, edge=0)[1] #Transmissive
186        assert domain.quantities['xmomentum'].boundary_values[3] ==\
187               domain.get_conserved_quantities(2, edge=1)[1] #Transmissive
188        assert domain.quantities['xmomentum'].boundary_values[4] ==\
189               domain.get_conserved_quantities(3, edge=1)[1] #Transmissive
190        assert domain.quantities['xmomentum'].boundary_values[5] == -4.0  #Reflective
191       
192        #Ymomentum
193        assert domain.quantities['ymomentum'].boundary_values[0] == -10.0 #Reflective
194        assert domain.quantities['ymomentum'].boundary_values[1] == 1.  #Dirichlet
195        assert domain.quantities['ymomentum'].boundary_values[2] == 30. #Transmissive
196        assert domain.quantities['ymomentum'].boundary_values[3] == 30. #Transmissive
197        assert domain.quantities['ymomentum'].boundary_values[4] == 40. #Transmissive
198        assert domain.quantities['ymomentum'].boundary_values[5] == 40. #Reflective
199
200
201       
202
203    def test_compute_fluxes0(self):
204        #Do a full triangle and check that fluxes cancel out for
205        #the constant level case
206       
207        a = [0.0, 0.0]
208        b = [0.0, 2.0]
209        c = [2.0,0.0]
210        d = [0.0, 4.0]
211        e = [2.0, 2.0]
212        f = [4.0,0.0]
213
214        points = [a, b, c, d, e, f]
215        #bac, bce, ecf, dbe
216        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
217       
218        domain = Domain(points, vertices)       
219        domain.set_quantity('level', [[2,2,2], [2,2,2],
220                                      [2,2,2], [2,2,2]])
221        domain.check_integrity()
222
223        zl=zr=0. #Assume flat bed
224       
225        #Flux across right edge of volume 1
226        normal = domain.get_normal(1,0)       
227        ql = domain.get_conserved_quantities(vol_id=1, edge=0)
228        qr = domain.get_conserved_quantities(vol_id=2, edge=2)
229        flux0, max_speed = flux_function(normal, ql, qr, zl, zr)
230       
231        #Check that flux seen from other triangles is inverse
232        tmp = qr; qr=ql; ql=tmp
233        normal = domain.get_normal(2,2)
234        flux, max_speed = flux_function(normal, ql, qr, zl, zr)       
235        assert allclose(flux + flux0, 0.)
236   
237        #Flux across upper edge of volume 1
238        normal = domain.get_normal(1,1)       
239        ql = domain.get_conserved_quantities(vol_id=1, edge=1)
240        qr = domain.get_conserved_quantities(vol_id=3, edge=0)
241        flux1, max_speed = flux_function(normal, ql, qr, zl, zr)
242       
243        #Check that flux seen from other triangles is inverse
244        tmp = qr; qr=ql; ql=tmp
245        normal = domain.get_normal(3,0)
246        flux, max_speed = flux_function(normal, ql, qr, zl, zr)       
247        assert allclose(flux + flux1, 0.)
248       
249        #Flux across lower left hypotenuse of volume 1
250        normal = domain.get_normal(1,2)       
251        ql = domain.get_conserved_quantities(vol_id=1, edge=2)
252        qr = domain.get_conserved_quantities(vol_id=0, edge=1)
253        flux2, max_speed = flux_function(normal, ql, qr, zl, zr)
254       
255        #Check that flux seen from other triangles is inverse
256        tmp = qr; qr=ql; ql=tmp
257        normal = domain.get_normal(0,1)
258        flux, max_speed = flux_function(normal, ql, qr, zl, zr)       
259        assert allclose(flux + flux2, 0.)
260
261
262        #Scale by edgelengths, add up anc check that total flux is zero
263        e0 = domain.edgelengths[1, 0]
264        e1 = domain.edgelengths[1, 1]
265        e2 = domain.edgelengths[1, 2]
266       
267        assert allclose(e0*flux0+e1*flux1+e2*flux2, 0.)
268
269        #Now check that compute_flux yields zeros as well
270        domain.compute_fluxes()
271
272        for name in ['level', 'xmomentum', 'ymomentum']:
273            #print name, domain.quantities[name].explicit_update
274            assert allclose(domain.quantities[name].explicit_update[1], 0)
275
276
277
278    def test_compute_fluxes1(self):
279        #Use values from previous version
280       
281        a = [0.0, 0.0]
282        b = [0.0, 2.0]
283        c = [2.0,0.0]
284        d = [0.0, 4.0]
285        e = [2.0, 2.0]
286        f = [4.0,0.0]
287
288        points = [a, b, c, d, e, f]
289        #bac, bce, ecf, dbe
290        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
291       
292        domain = Domain(points, vertices)
293        val0 = 2.+2.0/3
294        val1 = 4.+4.0/3
295        val2 = 8.+2.0/3
296        val3 = 2.+8.0/3
297       
298        domain.set_quantity('level', [[val0, val0, val0], [val1, val1, val1],
299                                      [val2, val2, val2], [val3, val3, val3]])
300        domain.check_integrity()
301
302        zl=zr=0. #Assume flat bed
303
304        #Flux across right edge of volume 1
305        normal = domain.get_normal(1,0)       
306        ql = domain.get_conserved_quantities(vol_id=1, edge=0)
307        qr = domain.get_conserved_quantities(vol_id=2, edge=2)
308        flux0, max_speed = flux_function(normal, ql, qr, zl, zr)
309        assert allclose(flux0, [-15.3598804, 253.71111111, 0.])
310        assert allclose(max_speed, 9.21592824046)
311
312
313        #Flux across upper edge of volume 1
314        normal = domain.get_normal(1,1)       
315        ql = domain.get_conserved_quantities(vol_id=1, edge=1)
316        qr = domain.get_conserved_quantities(vol_id=3, edge=0)
317        flux1, max_speed = flux_function(normal, ql, qr, zl, zr)
318        assert allclose(flux1, [2.4098563, 0., 123.04444444])
319        assert allclose(max_speed, 7.22956891292)
320
321        #Flux across lower left hypotenuse of volume 1
322        normal = domain.get_normal(1,2)       
323        ql = domain.get_conserved_quantities(vol_id=1, edge=2)
324        qr = domain.get_conserved_quantities(vol_id=0, edge=1)
325        flux2, max_speed = flux_function(normal, ql, qr, zl, zr)
326
327        assert allclose(flux2, [9.63942522, -61.59685738, -61.59685738])
328        assert allclose(max_speed, 7.22956891292)
329       
330        #Scale, add up and check that compute_fluxes is correct for vol 1
331        e0 = domain.edgelengths[1, 0]
332        e1 = domain.edgelengths[1, 1]
333        e2 = domain.edgelengths[1, 2]
334
335        total_flux = -(e0*flux0+e1*flux1+e2*flux2)/domain.areas[1]
336        assert allclose(total_flux, [-0.68218178, -166.6, -35.93333333])
337
338
339        domain.compute_fluxes()
340
341        #assert allclose(total_flux, domain.explicit_update[1,:])
342        for i, name in enumerate(['level', 'xmomentum', 'ymomentum']):
343            assert allclose(total_flux[i],
344                            domain.quantities[name].explicit_update[1])
345
346        #assert allclose(domain.explicit_update, [
347        #    [0., -69.68888889, -69.68888889],
348        #    [-0.68218178, -166.6, -35.93333333],
349        #    [-111.77316251, 69.68888889, 0.],
350        #    [-35.68522449, 0., 69.68888889]])
351
352        assert allclose(domain.quantities['level'].explicit_update,
353                        [0., -0.68218178, -111.77316251, -35.68522449])
354        assert allclose(domain.quantities['xmomentum'].explicit_update,
355                        [-69.68888889, -166.6, 69.68888889, 0])
356        assert allclose(domain.quantities['ymomentum'].explicit_update,
357                        [-69.68888889, -35.93333333, 0., 69.68888889])       
358
359       
360        #assert allclose(domain.quantities[name].explicit_update
361
362       
363       
364
365
366    def test_compute_fluxes2(self):
367        #Random values, incl momentum
368       
369        a = [0.0, 0.0]
370        b = [0.0, 2.0]
371        c = [2.0,0.0]
372        d = [0.0, 4.0]
373        e = [2.0, 2.0]
374        f = [4.0,0.0]
375
376        points = [a, b, c, d, e, f]
377        #bac, bce, ecf, dbe
378        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
379       
380        domain = Domain(points, vertices)
381        val0 = 2.+2.0/3
382        val1 = 4.+4.0/3
383        val2 = 8.+2.0/3
384        val3 = 2.+8.0/3
385
386        zl=zr=0 #Assume flat zero bed
387
388        domain.set_quantity('elevation', zl*ones( (4,3) ))
389       
390       
391        domain.set_quantity('level', [[val0, val0-1, val0-2],
392                                      [val1, val1+1, val1],
393                                      [val2, val2-2, val2],
394                                      [val3-0.5, val3, val3]])
395
396        domain.set_quantity('xmomentum',
397                            [[1, 2, 3], [3, 4, 5],
398                             [1, -1, 0], [0, -2, 2]])
399
400        domain.set_quantity('ymomentum',
401                            [[1, -1, 0], [0, -3, 2],
402                             [0, 1, 0], [-1, 2, 2]])               
403
404       
405        domain.check_integrity()
406
407
408
409        #Flux across right edge of volume 1
410        normal = domain.get_normal(1,0)       
411        ql = domain.get_conserved_quantities(vol_id=1, edge=0)
412        qr = domain.get_conserved_quantities(vol_id=2, edge=2)
413        flux0, max_speed = flux_function(normal, ql, qr, zl, zr)
414
415        #Flux across upper edge of volume 1
416        normal = domain.get_normal(1,1)       
417        ql = domain.get_conserved_quantities(vol_id=1, edge=1)
418        qr = domain.get_conserved_quantities(vol_id=3, edge=0)
419        flux1, max_speed = flux_function(normal, ql, qr, zl, zr)
420
421        #Flux across lower left hypotenuse of volume 1
422        normal = domain.get_normal(1,2)       
423        ql = domain.get_conserved_quantities(vol_id=1, edge=2)
424        qr = domain.get_conserved_quantities(vol_id=0, edge=1)
425        flux2, max_speed = flux_function(normal, ql, qr, zl, zr)
426       
427        #Scale, add up and check that compute_fluxes is correct for vol 1
428        e0 = domain.edgelengths[1, 0]
429        e1 = domain.edgelengths[1, 1]
430        e2 = domain.edgelengths[1, 2]
431
432        total_flux = -(e0*flux0+e1*flux1+e2*flux2)/domain.areas[1]
433
434       
435        domain.compute_fluxes()
436        for i, name in enumerate(['level', 'xmomentum', 'ymomentum']):
437            assert allclose(total_flux[i],
438                            domain.quantities[name].explicit_update[1])       
439        #assert allclose(total_flux, domain.explicit_update[1,:])
440       
441
442    def test_compute_fluxes3(self):
443        #Random values, incl momentum
444       
445        a = [0.0, 0.0]
446        b = [0.0, 2.0]
447        c = [2.0,0.0]
448        d = [0.0, 4.0]
449        e = [2.0, 2.0]
450        f = [4.0,0.0]
451
452        points = [a, b, c, d, e, f]
453        #bac, bce, ecf, dbe
454        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
455       
456        domain = Domain(points, vertices)
457        val0 = 2.+2.0/3
458        val1 = 4.+4.0/3
459        val2 = 8.+2.0/3
460        val3 = 2.+8.0/3
461
462        zl=zr=-3.75 #Assume constant bed (must be less than level)
463        domain.set_quantity('elevation', zl*ones( (4,3) ))
464       
465       
466        domain.set_quantity('level', [[val0, val0-1, val0-2],
467                                      [val1, val1+1, val1],
468                                      [val2, val2-2, val2],
469                                      [val3-0.5, val3, val3]])
470
471        domain.set_quantity('xmomentum',
472                            [[1, 2, 3], [3, 4, 5],
473                             [1, -1, 0], [0, -2, 2]])
474
475        domain.set_quantity('ymomentum',
476                            [[1, -1, 0], [0, -3, 2],
477                             [0, 1, 0], [-1, 2, 2]])               
478
479       
480        domain.check_integrity()
481
482
483
484        #Flux across right edge of volume 1
485        normal = domain.get_normal(1,0)       
486        ql = domain.get_conserved_quantities(vol_id=1, edge=0)
487        qr = domain.get_conserved_quantities(vol_id=2, edge=2)
488        flux0, max_speed = flux_function(normal, ql, qr, zl, zr)
489
490        #Flux across upper edge of volume 1
491        normal = domain.get_normal(1,1)       
492        ql = domain.get_conserved_quantities(vol_id=1, edge=1)
493        qr = domain.get_conserved_quantities(vol_id=3, edge=0)
494        flux1, max_speed = flux_function(normal, ql, qr, zl, zr)
495
496        #Flux across lower left hypotenuse of volume 1
497        normal = domain.get_normal(1,2)       
498        ql = domain.get_conserved_quantities(vol_id=1, edge=2)
499        qr = domain.get_conserved_quantities(vol_id=0, edge=1)
500        flux2, max_speed = flux_function(normal, ql, qr, zl, zr)
501       
502        #Scale, add up and check that compute_fluxes is correct for vol 1
503        e0 = domain.edgelengths[1, 0]
504        e1 = domain.edgelengths[1, 1]
505        e2 = domain.edgelengths[1, 2]
506
507        total_flux = -(e0*flux0+e1*flux1+e2*flux2)/domain.areas[1]
508       
509        domain.compute_fluxes()
510        for i, name in enumerate(['level', 'xmomentum', 'ymomentum']):
511            assert allclose(total_flux[i],
512                            domain.quantities[name].explicit_update[1])
513
514
515
516    def test_catching_negative_heights(self):
517       
518        a = [0.0, 0.0]
519        b = [0.0, 2.0]
520        c = [2.0,0.0]
521        d = [0.0, 4.0]
522        e = [2.0, 2.0]
523        f = [4.0,0.0]
524
525        points = [a, b, c, d, e, f]
526        #bac, bce, ecf, dbe
527        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
528       
529        domain = Domain(points, vertices)
530        val0 = 2.+2.0/3
531        val1 = 4.+4.0/3
532        val2 = 8.+2.0/3
533        val3 = 2.+8.0/3
534
535        zl=zr=4  #Too large
536        domain.set_quantity('elevation', zl*ones( (4,3) ))
537        domain.set_quantity('level', [[val0, val0-1, val0-2],
538                                      [val1, val1+1, val1],
539                                      [val2, val2-2, val2],
540                                      [val3-0.5, val3, val3]])
541
542        #Should fail
543        try:
544            domain.check_integrity()
545        except:
546            pass
547
548    #####################################################   
549    def test_gravity(self):
550        #Assuming no friction
551
552        from config import g
553       
554        a = [0.0, 0.0]
555        b = [0.0, 2.0]
556        c = [2.0, 0.0]
557        d = [0.0, 4.0]
558        e = [2.0, 2.0]
559        f = [4.0, 0.0]
560
561        points = [a, b, c, d, e, f]
562        #bac, bce, ecf, dbe
563        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
564       
565        domain = Domain(points, vertices)
566
567        #Set up for a gradient of (3,0) at mid triangle         
568        def slope(x, y):
569            return 3*x
570
571        h = 0.1
572        def level(x,y):
573            return slope(x,y)+h
574
575        domain.set_quantity('elevation', slope)
576        domain.set_quantity('level', level)
577
578        for name in domain.conserved_quantities:
579            assert allclose(domain.quantities[name].explicit_update, 0)
580            assert allclose(domain.quantities[name].semi_implicit_update, 0)
581           
582        domain.compute_forcing_terms()                                   
583
584        assert allclose(domain.quantities['level'].explicit_update, 0)
585        assert allclose(domain.quantities['xmomentum'].explicit_update, -g*h*3)
586        assert allclose(domain.quantities['ymomentum'].explicit_update, 0)
587       
588       
589
590    #####################################################
591    def test_first_order_extrapolator_const_z(self):
592       
593        a = [0.0, 0.0]
594        b = [0.0, 2.0]
595        c = [2.0, 0.0]
596        d = [0.0, 4.0]
597        e = [2.0, 2.0]
598        f = [4.0, 0.0]
599
600        points = [a, b, c, d, e, f]
601        #bac, bce, ecf, dbe
602        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
603       
604        domain = Domain(points, vertices)
605        val0 = 2.+2.0/3
606        val1 = 4.+4.0/3
607        val2 = 8.+2.0/3
608        val3 = 2.+8.0/3
609
610        zl=zr=-3.75 #Assume constant bed (must be less than level)       
611        domain.set_quantity('elevation', zl*ones( (4,3) ))
612        domain.set_quantity('level', [[val0, val0-1, val0-2],
613                                      [val1, val1+1, val1],
614                                      [val2, val2-2, val2],
615                                      [val3-0.5, val3, val3]])
616
617
618
619        domain.order = 1
620        domain.distribute_to_vertices_and_edges()
621
622        #Check that centroid values were distributed to vertices
623        C = domain.quantities['level'].centroid_values
624        for i in range(3): 
625            assert allclose( domain.quantities['level'].vertex_values[:,i], C)
626       
627
628    def test_first_order_limiter_variable_z(self):
629        #Check that first order limiter follows bed_slope
630        from Numeric import alltrue, greater_equal
631        from config import epsilon
632       
633        a = [0.0, 0.0]
634        b = [0.0, 2.0]
635        c = [2.0,0.0]
636        d = [0.0, 4.0]
637        e = [2.0, 2.0]
638        f = [4.0,0.0]
639
640        points = [a, b, c, d, e, f]
641        #bac, bce, ecf, dbe
642        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
643       
644        domain = Domain(points, vertices)
645        val0 = 2.+2.0/3
646        val1 = 4.+4.0/3
647        val2 = 8.+2.0/3
648        val3 = 2.+8.0/3
649
650        domain.set_quantity('elevation', [[0,0,0], [6,0,0],
651                                          [6,6,6], [6,6,6]])
652        domain.set_quantity('level', [[val0, val0, val0],
653                                      [val1, val1, val1],
654                                      [val2, val2, val2],
655                                      [val3, val3, val3]])
656
657        E = domain.quantities['elevation'].vertex_values
658        L = domain.quantities['level'].vertex_values       
659
660       
661        #Check that some levels are not above elevation (within eps)
662        #- so that the limiter has something to work with
663        assert not alltrue(alltrue(greater_equal(L,E-epsilon)))               
664
665        domain.order = 1
666        domain.distribute_to_vertices_and_edges()                           
667
668        #Check that all levels are above elevation (within eps)
669        assert alltrue(alltrue(greater_equal(L,E-epsilon)))
670
671
672    #####################################################   
673    def test_distribute_basic(self):
674        #Using test data generated by pyvolution-2
675        #Assuming no friction and flat bed (0.0)
676       
677        a = [0.0, 0.0]
678        b = [0.0, 2.0]
679        c = [2.0, 0.0]
680        d = [0.0, 4.0]
681        e = [2.0, 2.0]
682        f = [4.0, 0.0]
683
684        points = [a, b, c, d, e, f]
685        #bac, bce, ecf, dbe
686        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
687       
688        domain = Domain(points, vertices)
689       
690        val0 = 2.
691        val1 = 4.
692        val2 = 8.
693        val3 = 2.
694
695        domain.set_quantity('level', [val0, val1, val2, val3], 'centroids')
696        L = domain.quantities['level'].vertex_values
697       
698        #First order
699        domain.order = 1       
700        domain.distribute_to_vertices_and_edges()
701        assert allclose(L[1], val1)
702       
703        #Second order
704        domain.order = 2       
705        domain.distribute_to_vertices_and_edges()
706        assert allclose(L[1], [2.2, 4.9, 4.9])
707               
708
709       
710    def test_distribute_away_from_bed(self):           
711        #Using test data generated by pyvolution-2
712        #Assuming no friction and flat bed (0.0)
713       
714        a = [0.0, 0.0]
715        b = [0.0, 2.0]
716        c = [2.0, 0.0]
717        d = [0.0, 4.0]
718        e = [2.0, 2.0]
719        f = [4.0, 0.0]
720
721        points = [a, b, c, d, e, f]
722        #bac, bce, ecf, dbe
723        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
724       
725        domain = Domain(points, vertices)
726        L = domain.quantities['level'].vertex_values   
727       
728        def level(x,y):
729            return x**2
730       
731        domain.set_quantity('level', level, 'centroids')
732       
733        a, b = domain.quantities['level'].compute_gradients()           
734        assert allclose(a[1], 3.33333334)
735        assert allclose(b[1], 0.0)
736       
737        domain.order = 1
738        domain.distribute_to_vertices_and_edges()
739        assert allclose(L[1], 1.77777778)
740       
741        domain.order = 2
742        domain.distribute_to_vertices_and_edges()
743        assert allclose(L[1], [0.57777777, 2.37777778, 2.37777778])     
744
745
746       
747    def test_distribute_away_from_bed1(self):           
748        #Using test data generated by pyvolution-2
749        #Assuming no friction and flat bed (0.0)
750       
751        a = [0.0, 0.0]
752        b = [0.0, 2.0]
753        c = [2.0, 0.0]
754        d = [0.0, 4.0]
755        e = [2.0, 2.0]
756        f = [4.0, 0.0]
757
758        points = [a, b, c, d, e, f]
759        #bac, bce, ecf, dbe
760        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
761       
762        domain = Domain(points, vertices)
763        L = domain.quantities['level'].vertex_values   
764       
765        def level(x,y):
766            return x**4+y**2
767       
768        domain.set_quantity('level', level, 'centroids')
769        #print domain.quantities['level'].centroid_values               
770       
771        a, b = domain.quantities['level'].compute_gradients()           
772        assert allclose(a[1], 25.18518519)
773        assert allclose(b[1], 3.33333333)
774       
775        domain.order = 1
776        domain.distribute_to_vertices_and_edges()
777        assert allclose(L[1], 4.9382716)
778       
779        domain.order = 2
780        domain.distribute_to_vertices_and_edges()
781        assert allclose(L[1], [1.07160494, 6.46058131, 7.28262855])     
782
783   
784   
785    def test_distribute_near_bed(self):
786        #Using test data generated by pyvolution-2
787        #Assuming no friction and flat bed (0.0)
788
789        a = [0.0, 0.0]
790        b = [0.0, 2.0]
791        c = [2.0, 0.0]
792        d = [0.0, 4.0]
793        e = [2.0, 2.0]
794        f = [4.0, 0.0]
795
796        points = [a, b, c, d, e, f]
797        #bac, bce, ecf, dbe
798        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
799       
800        domain = Domain(points, vertices)
801
802
803        #Set up for a gradient of (3,0) at mid triangle         
804        def slope(x, y):
805            return 10*x
806
807        h = 0.1
808        def level(x,y):
809            return slope(x,y)+h
810
811        domain.set_quantity('elevation', slope)
812        domain.set_quantity('level', level, 'centroids')
813
814        #print domain.quantities['elevation'].centroid_values   
815        #print domain.quantities['level'].centroid_values
816       
817        E = domain.quantities['elevation'].vertex_values
818        L = domain.quantities['level'].vertex_values
819
820        #print E
821        domain.order = 1       
822        domain.distribute_to_vertices_and_edges()
823        ##assert allclose(L[1], [0.19999999, 20.05, 20.05])
824        assert allclose(L[1], [0.1, 20.1, 20.1])
825       
826        domain.order = 2       
827        domain.distribute_to_vertices_and_edges()
828        assert allclose(L[1], [0.1, 20.1, 20.1])
829
830    def test_distribute_near_bed1(self):
831        #Using test data generated by pyvolution-2
832        #Assuming no friction and flat bed (0.0)
833
834        a = [0.0, 0.0]
835        b = [0.0, 2.0]
836        c = [2.0, 0.0]
837        d = [0.0, 4.0]
838        e = [2.0, 2.0]
839        f = [4.0, 0.0]
840
841        points = [a, b, c, d, e, f]
842        #bac, bce, ecf, dbe
843        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
844       
845        domain = Domain(points, vertices)
846
847
848        #Set up for a gradient of (3,0) at mid triangle         
849        def slope(x, y):
850            return x**4+y**2
851
852        h = 0.1
853        def level(x,y):
854            return slope(x,y)+h
855
856        domain.set_quantity('elevation', slope)
857        domain.set_quantity('level', level)
858
859        #print domain.quantities['elevation'].centroid_values   
860        #print domain.quantities['level'].centroid_values
861       
862        E = domain.quantities['elevation'].vertex_values
863        L = domain.quantities['level'].vertex_values
864
865        #print E
866        domain.order = 1
867        domain.distribute_to_vertices_and_edges()
868        ##assert allclose(L[1], [4.19999999, 16.07142857, 20.02857143])
869        assert allclose(L[1], [4.1, 16.1, 20.1])       
870
871        domain.order = 2       
872        domain.distribute_to_vertices_and_edges()
873        assert allclose(L[1], [4.1, 16.1, 20.1])
874
875 
876
877    def test_second_order_distribute_real_data(self):
878        #Using test data generated by pyvolution-2
879        #Assuming no friction and flat bed (0.0)
880
881        a = [0.0, 0.0]
882        b = [0.0, 1.0/5]
883        c = [0.0, 2.0/5]
884        d = [1.0/5, 0.0]
885        e = [1.0/5, 1.0/5]
886        f = [1.0/5, 2.0/5]
887        g = [2.0/5, 2.0/5]
888
889        points = [a, b, c, d, e, f, g]
890        #bae, efb, cbf, feg
891        vertices = [ [1,0,4], [4,5,1], [2,1,5], [5,4,6]]
892       
893        domain = Domain(points, vertices)
894
895        def slope(x, y):
896            return -x/3
897
898        domain.set_quantity('elevation', slope)
899        domain.set_quantity('level', 
900                            [0.01298164, 0.00365611, 0.01440365, -0.0381856437096], 
901                            'centroids')
902        domain.set_quantity('xmomentum', 
903                            [0.00670439, 0.01263789, 0.00647805, 0.0178180740668],
904                            'centroids')
905        domain.set_quantity('ymomentum',                           
906                            [-7.23510980e-004, -6.30413883e-005, 6.30413883e-005, 0.000200907255866],
907                            'centroids')                           
908       
909        E = domain.quantities['elevation'].vertex_values
910        L = domain.quantities['level'].vertex_values
911        X = domain.quantities['xmomentum'].vertex_values       
912        Y = domain.quantities['ymomentum'].vertex_values               
913
914        #print E
915        domain.order = 2       
916        domain.distribute_to_vertices_and_edges()
917
918        #print L[1,:]
919        #print X[1,:]
920        #print Y[1,:]
921
922        assert allclose(L[1,:], [-0.00825735775384, -0.00801881482869, 0.0272445025825])
923        assert allclose(X[1,:], [0.0143507718962, 0.0142502147066, 0.00931268339717])
924        assert allclose(Y[1,:], [-0.000117062180693, 7.94434448109e-005, -0.000151505429018])
925       
926       
927    def test_second_order_flat_bed_onestep(self):
928
929        from mesh_factory import rectangular
930        from shallow_water import Domain, Reflective_boundary,\
931             Dirichlet_boundary, Constant_height
932        from Numeric import array
933
934        #Create basic mesh
935        points, vertices, boundary = rectangular(6, 6)
936
937        #Create shallow water domain
938        domain = Domain(points, vertices, boundary)
939        domain.smooth = False
940        domain.default_order=2
941
942        # Boundary conditions
943        Br = Reflective_boundary(domain)
944        Bd = Dirichlet_boundary([0.1, 0., 0.])
945        domain.set_boundary({'left': Bd, 'right': Br, 'top': Br, 'bottom': Br})
946
947        domain.check_integrity()
948
949        #Evolution
950        for t in domain.evolve(yieldstep = 0.05, finaltime = 0.05):
951            pass# domain.write_time()
952
953        #Data from earlier version of pyvolution
954        assert allclose(domain.min_timestep, 0.0396825396825)
955        assert allclose(domain.max_timestep, 0.0396825396825) 
956
957        assert allclose(domain.quantities['level'].centroid_values[:12],
958                        [0.00171396, 0.02656103, 0.00241523, 0.02656103,
959                        0.00241523, 0.02656103, 0.00241523, 0.02656103,
960                        0.00241523, 0.02656103, 0.00241523, 0.0272623])
961
962        domain.distribute_to_vertices_and_edges()
963        assert allclose(domain.quantities['level'].vertex_values[:12,0],
964                        [0.0001714, 0.02656103, 0.00024152,
965                        0.02656103, 0.00024152, 0.02656103,
966                        0.00024152, 0.02656103, 0.00024152,
967                        0.02656103, 0.00024152, 0.0272623])
968
969        assert allclose(domain.quantities['level'].vertex_values[:12,1],
970                        [0.00315012, 0.02656103, 0.00024152, 0.02656103,
971                         0.00024152, 0.02656103, 0.00024152, 0.02656103,
972                         0.00024152, 0.02656103, 0.00040506, 0.0272623])
973       
974        assert allclose(domain.quantities['level'].vertex_values[:12,2],
975                        [0.00182037, 0.02656103, 0.00676264,
976                         0.02656103, 0.00676264, 0.02656103,
977                         0.00676264, 0.02656103, 0.00676264,
978                         0.02656103, 0.0065991, 0.0272623])
979
980        assert allclose(domain.quantities['xmomentum'].centroid_values[:12],
981                        [0.00113961, 0.01302432, 0.00148672,
982                         0.01302432, 0.00148672, 0.01302432,
983                         0.00148672, 0.01302432, 0.00148672 ,
984                         0.01302432, 0.00148672, 0.01337143])
985
986        assert allclose(domain.quantities['ymomentum'].centroid_values[:12],
987                        [-2.91240050e-004 , 1.22721531e-004,
988                         -1.22721531e-004, 1.22721531e-004 ,
989                         -1.22721531e-004, 1.22721531e-004,
990                         -1.22721531e-004 , 1.22721531e-004,
991                         -1.22721531e-004, 1.22721531e-004,
992                         -1.22721531e-004, -4.57969873e-005])
993
994
995       
996    def test_second_order_flat_bed_moresteps(self):
997
998        from mesh_factory import rectangular
999        from shallow_water import Domain, Reflective_boundary,\
1000             Dirichlet_boundary, Constant_height
1001        from Numeric import array
1002
1003        #Create basic mesh
1004        points, vertices, boundary = rectangular(6, 6)
1005
1006        #Create shallow water domain
1007        domain = Domain(points, vertices, boundary)
1008        domain.smooth = False
1009        domain.default_order=2
1010
1011        # Boundary conditions
1012        Br = Reflective_boundary(domain)
1013        Bd = Dirichlet_boundary([0.1, 0., 0.])
1014        domain.set_boundary({'left': Bd, 'right': Br, 'top': Br, 'bottom': Br})
1015
1016        domain.check_integrity()
1017
1018        #Evolution
1019        for t in domain.evolve(yieldstep = 0.05, finaltime = 0.1):
1020            pass
1021
1022        #Data from earlier version of pyvolution
1023        #assert allclose(domain.min_timestep, 0.0396825396825)
1024        #assert allclose(domain.max_timestep, 0.0396825396825)
1025        #print domain.quantities['level'].centroid_values
1026
1027
1028    def test_flatbed_first_order(self):
1029        from mesh_factory import rectangular
1030        from shallow_water import Domain,\
1031             Reflective_boundary, Dirichlet_boundary
1032
1033        from Numeric import array
1034
1035        #Create basic mesh
1036        N = 8
1037        points, vertices, boundary = rectangular(N, N)
1038
1039        #Create shallow water domain
1040        domain = Domain(points, vertices, boundary)
1041        domain.smooth = False
1042        domain.visualise = False
1043        domain.default_order=1
1044
1045        # Boundary conditions
1046        Br = Reflective_boundary(domain)
1047        Bd = Dirichlet_boundary([0.2,0.,0.])
1048       
1049        domain.set_boundary({'left': Bd, 'right': Br, 'top': Br, 'bottom': Br})
1050        domain.check_integrity()
1051
1052 
1053        #Evolution
1054        for t in domain.evolve(yieldstep = 0.02, finaltime = 0.5):
1055            pass
1056            #domain.write_time()
1057
1058        #FIXME: These numbers were from version before 25/10       
1059        #assert allclose(domain.min_timestep, 0.0140413643926)
1060        #assert allclose(domain.max_timestep, 0.0140947355753)
1061
1062        for i in range(3):
1063            #assert allclose(domain.quantities['level'].edge_values[:4,i],
1064            #                [0.10730244,0.12337617,0.11200126,0.12605666])
1065
1066            assert allclose(domain.quantities['xmomentum'].edge_values[:4,i],
1067                            [0.07610894,0.06901572,0.07284461,0.06819712])
1068
1069            #assert allclose(domain.quantities['ymomentum'].edge_values[:4,i],
1070            #                [-0.0060238, -0.00157404, -0.00309633, -0.0001637])                         
1071       
1072    def test_flatbed_second_order(self):
1073        from mesh_factory import rectangular
1074        from shallow_water import Domain,\
1075             Reflective_boundary, Dirichlet_boundary
1076
1077        from Numeric import array
1078
1079        #Create basic mesh
1080        N = 8
1081        points, vertices, boundary = rectangular(N, N)
1082
1083        #Create shallow water domain
1084        domain = Domain(points, vertices, boundary)
1085        domain.smooth = False
1086        domain.visualise = False
1087        domain.default_order=2
1088        #domain.minimum_allowed_height = 0.0 #Makes it like the 'oldstyle' balance
1089
1090        # Boundary conditions
1091        Br = Reflective_boundary(domain)
1092        Bd = Dirichlet_boundary([0.2,0.,0.])
1093
1094        domain.set_boundary({'left': Bd, 'right': Br, 'top': Br, 'bottom': Br})
1095        domain.check_integrity()
1096 
1097        #Evolution
1098        for t in domain.evolve(yieldstep = 0.01, finaltime = 0.03):
1099            pass
1100
1101
1102        assert allclose(domain.min_timestep, 0.0210448446782)
1103        assert allclose(domain.max_timestep, 0.0210448446782) 
1104
1105        #print domain.quantities['level'].vertex_values[:4,0]
1106        #print domain.quantities['xmomentum'].vertex_values[:4,0]
1107        #print domain.quantities['ymomentum'].vertex_values[:4,0]
1108
1109        #FIXME: These numbers were from version before 25/10
1110        #assert allclose(domain.quantities['level'].vertex_values[:4,0],
1111        #                [0.00101913,0.05352143,0.00104852,0.05354394])
1112       
1113        assert allclose(domain.quantities['xmomentum'].vertex_values[:4,0],
1114                        [ 0.00064835, 0.03685719, 0.00085073, 0.03687313])
1115
1116        #assert allclose(domain.quantities['xmomentum'].vertex_values[:4,0],
1117        #                [0.00090581,0.03685719,0.00088303,0.03687313])
1118       
1119        assert allclose(domain.quantities['ymomentum'].vertex_values[:4,0],
1120                        [-0.00139463,0.0006156,-0.00060364,0.00061827])
1121
1122
1123       
1124    def test_flatbed_second_order_distribute(self):
1125        #Use real data from pyvolution 2
1126        #painfully setup and extracted.
1127        from mesh_factory import rectangular
1128        from shallow_water import Domain,\
1129             Reflective_boundary, Dirichlet_boundary
1130
1131        from Numeric import array
1132
1133        #Create basic mesh
1134        N = 8
1135        points, vertices, boundary = rectangular(N, N)
1136
1137        #Create shallow water domain
1138        domain = Domain(points, vertices, boundary)
1139        domain.smooth = False
1140        domain.visualise = False
1141        domain.default_order=domain.order=2
1142
1143        # Boundary conditions
1144        Br = Reflective_boundary(domain)
1145        Bd = Dirichlet_boundary([0.2,0.,0.])
1146
1147        domain.set_boundary({'left': Bd, 'right': Br, 'top': Br, 'bottom': Br})
1148        domain.check_integrity()
1149
1150
1151
1152        for V in [False, True]:
1153            if V:
1154                #Set centroids as if system had been evolved           
1155                L = zeros(2*N*N, Float)
1156                L[:32] = [7.21205592e-003, 5.35214298e-002, 1.00910824e-002,
1157                          5.35439433e-002, 1.00910824e-002, 5.35439433e-002,
1158                          1.00910824e-002, 5.35439433e-002, 1.00910824e-002,
1159                          5.35439433e-002, 1.00910824e-002, 5.35439433e-002,
1160                          1.00910824e-002, 5.35393928e-002, 1.02344264e-002,
1161                          5.59605058e-002, 0.00000000e+000, 3.31027800e-004,
1162                          0.00000000e+000, 4.37962142e-005, 0.00000000e+000,
1163                          4.37962142e-005, 0.00000000e+000, 4.37962142e-005,
1164                          0.00000000e+000, 4.37962142e-005, 0.00000000e+000,
1165                          4.37962142e-005, 0.00000000e+000, 4.37962142e-005,
1166                          0.00000000e+000, 5.57305948e-005]
1167       
1168                X = zeros(2*N*N, Float)
1169                X[:32] = [6.48351607e-003, 3.68571894e-002, 8.50733285e-003,
1170                          3.68731327e-002, 8.50733285e-003, 3.68731327e-002,
1171                          8.50733285e-003, 3.68731327e-002, 8.50733285e-003,
1172                          3.68731327e-002, 8.50733285e-003, 3.68731327e-002,
1173                          8.50733285e-003, 3.68693861e-002, 8.65220973e-003,
1174                          3.85055387e-002, 0.00000000e+000, 2.86060840e-004,
1175                          0.00000000e+000, 3.58905503e-005, 0.00000000e+000,
1176                          3.58905503e-005, 0.00000000e+000, 3.58905503e-005,
1177                          0.00000000e+000, 3.58905503e-005, 0.00000000e+000,
1178                          3.58905503e-005, 0.00000000e+000, 3.58905503e-005,
1179                          0.00000000e+000, 4.57662812e-005]
1180
1181                Y = zeros(2*N*N, Float)
1182                Y[:32]=[-1.39463104e-003, 6.15600298e-004, -6.03637382e-004,
1183                        6.18272251e-004, -6.03637382e-004, 6.18272251e-004,
1184                        -6.03637382e-004, 6.18272251e-004, -6.03637382e-004,
1185                        6.18272251e-004, -6.03637382e-004, 6.18272251e-004,
1186                        -6.03637382e-004, 6.18599320e-004, -6.74622797e-004,
1187                        -1.48934756e-004, 0.00000000e+000, -5.35079969e-005,
1188                        0.00000000e+000, -2.57264987e-005, 0.00000000e+000,
1189                        -2.57264987e-005, 0.00000000e+000, -2.57264987e-005,
1190                        0.00000000e+000, -2.57264987e-005, 0.00000000e+000,
1191                        -2.57264987e-005, 0.00000000e+000, -2.57264987e-005,
1192                        0.00000000e+000, -2.57635178e-005]
1193
1194       
1195                domain.set_quantity('level', L, 'centroids')
1196                domain.set_quantity('xmomentum', X, 'centroids')
1197                domain.set_quantity('ymomentum', Y, 'centroids')
1198       
1199                domain.check_integrity()
1200            else:   
1201                #Evolution
1202                for t in domain.evolve(yieldstep = 0.01, finaltime = 0.03):
1203                    pass
1204                assert allclose(domain.min_timestep, 0.0210448446782)
1205                assert allclose(domain.max_timestep, 0.0210448446782)
1206       
1207
1208            #Centroids were correct but not vertices.
1209            #Hence the check of distribute below.
1210            assert allclose(domain.quantities['level'].centroid_values[:4],
1211                            [0.00721206,0.05352143,0.01009108,0.05354394])
1212
1213            assert allclose(domain.quantities['xmomentum'].centroid_values[:4],
1214                            [0.00648352,0.03685719,0.00850733,0.03687313])
1215
1216            assert allclose(domain.quantities['ymomentum'].centroid_values[:4],
1217                            [-0.00139463,0.0006156,-0.00060364,0.00061827])
1218
1219            #print 'C17=', domain.quantities['xmomentum'].centroid_values[17]
1220            #print 'C19=', domain.quantities['xmomentum'].centroid_values[19]
1221
1222            #assert allclose(domain.quantities['xmomentum'].centroid_values[17],0.00028606084)
1223            ##print domain.quantities['xmomentum'].centroid_values[17], V
1224            ##print
1225            if not V:           
1226                assert allclose(domain.quantities['xmomentum'].centroid_values[17], 0.0)
1227            else: 
1228                assert allclose(domain.quantities['xmomentum'].centroid_values[17], 0.00028606084)           
1229
1230            import copy
1231            XX = copy.copy(domain.quantities['xmomentum'].centroid_values)
1232            assert allclose(domain.quantities['xmomentum'].centroid_values, XX)
1233       
1234            domain.distribute_to_vertices_and_edges()
1235
1236            #assert allclose(domain.quantities['xmomentum'].centroid_values, XX)
1237
1238            assert allclose(domain.quantities['xmomentum'].centroid_values[17],
1239                            0.0)
1240
1241
1242            #FIXME: These numbers were from version before 25/10
1243            #assert allclose(domain.quantities['level'].vertex_values[:4,0],
1244            #                [0.00101913,0.05352143,0.00104852,0.05354394])
1245
1246            assert allclose(domain.quantities['ymomentum'].vertex_values[:4,0],
1247                            [-0.00139463,0.0006156,-0.00060364,0.00061827])
1248
1249
1250            assert allclose(domain.quantities['xmomentum'].vertex_values[:4,0],
1251                            [0.00064835,0.03685719,0.00085073,0.03687313])
1252
1253
1254            #NB NO longer relvant:
1255
1256            #This was the culprit. First triangles vertex 0 had an
1257            #x-momentum of 0.0064835 instead of 0.00090581 and
1258            #third triangle had 0.00850733 instead of 0.00088303
1259            #print domain.quantities['xmomentum'].vertex_values[:4,0]
1260
1261            #print domain.quantities['xmomentum'].vertex_values[:4,0]
1262            #assert allclose(domain.quantities['xmomentum'].vertex_values[:4,0],
1263            #                [0.00090581,0.03685719,0.00088303,0.03687313])       
1264
1265
1266
1267       
1268       
1269    def test_bedslope_problem_first_order(self):
1270
1271        from mesh_factory import rectangular
1272        from shallow_water import Domain, Reflective_boundary, Constant_height
1273        from Numeric import array
1274
1275        #Create basic mesh
1276        points, vertices, boundary = rectangular(6, 6)
1277
1278        #Create shallow water domain
1279        domain = Domain(points, vertices, boundary)
1280        domain.smooth = False
1281        domain.default_order=1
1282
1283        #Bed-slope and friction
1284        def x_slope(x, y):
1285            return -x/3
1286
1287        domain.set_quantity('elevation', x_slope)
1288
1289        # Boundary conditions
1290        Br = Reflective_boundary(domain)
1291        domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
1292
1293        #Initial condition
1294        domain.set_quantity('level', Constant_height(x_slope, 0.05))
1295        domain.check_integrity()
1296
1297        #Evolution
1298        for t in domain.evolve(yieldstep = 0.05, finaltime = 0.05):
1299            pass# domain.write_time()
1300
1301        assert allclose(domain.min_timestep, 0.050010003001)
1302        assert allclose(domain.max_timestep, 0.050010003001)
1303
1304       
1305    def test_bedslope_problem_first_order_moresteps(self):
1306
1307        from mesh_factory import rectangular
1308        from shallow_water import Domain, Reflective_boundary, Constant_height
1309        from Numeric import array
1310
1311        #Create basic mesh
1312        points, vertices, boundary = rectangular(6, 6)
1313
1314        #Create shallow water domain
1315        domain = Domain(points, vertices, boundary)
1316        domain.smooth = False
1317        domain.default_order=1
1318
1319        #Bed-slope and friction
1320        def x_slope(x, y):
1321            return -x/3
1322
1323        domain.set_quantity('elevation', x_slope)
1324
1325        # Boundary conditions
1326        Br = Reflective_boundary(domain)
1327        domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
1328
1329        #Initial condition
1330        domain.set_quantity('level', Constant_height(x_slope, 0.05))
1331        domain.check_integrity()
1332
1333        #Evolution
1334        for t in domain.evolve(yieldstep = 0.05, finaltime = 0.5):
1335            pass# domain.write_time()
1336
1337        #Data from earlier version of pyvolution
1338        #print domain.quantities['level'].centroid_values
1339       
1340        assert allclose(domain.quantities['level'].centroid_values,
1341                        [-0.02998628, -0.01520652, -0.03043492,
1342                        -0.0149132, -0.03004706, -0.01476251,
1343                        -0.0298215, -0.01467976, -0.02988158,
1344                        -0.01474662, -0.03036161, -0.01442995,
1345                        -0.07624583, -0.06297061, -0.07733792,
1346                        -0.06342237, -0.07695439, -0.06289595,
1347                        -0.07635559, -0.0626065, -0.07633628,
1348                        -0.06280072, -0.07739632, -0.06386738,
1349                        -0.12161738, -0.11028239, -0.1223796,
1350                        -0.11095953, -0.12189744, -0.11048616,
1351                        -0.12074535, -0.10987605, -0.12014311,
1352                        -0.10976691, -0.12096859, -0.11087692,
1353                        -0.16868259, -0.15868061, -0.16801135,
1354                        -0.1588003, -0.16674343, -0.15813323,
1355                        -0.16457595, -0.15693826, -0.16281096,
1356                        -0.15585154, -0.16283873, -0.15540068,
1357                        -0.17450362, -0.19919913, -0.18062882,
1358                        -0.19764131, -0.17783111, -0.19407213,
1359                        -0.1736915, -0.19053624, -0.17228678,
1360                        -0.19105634, -0.17920133, -0.1968828,
1361                        -0.14244395, -0.14604641, -0.14473537,
1362                        -0.1506107, -0.14510055, -0.14919522,
1363                        -0.14175896, -0.14560798, -0.13911658,
1364                        -0.14439383, -0.13924047, -0.14829043])
1365                       
1366       
1367    def test_bedslope_problem_second_order_one_step(self):
1368
1369        from mesh_factory import rectangular
1370        from shallow_water import Domain, Reflective_boundary, Constant_height
1371        from Numeric import array
1372
1373        #Create basic mesh
1374        points, vertices, boundary = rectangular(6, 6)
1375
1376        #Create shallow water domain
1377        domain = Domain(points, vertices, boundary)
1378        domain.smooth = False
1379        domain.default_order=2
1380
1381        #Bed-slope and friction at vertices (and interpolated elsewhere)
1382        def x_slope(x, y):
1383            return -x/3
1384
1385        domain.set_quantity('elevation', x_slope)
1386
1387        # Boundary conditions
1388        Br = Reflective_boundary(domain)
1389        domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
1390
1391        #Initial condition
1392        domain.set_quantity('level', Constant_height(x_slope, 0.05))
1393        domain.check_integrity()
1394
1395        assert allclose(domain.quantities['level'].centroid_values,
1396                        [0.01296296, 0.03148148, 0.01296296,
1397                        0.03148148, 0.01296296, 0.03148148,
1398                        0.01296296, 0.03148148, 0.01296296,
1399                        0.03148148, 0.01296296, 0.03148148,
1400                        -0.04259259, -0.02407407, -0.04259259,
1401                        -0.02407407, -0.04259259, -0.02407407,
1402                        -0.04259259, -0.02407407, -0.04259259,
1403                        -0.02407407, -0.04259259, -0.02407407,
1404                        -0.09814815, -0.07962963, -0.09814815,
1405                        -0.07962963, -0.09814815, -0.07962963,
1406                        -0.09814815, -0.07962963, -0.09814815,
1407                        -0.07962963, -0.09814815, -0.07962963,
1408                        -0.1537037 , -0.13518519, -0.1537037,
1409                        -0.13518519, -0.1537037, -0.13518519,
1410                        -0.1537037 , -0.13518519, -0.1537037,
1411                        -0.13518519, -0.1537037, -0.13518519,
1412                        -0.20925926, -0.19074074, -0.20925926,
1413                        -0.19074074, -0.20925926, -0.19074074,
1414                        -0.20925926, -0.19074074, -0.20925926,
1415                        -0.19074074, -0.20925926, -0.19074074,
1416                        -0.26481481, -0.2462963, -0.26481481,
1417                        -0.2462963, -0.26481481, -0.2462963,
1418                        -0.26481481, -0.2462963, -0.26481481,
1419                        -0.2462963, -0.26481481, -0.2462963])
1420
1421
1422        #print domain.quantities['level'].extrapolate_second_order()
1423        #domain.distribute_to_vertices_and_edges()
1424        #print domain.quantities['level'].vertex_values[:,0]       
1425       
1426        #Evolution
1427        for t in domain.evolve(yieldstep = 0.05, finaltime = 0.05):
1428            #domain.write_time()           
1429            pass
1430
1431
1432        #print domain.quantities['level'].centroid_values
1433        assert allclose(domain.quantities['level'].centroid_values,
1434                 [0.01290985, 0.02356019, 0.01619096, 0.02356019, 0.01619096,
1435                  0.02356019, 0.01619096, 0.02356019, 0.01619096, 0.02356019,
1436                  0.01619096, 0.0268413, -0.04411074, -0.0248011, -0.04186556,
1437                  -0.0248011, -0.04186556, -0.0248011, -0.04186556, -0.0248011,
1438                  -0.04186556, -0.0248011, -0.04186556, -0.02255593,
1439                  -0.09966629, -0.08035666, -0.09742112, -0.08035666,
1440                  -0.09742112, -0.08035666, -0.09742112, -0.08035666,
1441                  -0.09742112, -0.08035666, -0.09742112, -0.07811149,
1442                  -0.15522185, -0.13591222, -0.15297667, -0.13591222,
1443                  -0.15297667, -0.13591222, -0.15297667, -0.13591222,
1444                  -0.15297667, -0.13591222, -0.15297667, -0.13366704,
1445                  -0.2107774, -0.19146777, -0.20853223, -0.19146777,
1446                  -0.20853223, -0.19146777, -0.20853223, -0.19146777,
1447                  -0.20853223, -0.19146777, -0.20853223, -0.1892226,
1448                  -0.26120669, -0.24776246, -0.25865535, -0.24776246,
1449                  -0.25865535, -0.24776246, -0.25865535, -0.24776246,
1450                  -0.25865535, -0.24776246, -0.25865535, -0.24521113])
1451
1452       
1453
1454    def test_bedslope_problem_second_order_two_steps(self):
1455
1456        from mesh_factory import rectangular
1457        from shallow_water import Domain, Reflective_boundary, Constant_height
1458        from Numeric import array
1459
1460        #Create basic mesh
1461        points, vertices, boundary = rectangular(6, 6)
1462
1463        #Create shallow water domain
1464        domain = Domain(points, vertices, boundary)
1465        domain.smooth = False
1466        domain.default_order=2
1467
1468        #Bed-slope and friction at vertices (and interpolated elsewhere)
1469        def x_slope(x, y):
1470            return -x/3
1471
1472        domain.set_quantity('elevation', x_slope)
1473
1474        # Boundary conditions
1475        Br = Reflective_boundary(domain)
1476        domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
1477
1478        #Initial condition
1479        domain.set_quantity('level', Constant_height(x_slope, 0.05))
1480        domain.check_integrity()
1481
1482        assert allclose(domain.quantities['level'].centroid_values,
1483                        [0.01296296, 0.03148148, 0.01296296,
1484                        0.03148148, 0.01296296, 0.03148148,
1485                        0.01296296, 0.03148148, 0.01296296,
1486                        0.03148148, 0.01296296, 0.03148148,
1487                        -0.04259259, -0.02407407, -0.04259259,
1488                        -0.02407407, -0.04259259, -0.02407407,
1489                        -0.04259259, -0.02407407, -0.04259259,
1490                        -0.02407407, -0.04259259, -0.02407407,
1491                        -0.09814815, -0.07962963, -0.09814815,
1492                        -0.07962963, -0.09814815, -0.07962963,
1493                        -0.09814815, -0.07962963, -0.09814815,
1494                        -0.07962963, -0.09814815, -0.07962963,
1495                        -0.1537037 , -0.13518519, -0.1537037,
1496                        -0.13518519, -0.1537037, -0.13518519,
1497                        -0.1537037 , -0.13518519, -0.1537037,
1498                        -0.13518519, -0.1537037, -0.13518519,
1499                        -0.20925926, -0.19074074, -0.20925926,
1500                        -0.19074074, -0.20925926, -0.19074074,
1501                        -0.20925926, -0.19074074, -0.20925926,
1502                        -0.19074074, -0.20925926, -0.19074074,
1503                        -0.26481481, -0.2462963, -0.26481481,
1504                        -0.2462963, -0.26481481, -0.2462963,
1505                        -0.26481481, -0.2462963, -0.26481481,
1506                        -0.2462963, -0.26481481, -0.2462963])
1507
1508
1509        #print domain.quantities['level'].extrapolate_second_order()
1510        #domain.distribute_to_vertices_and_edges()
1511        #print domain.quantities['level'].vertex_values[:,0]       
1512       
1513        #Evolution
1514        for t in domain.evolve(yieldstep = 0.05, finaltime = 0.1):
1515            pass
1516
1517       
1518        #Data from earlier version of pyvolution ft=0.1
1519        assert allclose(domain.min_timestep, 0.0376895634803) 
1520        assert allclose(domain.max_timestep, 0.0415635655309)
1521
1522
1523        assert allclose(domain.quantities['level'].centroid_values,
1524                        [0.00855788, 0.01575204, 0.00994606, 0.01717072,
1525                         0.01005985, 0.01716362, 0.01005985, 0.01716299,
1526                         0.01007098, 0.01736248, 0.01216452, 0.02026776,
1527                         -0.04462374, -0.02479045, -0.04199789, -0.0229465,
1528                         -0.04184033, -0.02295693, -0.04184013, -0.02295675,
1529                         -0.04184486, -0.0228168, -0.04028876, -0.02036486,
1530                         -0.10029444, -0.08170809, -0.09772846, -0.08021704,
1531                         -0.09760006, -0.08022143, -0.09759984, -0.08022124,
1532                         -0.09760261, -0.08008893, -0.09603914, -0.07758209,
1533                         -0.15584152, -0.13723138, -0.15327266, -0.13572906,
1534                         -0.15314427, -0.13573349, -0.15314405, -0.13573331,
1535                         -0.15314679, -0.13560104, -0.15158523, -0.13310701,
1536                         -0.21208605, -0.19283913, -0.20955631, -0.19134189,
1537                         -0.20942821, -0.19134598, -0.20942799, -0.1913458,
1538                         -0.20943005, -0.19120952, -0.20781177, -0.18869401,
1539                         -0.25384082, -0.2463294, -0.25047649, -0.24464654,
1540                         -0.25031159, -0.24464253, -0.25031112, -0.24464253,
1541                         -0.25031463, -0.24454764, -0.24885323, -0.24286438])
1542
1543       
1544       
1545               
1546    def test_bedslope_problem_second_order_two_yieldsteps(self):
1547
1548        from mesh_factory import rectangular
1549        from shallow_water import Domain, Reflective_boundary, Constant_height
1550        from Numeric import array
1551
1552        #Create basic mesh
1553        points, vertices, boundary = rectangular(6, 6)
1554
1555        #Create shallow water domain
1556        domain = Domain(points, vertices, boundary)
1557        domain.smooth = False
1558        domain.default_order=2
1559
1560        #Bed-slope and friction at vertices (and interpolated elsewhere)
1561        def x_slope(x, y):
1562            return -x/3
1563
1564        domain.set_quantity('elevation', x_slope)
1565
1566        # Boundary conditions
1567        Br = Reflective_boundary(domain)
1568        domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
1569
1570        #Initial condition
1571        domain.set_quantity('level', Constant_height(x_slope, 0.05))
1572        domain.check_integrity()
1573
1574        assert allclose(domain.quantities['level'].centroid_values,
1575                        [0.01296296, 0.03148148, 0.01296296,
1576                        0.03148148, 0.01296296, 0.03148148,
1577                        0.01296296, 0.03148148, 0.01296296,
1578                        0.03148148, 0.01296296, 0.03148148,
1579                        -0.04259259, -0.02407407, -0.04259259,
1580                        -0.02407407, -0.04259259, -0.02407407,
1581                        -0.04259259, -0.02407407, -0.04259259,
1582                        -0.02407407, -0.04259259, -0.02407407,
1583                        -0.09814815, -0.07962963, -0.09814815,
1584                        -0.07962963, -0.09814815, -0.07962963,
1585                        -0.09814815, -0.07962963, -0.09814815,
1586                        -0.07962963, -0.09814815, -0.07962963,
1587                        -0.1537037 , -0.13518519, -0.1537037,
1588                        -0.13518519, -0.1537037, -0.13518519,
1589                        -0.1537037 , -0.13518519, -0.1537037,
1590                        -0.13518519, -0.1537037, -0.13518519,
1591                        -0.20925926, -0.19074074, -0.20925926,
1592                        -0.19074074, -0.20925926, -0.19074074,
1593                        -0.20925926, -0.19074074, -0.20925926,
1594                        -0.19074074, -0.20925926, -0.19074074,
1595                        -0.26481481, -0.2462963, -0.26481481,
1596                        -0.2462963, -0.26481481, -0.2462963,
1597                        -0.26481481, -0.2462963, -0.26481481,
1598                        -0.2462963, -0.26481481, -0.2462963])
1599
1600
1601        #print domain.quantities['level'].extrapolate_second_order()
1602        #domain.distribute_to_vertices_and_edges()
1603        #print domain.quantities['level'].vertex_values[:,0]       
1604       
1605        #Evolution
1606        for t in domain.evolve(yieldstep = 0.05, finaltime = 0.1):   #0.05??
1607            #domain.write_time()           
1608            pass
1609
1610
1611       
1612        assert allclose(domain.quantities['level'].centroid_values, 
1613                 [0.00855788, 0.01575204, 0.00994606, 0.01717072, 0.01005985,
1614                  0.01716362, 0.01005985, 0.01716299, 0.01007098, 0.01736248,
1615                  0.01216452, 0.02026776, -0.04462374, -0.02479045, -0.04199789,
1616                  -0.0229465, -0.04184033, -0.02295693, -0.04184013,
1617                  -0.02295675, -0.04184486, -0.0228168, -0.04028876,
1618                  -0.02036486, -0.10029444, -0.08170809, -0.09772846,
1619                  -0.08021704, -0.09760006, -0.08022143, -0.09759984,
1620                  -0.08022124, -0.09760261, -0.08008893, -0.09603914,
1621                  -0.07758209, -0.15584152, -0.13723138, -0.15327266,
1622                  -0.13572906, -0.15314427, -0.13573349, -0.15314405,
1623                  -0.13573331, -0.15314679, -0.13560104, -0.15158523,
1624                  -0.13310701, -0.21208605, -0.19283913, -0.20955631,
1625                  -0.19134189, -0.20942821, -0.19134598, -0.20942799,
1626                  -0.1913458, -0.20943005, -0.19120952, -0.20781177,
1627                  -0.18869401, -0.25384082, -0.2463294, -0.25047649,
1628                  -0.24464654, -0.25031159, -0.24464253, -0.25031112,
1629                  -0.24464253, -0.25031463, -0.24454764, -0.24885323,
1630                  -0.24286438])
1631
1632
1633
1634    def test_bedslope_problem_second_order_more_steps(self):
1635
1636        from mesh_factory import rectangular
1637        from shallow_water import Domain, Reflective_boundary, Constant_height
1638        from Numeric import array
1639
1640        #Create basic mesh
1641        points, vertices, boundary = rectangular(6, 6)
1642
1643        #Create shallow water domain
1644        domain = Domain(points, vertices, boundary)
1645        domain.smooth = False
1646        domain.default_order=2
1647
1648        #Bed-slope and friction at vertices (and interpolated elsewhere)
1649        def x_slope(x, y):
1650            return -x/3
1651
1652        domain.set_quantity('elevation', x_slope)
1653
1654        # Boundary conditions
1655        Br = Reflective_boundary(domain)
1656        domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
1657
1658        #Initial condition
1659        domain.set_quantity('level', Constant_height(x_slope, 0.05))
1660        domain.check_integrity()
1661
1662        assert allclose(domain.quantities['level'].centroid_values,
1663                        [0.01296296, 0.03148148, 0.01296296,
1664                        0.03148148, 0.01296296, 0.03148148,
1665                        0.01296296, 0.03148148, 0.01296296,
1666                        0.03148148, 0.01296296, 0.03148148,
1667                        -0.04259259, -0.02407407, -0.04259259,
1668                        -0.02407407, -0.04259259, -0.02407407,
1669                        -0.04259259, -0.02407407, -0.04259259,
1670                        -0.02407407, -0.04259259, -0.02407407,
1671                        -0.09814815, -0.07962963, -0.09814815,
1672                        -0.07962963, -0.09814815, -0.07962963,
1673                        -0.09814815, -0.07962963, -0.09814815,
1674                        -0.07962963, -0.09814815, -0.07962963,
1675                        -0.1537037 , -0.13518519, -0.1537037,
1676                        -0.13518519, -0.1537037, -0.13518519,
1677                        -0.1537037 , -0.13518519, -0.1537037,
1678                        -0.13518519, -0.1537037, -0.13518519,
1679                        -0.20925926, -0.19074074, -0.20925926,
1680                        -0.19074074, -0.20925926, -0.19074074,
1681                        -0.20925926, -0.19074074, -0.20925926,
1682                        -0.19074074, -0.20925926, -0.19074074,
1683                        -0.26481481, -0.2462963, -0.26481481,
1684                        -0.2462963, -0.26481481, -0.2462963,
1685                        -0.26481481, -0.2462963, -0.26481481,
1686                        -0.2462963, -0.26481481, -0.2462963])
1687
1688
1689        #print domain.quantities['level'].extrapolate_second_order()
1690        #domain.distribute_to_vertices_and_edges()
1691        #print domain.quantities['level'].vertex_values[:,0]       
1692       
1693        #Evolution
1694        for t in domain.evolve(yieldstep = 0.05, finaltime = 0.5):
1695            pass
1696
1697       
1698        assert allclose(domain.quantities['level'].centroid_values,
1699     [-0.02907028, -0.01475478, -0.02973417, -0.01447186, -0.02932665, -0.01428285,
1700      -0.02901975, -0.0141361,  -0.02898816, -0.01418135, -0.02961409, -0.01403487,
1701      -0.07597998, -0.06252591, -0.07664854, -0.06312532, -0.07638287, -0.06265139,
1702      -0.07571145, -0.06235231, -0.0756817,  -0.06245309, -0.07652292, -0.06289946,
1703      -0.12367464, -0.11088981, -0.12237277, -0.11115338, -0.1218934,  -0.1107174, 
1704      -0.12081485, -0.11000491, -0.12038451, -0.11010335, -0.12102113, -0.11012105,
1705      -0.16909116, -0.15831543, -0.16730214, -0.15786249, -0.1665493,  -0.15697919,
1706      -0.16496618, -0.15559852, -0.16338679, -0.15509088, -0.16364092, -0.15424423,
1707      -0.18771107, -0.19903904, -0.18903759, -0.19858437, -0.18701552, -0.19697797,
1708      -0.1833593,  -0.19505871, -0.1818806,  -0.19418042, -0.18586159, -0.19576946,
1709      -0.13986873, -0.14170053, -0.14132188, -0.14560674, -0.14095617, -0.14373292,
1710      -0.13785933, -0.14033364, -0.13592955, -0.13936356, -0.13596008, -0.14216296])
1711
1712        assert allclose(domain.quantities['xmomentum'].centroid_values,
1713     [ 0.00831121,  0.00317948,  0.00731797,  0.00334939,  0.00764717,  0.00348053,
1714       0.00788729,  0.00356522,  0.00780649,  0.00341919,  0.00693525,  0.00310375,
1715       0.02166196,  0.01421475,  0.02017737,  0.01316839,  0.02037015,  0.01368659,
1716       0.02106,     0.01399161,  0.02074514,  0.01354935,  0.01887407,  0.0123113, 
1717       0.03775083,  0.02855197,  0.03689337,  0.02759782,  0.03732848,  0.02812072,
1718       0.03872545,  0.02913348,  0.03880939,  0.02803804,  0.03546499,  0.0260039, 
1719       0.0632131,   0.04730634,  0.0576324,   0.04592336,  0.05790921,  0.04690514,
1720       0.05986467,  0.04871165,  0.06170068,  0.04811572,  0.05657041,  0.04416292,
1721       0.08489642,  0.07188097,  0.07835261,  0.06843406,  0.07986412,  0.0698247, 
1722       0.08201071,  0.07216756,  0.08378418,  0.07273624,  0.080399,    0.06645841,
1723       0.01631548,  0.04691608,  0.0206632,   0.044409,    0.02115518,  0.04560305,
1724       0.02160608,  0.04663725,  0.02174734,  0.04795559,  0.02281427,  0.05667111])
1725
1726       
1727        assert allclose(domain.quantities['ymomentum'].centroid_values, 
1728     [ 1.45876601e-004, -3.24627393e-004, -1.57572719e-004, -2.92790187e-004,
1729      -9.90988382e-005, -3.06677335e-004, -1.62493106e-004, -3.71310004e-004,
1730      -1.99445058e-004, -3.28493467e-004,  6.68217349e-005, -8.42042805e-006,
1731       5.05093371e-004, -1.42842214e-004, -6.81454718e-005, -5.02084057e-004,
1732      -8.50583861e-005, -4.65443981e-004, -1.96406564e-004, -5.88889562e-004,
1733      -2.70160173e-004, -5.35485454e-004,  2.60780997e-004,  3.12145471e-005,
1734       5.16189608e-004,  1.07069062e-004,  9.29989252e-005, -3.71211119e-004,
1735       1.16350246e-004, -3.82407830e-004, -1.62077969e-004, -6.30906636e-004,
1736      -4.74025708e-004, -6.94463009e-004,  6.15092843e-005,  2.22106820e-004,
1737      -6.29589294e-004,  2.43611937e-004, -5.88125094e-004, -6.94293192e-005,
1738      -4.17914641e-004,  6.64609019e-005, -7.68334577e-004, -3.40232101e-004,
1739      -1.67424308e-003, -7.39485066e-004, -1.59966988e-003,  5.68262838e-005,
1740      -1.48470633e-003, -1.84554882e-003, -2.27200099e-003, -1.67506848e-003,
1741      -1.95610258e-003, -1.47638801e-003, -1.73779477e-003, -1.85498791e-003,
1742      -2.01357843e-003, -2.17675471e-003, -1.65783870e-003, -1.15818681e-003,
1743      -1.18663036e-003, -2.94229849e-003, -3.59309018e-003, -5.13496584e-003,
1744      -6.17359400e-003, -5.98761937e-003, -6.00540116e-003, -5.01121966e-003,
1745      -4.50964850e-003, -3.06319963e-003,  6.08950810e-004, -4.79537921e-004])
1746       
1747
1748       
1749
1750    def test_temp_play(self):
1751
1752        from mesh_factory import rectangular
1753        from shallow_water import Domain, Reflective_boundary, Constant_height
1754        from Numeric import array
1755
1756        #Create basic mesh
1757        points, vertices, boundary = rectangular(5, 5)
1758
1759        #Create shallow water domain
1760        domain = Domain(points, vertices, boundary)
1761        domain.smooth = False
1762        domain.default_order=2
1763
1764        #Bed-slope and friction at vertices (and interpolated elsewhere)
1765        def x_slope(x, y):
1766            return -x/3
1767
1768        domain.set_quantity('elevation', x_slope)
1769
1770        # Boundary conditions
1771        Br = Reflective_boundary(domain)
1772        domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
1773
1774        #Initial condition
1775        domain.set_quantity('level', Constant_height(x_slope, 0.05))
1776        domain.check_integrity()
1777
1778        #Evolution
1779        for t in domain.evolve(yieldstep = 0.05, finaltime = 0.1):
1780            pass
1781
1782        assert allclose(domain.quantities['level'].centroid_values[:4],
1783                        [0.00206836, 0.01296714, 0.00363415, 0.01438924])
1784        assert allclose(domain.quantities['xmomentum'].centroid_values[:4],                     
1785                        [0.01360154, 0.00671133, 0.01264578, 0.00648503])
1786        assert allclose(domain.quantities['ymomentum'].centroid_values[:4],
1787                        [-1.19201077e-003, -7.23647546e-004, 
1788                        -6.39083123e-005, 6.29815168e-005])
1789       
1790       
1791    def test_complex_bed(self):
1792        #No friction is tested here
1793       
1794        from shallow_water import Domain, Reflective_boundary, Dirichlet_boundary,\
1795             Transmissive_boundary, Time_boundary,\
1796             Weir_simple as Weir, Constant_height
1797
1798        from mesh_factory import rectangular
1799        from Numeric import array
1800   
1801        N = 12
1802        points, vertices, boundary = rectangular(N, N/2, len1=1.2,len2=0.6,
1803                                                 origin=(-0.07, 0))
1804
1805
1806        domain = Domain(points, vertices, boundary)
1807        domain.smooth = False
1808        domain.visualise = False
1809        domain.default_order=2
1810
1811       
1812        inflow_stage = 0.1
1813        Z = Weir(inflow_stage)
1814        domain.set_quantity('elevation', Z)
1815
1816        Br = Reflective_boundary(domain)
1817        Bd = Dirichlet_boundary([inflow_stage, 0.0, 0.0])
1818        domain.set_boundary({'left': Bd, 'right': Br, 'bottom': Br, 'top': Br})
1819
1820        domain.set_quantity('level', Constant_height(Z, 0.))
1821
1822        for t in domain.evolve(yieldstep = 0.02, finaltime = 0.2):
1823            pass
1824
1825
1826        #print domain.quantities['level'].centroid_values
1827
1828        #FIXME: These numbers were from version before 25/10       
1829        #assert allclose(domain.quantities['level'].centroid_values,
1830# [3.95822638e-002,  5.61022588e-002,  4.66437868e-002,  5.73081011e-002,
1831#  4.72394613e-002,  5.74684939e-002,  4.74309483e-002,  5.77458084e-002,
1832#  4.80628177e-002,  5.85656225e-002,  4.90498542e-002,  6.02609831e-002,
1833#  1.18470315e-002,  1.75136443e-002,  1.18035266e-002,  2.15565695e-002,
1834#  1.31620268e-002,  2.14351640e-002,  1.32351076e-002,  2.15450687e-002,
1835#  1.36414028e-002,  2.24274619e-002,  1.51689511e-002,  2.21789655e-002,
1836# -7.54337535e-003, -6.86362021e-004, -7.74146760e-003, -1.83756530e-003,
1837# -8.16773628e-003, -4.49916813e-004, -8.08202599e-003, -3.91118720e-004,
1838# -8.10292716e-003, -3.88584984e-004, -7.35226124e-003,  2.73985295e-004,
1839#  1.86166683e-001,  8.74070369e-002,  1.86166712e-001,  8.74035875e-002,
1840#  6.11666935e-002, -3.76173225e-002, -6.38333276e-002, -3.76147365e-002,
1841#  6.11666725e-002,  8.73846774e-002,  1.86166697e-001,  8.74171550e-002,
1842# -4.83333333e-002,  1.18333333e-001, -4.83333333e-002,  1.18333333e-001,
1843# -4.83333333e-002, -6.66666667e-003, -1.73333333e-001, -1.31666667e-001,
1844# -1.73333333e-001, -6.66666667e-003, -4.83333333e-002,  1.18333333e-001,
1845# -2.48333333e-001, -2.31666667e-001, -2.48333333e-001, -2.31666667e-001,
1846# -2.48333333e-001, -2.31666667e-001, -2.48333333e-001, -2.31666667e-001,
1847# -2.48333333e-001, -2.31666667e-001, -2.48333333e-001, -2.31666667e-001,
1848# -4.65000000e-001, -3.65000000e-001, -4.65000000e-001, -3.65000000e-001,
1849# -4.65000000e-001, -3.65000000e-001, -4.65000000e-001, -3.65000000e-001,
1850# -4.65000000e-001, -3.65000000e-001, -4.65000000e-001, -3.65000000e-001,
1851# -5.98333333e-001, -5.81666667e-001, -5.98333333e-001, -5.81666667e-001,
1852# -5.98333333e-001, -5.81666667e-001, -5.98333333e-001, -5.81666667e-001,
1853# -5.98333333e-001, -5.81666667e-001, -5.98333333e-001, -5.81666667e-001,
1854# -6.48333333e-001, -6.31666667e-001, -6.48333333e-001, -6.31666667e-001,
1855# -6.48333333e-001, -6.31666667e-001, -6.48333333e-001, -6.31666667e-001,
1856# -6.48333333e-001, -6.31666667e-001, -6.48333333e-001, -6.31666667e-001,
1857# -5.31666667e-001, -5.98333333e-001, -5.31666667e-001, -5.98333333e-001,
1858# -5.31666667e-001, -5.98333333e-001, -5.31666667e-001, -5.98333333e-001,
1859# -5.31666667e-001, -5.98333333e-001, -5.31666667e-001, -5.98333333e-001,
1860# -4.98333333e-001, -4.81666667e-001, -4.98333333e-001, -4.81666667e-001,
1861# -4.98333333e-001, -4.81666667e-001, -4.98333333e-001, -4.81666667e-001,
1862# -4.98333333e-001, -4.81666667e-001, -4.98333333e-001, -4.81666667e-001,
1863# -5.48333333e-001, -5.31666667e-001, -5.48333333e-001, -5.31666667e-001,
1864# -5.48333333e-001, -5.31666667e-001, -5.48333333e-001, -5.31666667e-001,
1865# -5.48333333e-001, -5.31666667e-001, -5.48333333e-001, -5.31666667e-001])
1866
1867
1868       
1869   
1870#-------------------------------------------------------------
1871if __name__ == "__main__":
1872    suite = unittest.makeSuite(TestCase,'test')
1873    runner = unittest.TextTestRunner()
1874    runner.run(suite)
1875
Note: See TracBrowser for help on using the repository browser.