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

Last change on this file since 668 was 648, checked in by ole, 20 years ago

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