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

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

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