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

Last change on this file since 1013 was 1013, checked in by ole, 19 years ago
File size: 103.8 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, take
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        #Check error check
97        try:
98            rotate(r, array([1,1,1]) )
99        except:
100            pass
101        else:
102            raise 'Should have raised an exception'
103
104    def test_flux_zero_case(self):
105        ql = zeros( 3, Float )
106        qr = zeros( 3, Float )
107        normal = zeros( 2, Float )
108        zl = zr = 0.
109        flux, max_speed = flux_function(normal, ql, qr, zl, zr)
110
111        assert allclose(flux, [0,0,0])
112        assert max_speed == 0.
113
114    def test_flux_constants(self):
115        w = 2.0
116       
117        normal = array([1.,0])       
118        ql = array([w, 0, 0])
119        qr = array([w, 0, 0])
120        zl = zr = 0.
121        h = w - (zl+zr)/2
122       
123        flux, max_speed = flux_function(normal, ql, qr, zl, zr)
124
125        assert allclose(flux, [0., 0.5*g*h**2, 0.])
126        assert max_speed == sqrt(g*h)
127       
128    #def test_flux_slope(self):
129    #    #FIXME: TODO
130    #    w = 2.0
131    #   
132    #    normal = array([1.,0])       
133    #    ql = array([w, 0, 0])
134    #    qr = array([w, 0, 0])
135    #    zl = zr = 0.
136    #    h = w - (zl+zr)/2
137    #   
138    #    flux, max_speed = flux_function(normal, ql, qr, zl, zr)
139    #
140    #    assert allclose(flux, [0., 0.5*g*h**2, 0.])
141    #    assert max_speed == sqrt(g*h)
142       
143
144    def test_flux1(self):
145        #Use data from previous version of pyvolution
146        normal = array([1.,0])       
147        ql = array([-0.2, 2, 3])
148        qr = array([-0.2, 2, 3])
149        zl = zr = -0.5
150        flux, max_speed = flux_function(normal, ql, qr, zl, zr)
151
152        assert allclose(flux, [2.,13.77433333, 20.])
153        assert allclose(max_speed, 8.38130948661)
154
155
156    def test_flux2(self):
157        #Use data from previous version of pyvolution
158        normal = array([0., -1.])       
159        ql = array([-0.075, 2, 3])
160        qr = array([-0.075, 2, 3])
161        zl = zr = -0.375
162        flux, max_speed = flux_function(normal, ql, qr, zl, zr)
163
164        assert allclose(flux, [-3.,-20.0, -30.441])
165        assert allclose(max_speed, 11.7146428199)
166
167    def test_flux3(self):
168        #Use data from previous version of pyvolution       
169        normal = array([-sqrt(2)/2, sqrt(2)/2])       
170        ql = array([-0.075, 2, 3])
171        qr = array([-0.075, 2, 3])
172        zl = zr = -0.375
173        flux, max_speed = flux_function(normal, ql, qr, zl, zr)
174
175        assert allclose(flux, [sqrt(2)/2, 4.40221112, 7.3829019])
176        assert allclose(max_speed, 4.0716654239)
177
178    def test_flux4(self):
179        #Use data from previous version of pyvolution               
180        normal = array([-sqrt(2)/2, sqrt(2)/2])       
181        ql = array([-0.34319278, 0.10254161, 0.07273855])
182        qr = array([-0.30683287, 0.1071986, 0.05930515])
183        zl = zr = -0.375
184        flux, max_speed = flux_function(normal, ql, qr, zl, zr)
185
186        assert allclose(flux, [-0.04072676, -0.07096636, -0.01604364])
187        assert allclose(max_speed, 1.31414103233)                               
188
189    def test_sw_domain_simple(self):
190        a = [0.0, 0.0]
191        b = [0.0, 2.0]
192        c = [2.0,0.0]
193        d = [0.0, 4.0]
194        e = [2.0, 2.0]
195        f = [4.0,0.0]
196
197        points = [a, b, c, d, e, f]
198        #bac, bce, ecf, dbe, daf, dae
199        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4], [3,0,5], [3,0,4]]
200       
201        domain = Domain(points, vertices)       
202        domain.check_integrity()
203
204        for name in ['stage', 'xmomentum', 'ymomentum',
205                     'elevation', 'friction']:
206            assert domain.quantities.has_key(name)
207
208
209        assert domain.get_conserved_quantities(0, edge=1) == 0.   
210
211
212    def test_boundary_conditions(self):
213       
214        a = [0.0, 0.0]
215        b = [0.0, 2.0]
216        c = [2.0,0.0]
217        d = [0.0, 4.0]
218        e = [2.0, 2.0]
219        f = [4.0,0.0]
220
221        points = [a, b, c, d, e, f]
222        #bac, bce, ecf, dbe
223        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4] ]
224        boundary = { (0, 0): 'Third',
225                     (0, 2): 'First',
226                     (2, 0): 'Second',
227                     (2, 1): 'Second',
228                     (3, 1): 'Second',
229                     (3, 2): 'Third'}                                         
230                     
231
232        domain = Domain(points, vertices, boundary)
233        domain.check_integrity()
234
235
236        domain.set_quantity('stage', [[1,2,3], [5,5,5],
237                                      [0,0,9], [-6, 3, 3]])
238
239        domain.set_quantity('xmomentum', [[1,1,1], [2,2,2],
240                                          [3,3,3], [4, 4, 4]])
241
242        domain.set_quantity('ymomentum', [[10,10,10], [20,20,20],
243                                          [30,30,30], [40, 40, 40]])       
244
245
246        D = Dirichlet_boundary([5,2,1])
247        T = Transmissive_boundary(domain)
248        R = Reflective_boundary(domain)
249        domain.set_boundary( {'First': D, 'Second': T, 'Third': R})
250       
251        domain.update_boundary()
252
253        #Stage
254        assert domain.quantities['stage'].boundary_values[0] == 2.5
255        assert domain.quantities['stage'].boundary_values[0] ==\
256               domain.get_conserved_quantities(0, edge=0)[0] #Reflective (2.5)
257        assert domain.quantities['stage'].boundary_values[1] == 5. #Dirichlet
258        assert domain.quantities['stage'].boundary_values[2] ==\
259               domain.get_conserved_quantities(2, edge=0)[0] #Transmissive (4.5)
260        assert domain.quantities['stage'].boundary_values[3] ==\
261               domain.get_conserved_quantities(2, edge=1)[0] #Transmissive (4.5)
262        assert domain.quantities['stage'].boundary_values[4] ==\
263               domain.get_conserved_quantities(3, edge=1)[0] #Transmissive (-1.5)
264        assert domain.quantities['stage'].boundary_values[5] ==\
265               domain.get_conserved_quantities(3, edge=2)[0] #Reflective (-1.5)
266
267        #Xmomentum
268        assert domain.quantities['xmomentum'].boundary_values[0] == 1.0 #Reflective
269        assert domain.quantities['xmomentum'].boundary_values[1] == 2. #Dirichlet
270        assert domain.quantities['xmomentum'].boundary_values[2] ==\
271               domain.get_conserved_quantities(2, edge=0)[1] #Transmissive
272        assert domain.quantities['xmomentum'].boundary_values[3] ==\
273               domain.get_conserved_quantities(2, edge=1)[1] #Transmissive
274        assert domain.quantities['xmomentum'].boundary_values[4] ==\
275               domain.get_conserved_quantities(3, edge=1)[1] #Transmissive
276        assert domain.quantities['xmomentum'].boundary_values[5] == -4.0  #Reflective
277       
278        #Ymomentum
279        assert domain.quantities['ymomentum'].boundary_values[0] == -10.0 #Reflective
280        assert domain.quantities['ymomentum'].boundary_values[1] == 1.  #Dirichlet
281        assert domain.quantities['ymomentum'].boundary_values[2] == 30. #Transmissive
282        assert domain.quantities['ymomentum'].boundary_values[3] == 30. #Transmissive
283        assert domain.quantities['ymomentum'].boundary_values[4] == 40. #Transmissive
284        assert domain.quantities['ymomentum'].boundary_values[5] == 40. #Reflective
285
286
287    def test_boundary_conditionsII(self):
288       
289        a = [0.0, 0.0]
290        b = [0.0, 2.0]
291        c = [2.0,0.0]
292        d = [0.0, 4.0]
293        e = [2.0, 2.0]
294        f = [4.0,0.0]
295
296        points = [a, b, c, d, e, f]
297        #bac, bce, ecf, dbe
298        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4] ]
299        boundary = { (0, 0): 'Third',
300                     (0, 2): 'First',
301                     (2, 0): 'Second',
302                     (2, 1): 'Second',
303                     (3, 1): 'Second',
304                     (3, 2): 'Third',
305                     (0, 1): 'Internal'}                                         
306                     
307
308        domain = Domain(points, vertices, boundary)
309        domain.check_integrity()
310
311
312        domain.set_quantity('stage', [[1,2,3], [5,5,5],
313                                      [0,0,9], [-6, 3, 3]])
314
315        domain.set_quantity('xmomentum', [[1,1,1], [2,2,2],
316                                          [3,3,3], [4, 4, 4]])
317
318        domain.set_quantity('ymomentum', [[10,10,10], [20,20,20],
319                                          [30,30,30], [40, 40, 40]])       
320
321
322        D = Dirichlet_boundary([5,2,1])
323        T = Transmissive_boundary(domain)
324        R = Reflective_boundary(domain)
325        domain.set_boundary( {'First': D, 'Second': T,
326                              'Third': R, 'Internal': None})
327
328        domain.update_boundary()
329        domain.check_integrity()
330
331
332    def test_compute_fluxes0(self):
333        #Do a full triangle and check that fluxes cancel out for
334        #the constant stage case
335       
336        a = [0.0, 0.0]
337        b = [0.0, 2.0]
338        c = [2.0,0.0]
339        d = [0.0, 4.0]
340        e = [2.0, 2.0]
341        f = [4.0,0.0]
342
343        points = [a, b, c, d, e, f]
344        #bac, bce, ecf, dbe
345        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
346       
347        domain = Domain(points, vertices)       
348        domain.set_quantity('stage', [[2,2,2], [2,2,2],
349                                      [2,2,2], [2,2,2]])
350        domain.check_integrity()
351
352        assert allclose(domain.neighbours, [[-1,1,-1], [2,3,0], [-1,-1,1],[1,-1,-1]])
353        assert allclose(domain.neighbour_edges, [[-1,2,-1], [2,0,1], [-1,-1,0],[1,-1,-1]])       
354
355        zl=zr=0. #Assume flat bed
356       
357        #Flux across right edge of volume 1
358        normal = domain.get_normal(1,0)       
359        ql = domain.get_conserved_quantities(vol_id=1, edge=0)
360        qr = domain.get_conserved_quantities(vol_id=2, edge=2)
361        flux0, 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(2,2)
366        flux, max_speed = flux_function(normal, ql, qr, zl, zr)       
367        assert allclose(flux + flux0, 0.)
368   
369        #Flux across upper edge of volume 1
370        normal = domain.get_normal(1,1)       
371        ql = domain.get_conserved_quantities(vol_id=1, edge=1)
372        qr = domain.get_conserved_quantities(vol_id=3, edge=0)
373        flux1, max_speed = flux_function(normal, ql, qr, zl, zr)
374       
375        #Check that flux seen from other triangles is inverse
376        tmp = qr; qr=ql; ql=tmp
377        normal = domain.get_normal(3,0)
378        flux, max_speed = flux_function(normal, ql, qr, zl, zr)       
379        assert allclose(flux + flux1, 0.)
380       
381        #Flux across lower left hypotenuse of volume 1
382        normal = domain.get_normal(1,2)       
383        ql = domain.get_conserved_quantities(vol_id=1, edge=2)
384        qr = domain.get_conserved_quantities(vol_id=0, edge=1)
385        flux2, max_speed = flux_function(normal, ql, qr, zl, zr)
386       
387        #Check that flux seen from other triangles is inverse
388        tmp = qr; qr=ql; ql=tmp
389        normal = domain.get_normal(0,1)
390        flux, max_speed = flux_function(normal, ql, qr, zl, zr)       
391        assert allclose(flux + flux2, 0.)
392
393
394        #Scale by edgelengths, add up anc check that total flux is zero
395        e0 = domain.edgelengths[1, 0]
396        e1 = domain.edgelengths[1, 1]
397        e2 = domain.edgelengths[1, 2]
398       
399        assert allclose(e0*flux0+e1*flux1+e2*flux2, 0.)
400
401        #Now check that compute_flux yields zeros as well
402        domain.compute_fluxes()
403
404        for name in ['stage', 'xmomentum', 'ymomentum']:
405            #print name, domain.quantities[name].explicit_update
406            assert allclose(domain.quantities[name].explicit_update[1], 0)
407
408
409
410    def test_compute_fluxes1(self):
411        #Use values from previous version
412       
413        a = [0.0, 0.0]
414        b = [0.0, 2.0]
415        c = [2.0,0.0]
416        d = [0.0, 4.0]
417        e = [2.0, 2.0]
418        f = [4.0,0.0]
419
420        points = [a, b, c, d, e, f]
421        #bac, bce, ecf, dbe
422        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
423       
424        domain = Domain(points, vertices)
425        val0 = 2.+2.0/3
426        val1 = 4.+4.0/3
427        val2 = 8.+2.0/3
428        val3 = 2.+8.0/3
429       
430        domain.set_quantity('stage', [[val0, val0, val0], [val1, val1, val1],
431                                      [val2, val2, val2], [val3, val3, val3]])
432        domain.check_integrity()
433
434        zl=zr=0. #Assume flat bed
435
436        #Flux across right edge of volume 1
437        normal = domain.get_normal(1,0)       
438        ql = domain.get_conserved_quantities(vol_id=1, edge=0)
439        qr = domain.get_conserved_quantities(vol_id=2, edge=2)
440        flux0, max_speed = flux_function(normal, ql, qr, zl, zr)
441        assert allclose(flux0, [-15.3598804, 253.71111111, 0.])
442        assert allclose(max_speed, 9.21592824046)
443
444
445        #Flux across upper edge of volume 1
446        normal = domain.get_normal(1,1)       
447        ql = domain.get_conserved_quantities(vol_id=1, edge=1)
448        qr = domain.get_conserved_quantities(vol_id=3, edge=0)
449        flux1, max_speed = flux_function(normal, ql, qr, zl, zr)
450        assert allclose(flux1, [2.4098563, 0., 123.04444444])
451        assert allclose(max_speed, 7.22956891292)
452
453        #Flux across lower left hypotenuse of volume 1
454        normal = domain.get_normal(1,2)       
455        ql = domain.get_conserved_quantities(vol_id=1, edge=2)
456        qr = domain.get_conserved_quantities(vol_id=0, edge=1)
457        flux2, max_speed = flux_function(normal, ql, qr, zl, zr)
458
459        assert allclose(flux2, [9.63942522, -61.59685738, -61.59685738])
460        assert allclose(max_speed, 7.22956891292)
461       
462        #Scale, add up and check that compute_fluxes is correct for vol 1
463        e0 = domain.edgelengths[1, 0]
464        e1 = domain.edgelengths[1, 1]
465        e2 = domain.edgelengths[1, 2]
466
467        total_flux = -(e0*flux0+e1*flux1+e2*flux2)/domain.areas[1]
468        assert allclose(total_flux, [-0.68218178, -166.6, -35.93333333])
469
470
471        domain.compute_fluxes()
472
473        #assert allclose(total_flux, domain.explicit_update[1,:])
474        for i, name in enumerate(['stage', 'xmomentum', 'ymomentum']):
475            assert allclose(total_flux[i],
476                            domain.quantities[name].explicit_update[1])
477
478        #assert allclose(domain.explicit_update, [
479        #    [0., -69.68888889, -69.68888889],
480        #    [-0.68218178, -166.6, -35.93333333],
481        #    [-111.77316251, 69.68888889, 0.],
482        #    [-35.68522449, 0., 69.68888889]])
483
484        assert allclose(domain.quantities['stage'].explicit_update,
485                        [0., -0.68218178, -111.77316251, -35.68522449])
486        assert allclose(domain.quantities['xmomentum'].explicit_update,
487                        [-69.68888889, -166.6, 69.68888889, 0])
488        assert allclose(domain.quantities['ymomentum'].explicit_update,
489                        [-69.68888889, -35.93333333, 0., 69.68888889])       
490
491       
492        #assert allclose(domain.quantities[name].explicit_update
493
494       
495       
496
497
498    def test_compute_fluxes2(self):
499        #Random values, incl momentum
500       
501        a = [0.0, 0.0]
502        b = [0.0, 2.0]
503        c = [2.0,0.0]
504        d = [0.0, 4.0]
505        e = [2.0, 2.0]
506        f = [4.0,0.0]
507
508        points = [a, b, c, d, e, f]
509        #bac, bce, ecf, dbe
510        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
511       
512        domain = Domain(points, vertices)
513        val0 = 2.+2.0/3
514        val1 = 4.+4.0/3
515        val2 = 8.+2.0/3
516        val3 = 2.+8.0/3
517
518        zl=zr=0 #Assume flat zero bed
519
520        domain.set_quantity('elevation', zl*ones( (4,3) ))
521       
522       
523        domain.set_quantity('stage', [[val0, val0-1, val0-2],
524                                      [val1, val1+1, val1],
525                                      [val2, val2-2, val2],
526                                      [val3-0.5, val3, val3]])
527
528        domain.set_quantity('xmomentum',
529                            [[1, 2, 3], [3, 4, 5],
530                             [1, -1, 0], [0, -2, 2]])
531
532        domain.set_quantity('ymomentum',
533                            [[1, -1, 0], [0, -3, 2],
534                             [0, 1, 0], [-1, 2, 2]])               
535
536       
537        domain.check_integrity()
538
539
540
541        #Flux across right edge of volume 1
542        normal = domain.get_normal(1,0)       
543        ql = domain.get_conserved_quantities(vol_id=1, edge=0)
544        qr = domain.get_conserved_quantities(vol_id=2, edge=2)
545        flux0, max_speed = flux_function(normal, ql, qr, zl, zr)
546
547        #Flux across upper edge of volume 1
548        normal = domain.get_normal(1,1)       
549        ql = domain.get_conserved_quantities(vol_id=1, edge=1)
550        qr = domain.get_conserved_quantities(vol_id=3, edge=0)
551        flux1, max_speed = flux_function(normal, ql, qr, zl, zr)
552
553        #Flux across lower left hypotenuse of volume 1
554        normal = domain.get_normal(1,2)       
555        ql = domain.get_conserved_quantities(vol_id=1, edge=2)
556        qr = domain.get_conserved_quantities(vol_id=0, edge=1)
557        flux2, max_speed = flux_function(normal, ql, qr, zl, zr)
558       
559        #Scale, add up and check that compute_fluxes is correct for vol 1
560        e0 = domain.edgelengths[1, 0]
561        e1 = domain.edgelengths[1, 1]
562        e2 = domain.edgelengths[1, 2]
563
564        total_flux = -(e0*flux0+e1*flux1+e2*flux2)/domain.areas[1]
565
566       
567        domain.compute_fluxes()
568        for i, name in enumerate(['stage', 'xmomentum', 'ymomentum']):
569            assert allclose(total_flux[i],
570                            domain.quantities[name].explicit_update[1])       
571        #assert allclose(total_flux, domain.explicit_update[1,:])
572       
573
574    def test_compute_fluxes3(self):
575        #Random values, incl momentum
576       
577        a = [0.0, 0.0]
578        b = [0.0, 2.0]
579        c = [2.0,0.0]
580        d = [0.0, 4.0]
581        e = [2.0, 2.0]
582        f = [4.0,0.0]
583
584        points = [a, b, c, d, e, f]
585        #bac, bce, ecf, dbe
586        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
587       
588        domain = Domain(points, vertices)
589        val0 = 2.+2.0/3
590        val1 = 4.+4.0/3
591        val2 = 8.+2.0/3
592        val3 = 2.+8.0/3
593
594        zl=zr=-3.75 #Assume constant bed (must be less than stage)
595        domain.set_quantity('elevation', zl*ones( (4,3) ))
596       
597       
598        domain.set_quantity('stage', [[val0, val0-1, val0-2],
599                                      [val1, val1+1, val1],
600                                      [val2, val2-2, val2],
601                                      [val3-0.5, val3, val3]])
602
603        domain.set_quantity('xmomentum',
604                            [[1, 2, 3], [3, 4, 5],
605                             [1, -1, 0], [0, -2, 2]])
606
607        domain.set_quantity('ymomentum',
608                            [[1, -1, 0], [0, -3, 2],
609                             [0, 1, 0], [-1, 2, 2]])               
610
611       
612        domain.check_integrity()
613
614
615
616        #Flux across right edge of volume 1
617        normal = domain.get_normal(1,0)       
618        ql = domain.get_conserved_quantities(vol_id=1, edge=0)
619        qr = domain.get_conserved_quantities(vol_id=2, edge=2)
620        flux0, max_speed = flux_function(normal, ql, qr, zl, zr)
621
622        #Flux across upper edge of volume 1
623        normal = domain.get_normal(1,1)       
624        ql = domain.get_conserved_quantities(vol_id=1, edge=1)
625        qr = domain.get_conserved_quantities(vol_id=3, edge=0)
626        flux1, max_speed = flux_function(normal, ql, qr, zl, zr)
627
628        #Flux across lower left hypotenuse of volume 1
629        normal = domain.get_normal(1,2)       
630        ql = domain.get_conserved_quantities(vol_id=1, edge=2)
631        qr = domain.get_conserved_quantities(vol_id=0, edge=1)
632        flux2, max_speed = flux_function(normal, ql, qr, zl, zr)
633       
634        #Scale, add up and check that compute_fluxes is correct for vol 1
635        e0 = domain.edgelengths[1, 0]
636        e1 = domain.edgelengths[1, 1]
637        e2 = domain.edgelengths[1, 2]
638
639        total_flux = -(e0*flux0+e1*flux1+e2*flux2)/domain.areas[1]
640       
641        domain.compute_fluxes()
642        for i, name in enumerate(['stage', 'xmomentum', 'ymomentum']):
643            assert allclose(total_flux[i],
644                            domain.quantities[name].explicit_update[1])
645
646
647
648    def test_catching_negative_heights(self):
649       
650        a = [0.0, 0.0]
651        b = [0.0, 2.0]
652        c = [2.0,0.0]
653        d = [0.0, 4.0]
654        e = [2.0, 2.0]
655        f = [4.0,0.0]
656
657        points = [a, b, c, d, e, f]
658        #bac, bce, ecf, dbe
659        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
660       
661        domain = Domain(points, vertices)
662        val0 = 2.+2.0/3
663        val1 = 4.+4.0/3
664        val2 = 8.+2.0/3
665        val3 = 2.+8.0/3
666
667        zl=zr=4  #Too large
668        domain.set_quantity('elevation', zl*ones( (4,3) ))
669        domain.set_quantity('stage', [[val0, val0-1, val0-2],
670                                      [val1, val1+1, val1],
671                                      [val2, val2-2, val2],
672                                      [val3-0.5, val3, val3]])
673
674        #Should fail
675        try:
676            domain.check_integrity()
677        except:
678            pass
679
680
681
682
683    #####################################################   
684    def test_initial_condition(self):
685        """Test that initial condition is output at time == 0
686        """
687
688        from config import g
689        import copy
690       
691        a = [0.0, 0.0]
692        b = [0.0, 2.0]
693        c = [2.0, 0.0]
694        d = [0.0, 4.0]
695        e = [2.0, 2.0]
696        f = [4.0, 0.0]
697
698        points = [a, b, c, d, e, f]
699        #bac, bce, ecf, dbe
700        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
701       
702        domain = Domain(points, vertices)
703
704        #Set up for a gradient of (3,0) at mid triangle         
705        def slope(x, y):
706            return 3*x
707
708        h = 0.1
709        def stage(x,y):
710            return slope(x,y)+h
711
712        domain.set_quantity('elevation', slope)
713        domain.set_quantity('stage', stage)
714
715        initial_stage = copy.copy(domain.quantities['stage'].vertex_values)
716
717        domain.set_boundary({'exterior': Reflective_boundary(domain)})
718
719        #Evolution
720        for t in domain.evolve(yieldstep = 1.0, finaltime = 2.0):
721            stage = domain.quantities['stage'].vertex_values
722
723            if t == 0.0:
724                assert allclose(stage, initial_stage)
725            else:
726                assert not allclose(stage, initial_stage)
727               
728     
729       
730
731    #####################################################   
732    def test_gravity(self):
733        #Assuming no friction
734
735        from config import g
736       
737        a = [0.0, 0.0]
738        b = [0.0, 2.0]
739        c = [2.0, 0.0]
740        d = [0.0, 4.0]
741        e = [2.0, 2.0]
742        f = [4.0, 0.0]
743
744        points = [a, b, c, d, e, f]
745        #bac, bce, ecf, dbe
746        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
747       
748        domain = Domain(points, vertices)
749
750        #Set up for a gradient of (3,0) at mid triangle         
751        def slope(x, y):
752            return 3*x
753
754        h = 0.1
755        def stage(x,y):
756            return slope(x,y)+h
757
758        domain.set_quantity('elevation', slope)
759        domain.set_quantity('stage', stage)
760
761        for name in domain.conserved_quantities:
762            assert allclose(domain.quantities[name].explicit_update, 0)
763            assert allclose(domain.quantities[name].semi_implicit_update, 0)
764           
765        domain.compute_forcing_terms()                                   
766
767        assert allclose(domain.quantities['stage'].explicit_update, 0)
768        assert allclose(domain.quantities['xmomentum'].explicit_update, -g*h*3)
769        assert allclose(domain.quantities['ymomentum'].explicit_update, 0)
770       
771       
772    def test_manning_friction(self):
773        from config import g
774       
775        a = [0.0, 0.0]
776        b = [0.0, 2.0]
777        c = [2.0, 0.0]
778        d = [0.0, 4.0]
779        e = [2.0, 2.0]
780        f = [4.0, 0.0]
781
782        points = [a, b, c, d, e, f]
783        #bac, bce, ecf, dbe
784        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
785       
786        domain = Domain(points, vertices)
787
788        #Set up for a gradient of (3,0) at mid triangle         
789        def slope(x, y):
790            return 3*x
791
792        h = 0.1
793        def stage(x,y):
794            return slope(x,y)+h
795
796        eta = 0.07
797        domain.set_quantity('elevation', slope)
798        domain.set_quantity('stage', stage)
799        domain.set_quantity('friction', eta)       
800
801        for name in domain.conserved_quantities:
802            assert allclose(domain.quantities[name].explicit_update, 0)
803            assert allclose(domain.quantities[name].semi_implicit_update, 0)
804           
805        domain.compute_forcing_terms()                                   
806
807        assert allclose(domain.quantities['stage'].explicit_update, 0)
808        assert allclose(domain.quantities['xmomentum'].explicit_update, -g*h*3)
809        assert allclose(domain.quantities['ymomentum'].explicit_update, 0)
810
811        assert allclose(domain.quantities['stage'].semi_implicit_update, 0)
812        assert allclose(domain.quantities['xmomentum'].semi_implicit_update, 0)
813        assert allclose(domain.quantities['ymomentum'].semi_implicit_update, 0)       
814
815        #Create some momentum for friction to work with
816        domain.set_quantity('xmomentum', 1)
817        S = -g * eta**2 / h**(7.0/3)
818
819        domain.compute_forcing_terms()                                   
820        assert allclose(domain.quantities['stage'].semi_implicit_update, 0)
821        assert allclose(domain.quantities['xmomentum'].semi_implicit_update, S)
822        assert allclose(domain.quantities['ymomentum'].semi_implicit_update, 0)       
823
824        #A more complex example
825        domain.quantities['stage'].semi_implicit_update[:] = 0.0
826        domain.quantities['xmomentum'].semi_implicit_update[:] = 0.0
827        domain.quantities['ymomentum'].semi_implicit_update[:] = 0.0       
828       
829        domain.set_quantity('xmomentum', 3)
830        domain.set_quantity('ymomentum', 4)       
831
832        S = -g * eta**2 * 5 / h**(7.0/3)
833
834
835        domain.compute_forcing_terms()
836
837        assert allclose(domain.quantities['stage'].semi_implicit_update, 0)
838        assert allclose(domain.quantities['xmomentum'].semi_implicit_update, 3*S)
839        assert allclose(domain.quantities['ymomentum'].semi_implicit_update, 4*S)       
840
841    def test_constant_wind_stress(self):
842        from config import rho_a, rho_w, eta_w
843        from math import pi, cos, sin, sqrt
844
845        a = [0.0, 0.0]
846        b = [0.0, 2.0]
847        c = [2.0, 0.0]
848        d = [0.0, 4.0]
849        e = [2.0, 2.0]
850        f = [4.0, 0.0]
851
852        points = [a, b, c, d, e, f]
853        #bac, bce, ecf, dbe
854        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
855
856       
857        domain = Domain(points, vertices)
858
859        #Flat surface with 1m of water
860        domain.set_quantity('elevation', 0)
861        domain.set_quantity('stage', 1.0)
862        domain.set_quantity('friction', 0)       
863
864        Br = Reflective_boundary(domain)
865        domain.set_boundary({'exterior': Br})
866
867        #Setup only one forcing term, constant wind stress
868        s = 100
869        phi = 135
870        domain.forcing_terms = []
871        domain.forcing_terms.append( Wind_stress(s, phi) )
872       
873        domain.compute_forcing_terms()
874
875
876        const = eta_w*rho_a/rho_w
877       
878        #Convert to radians
879        phi = phi*pi/180
880       
881        #Compute velocity vector (u, v)
882        u = s*cos(phi)
883        v = s*sin(phi)
884
885        #Compute wind stress
886        S = const * sqrt(u**2 + v**2)
887
888        assert allclose(domain.quantities['stage'].explicit_update, 0)
889        assert allclose(domain.quantities['xmomentum'].explicit_update, S*u)
890        assert allclose(domain.quantities['ymomentum'].explicit_update, S*v)
891
892
893    def test_variable_wind_stress(self):
894        from config import rho_a, rho_w, eta_w
895        from math import pi, cos, sin, sqrt
896
897        a = [0.0, 0.0]
898        b = [0.0, 2.0]
899        c = [2.0, 0.0]
900        d = [0.0, 4.0]
901        e = [2.0, 2.0]
902        f = [4.0, 0.0]
903
904        points = [a, b, c, d, e, f]
905        #bac, bce, ecf, dbe
906        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
907
908        domain = Domain(points, vertices)
909
910        #Flat surface with 1m of water
911        domain.set_quantity('elevation', 0)
912        domain.set_quantity('stage', 1.0)
913        domain.set_quantity('friction', 0)       
914
915        Br = Reflective_boundary(domain)
916        domain.set_boundary({'exterior': Br})
917
918
919        domain.time = 5.54 #Take a random time (not zero)
920
921        #Setup only one forcing term, constant wind stress
922        s = 100
923        phi = 135
924        domain.forcing_terms = []
925        domain.forcing_terms.append( Wind_stress(s = speed, phi = angle) )
926       
927        domain.compute_forcing_terms()
928
929        #Compute reference solution
930        const = eta_w*rho_a/rho_w
931       
932        N = domain.number_of_elements
933
934        xc = domain.get_centroid_coordinates()
935        t = domain.time
936
937        x = xc[:,0]
938        y = xc[:,1]
939        s_vec = speed(t,x,y)
940        phi_vec = angle(t,x,y)
941       
942
943        for k in range(N):
944            #Convert to radians
945            phi = phi_vec[k]*pi/180
946            s = s_vec[k]
947       
948            #Compute velocity vector (u, v)
949            u = s*cos(phi)
950            v = s*sin(phi)
951
952            #Compute wind stress
953            S = const * sqrt(u**2 + v**2)
954
955            assert allclose(domain.quantities['stage'].explicit_update[k], 0)
956            assert allclose(domain.quantities['xmomentum'].explicit_update[k], S*u)
957            assert allclose(domain.quantities['ymomentum'].explicit_update[k], S*v)
958
959
960
961
962    def test_windfield_from_file(self):
963        from config import rho_a, rho_w, eta_w
964        from math import pi, cos, sin, sqrt
965        from config import time_format
966        from util import file_function
967        import time
968       
969
970        a = [0.0, 0.0]
971        b = [0.0, 2.0]
972        c = [2.0, 0.0]
973        d = [0.0, 4.0]
974        e = [2.0, 2.0]
975        f = [4.0, 0.0]
976
977        points = [a, b, c, d, e, f]
978        #bac, bce, ecf, dbe
979        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
980
981        domain = Domain(points, vertices)
982
983        #Flat surface with 1m of water
984        domain.set_quantity('elevation', 0)
985        domain.set_quantity('stage', 1.0)
986        domain.set_quantity('friction', 0)       
987
988        Br = Reflective_boundary(domain)
989        domain.set_boundary({'exterior': Br})
990
991
992        domain.time = 7 #Take a time that is represented in file (not zero)
993
994        #Write wind stress file (ensure that domain.time is covered)
995        #Take x=1 and y=0
996        filename = 'test_windstress_from_file.txt'
997        start = time.mktime(time.strptime('2000', '%Y'))       
998        fid = open(filename, 'w')
999        dt = 1  #One second interval
1000        t = 0.0
1001        while t <= 10.0:
1002            t_string = time.strftime(time_format, time.gmtime(t+start))
1003
1004            fid.write('%s, %f %f\n' %(t_string,
1005                                      speed(t,[1],[0])[0],
1006                                      angle(t,[1],[0])[0]))
1007            t += dt
1008   
1009        fid.close()
1010
1011        #Setup wind stress
1012        F = file_function(filename)
1013        W = Wind_stress(F) 
1014        domain.forcing_terms = []
1015        domain.forcing_terms.append(W)
1016       
1017        domain.compute_forcing_terms()
1018
1019        #Compute reference solution
1020        const = eta_w*rho_a/rho_w
1021       
1022        N = domain.number_of_elements
1023
1024        t = domain.time
1025
1026        s = speed(t,[1],[0])[0]
1027        phi = angle(t,[1],[0])[0]
1028       
1029        #Convert to radians
1030        phi = phi*pi/180
1031
1032       
1033        #Compute velocity vector (u, v)
1034        u = s*cos(phi)
1035        v = s*sin(phi)
1036
1037        #Compute wind stress
1038        S = const * sqrt(u**2 + v**2)
1039
1040        for k in range(N):
1041            assert allclose(domain.quantities['stage'].explicit_update[k], 0)
1042            assert allclose(domain.quantities['xmomentum'].explicit_update[k], S*u)
1043            assert allclose(domain.quantities['ymomentum'].explicit_update[k], S*v)
1044
1045        os.remove(filename)
1046
1047    def test_wind_stress_error_condition(self):
1048        """Test that windstress reacts properly when forcing functions
1049        are wrong - e.g. returns a scalar
1050        """
1051       
1052        from config import rho_a, rho_w, eta_w
1053        from math import pi, cos, sin, sqrt
1054
1055        a = [0.0, 0.0]
1056        b = [0.0, 2.0]
1057        c = [2.0, 0.0]
1058        d = [0.0, 4.0]
1059        e = [2.0, 2.0]
1060        f = [4.0, 0.0]
1061
1062        points = [a, b, c, d, e, f]
1063        #bac, bce, ecf, dbe
1064        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
1065
1066        domain = Domain(points, vertices)
1067
1068        #Flat surface with 1m of water
1069        domain.set_quantity('elevation', 0)
1070        domain.set_quantity('stage', 1.0)
1071        domain.set_quantity('friction', 0)       
1072
1073        Br = Reflective_boundary(domain)
1074        domain.set_boundary({'exterior': Br})
1075
1076
1077        domain.time = 5.54 #Take a random time (not zero)
1078
1079        #Setup only one forcing term, bad func
1080        domain.forcing_terms = []
1081
1082        try:
1083            domain.forcing_terms.append(Wind_stress(s = scalar_func,
1084                                                    phi = angle))
1085        except AssertionError:
1086            pass
1087        else:
1088            msg = 'Should have raised exception'
1089            raise msg
1090       
1091
1092        try:
1093            domain.forcing_terms.append(Wind_stress(s = speed,
1094                                                    phi = scalar_func))
1095        except AssertionError:
1096            pass
1097        else:
1098            msg = 'Should have raised exception'
1099            raise msg
1100       
1101        try:
1102            domain.forcing_terms.append(Wind_stress(s = speed,
1103                                                    phi = 'xx'))
1104        except:
1105            pass
1106        else:
1107            msg = 'Should have raised exception'
1108            raise msg
1109       
1110
1111    #####################################################
1112    def test_first_order_extrapolator_const_z(self):
1113       
1114        a = [0.0, 0.0]
1115        b = [0.0, 2.0]
1116        c = [2.0, 0.0]
1117        d = [0.0, 4.0]
1118        e = [2.0, 2.0]
1119        f = [4.0, 0.0]
1120
1121        points = [a, b, c, d, e, f]
1122        #bac, bce, ecf, dbe
1123        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
1124       
1125        domain = Domain(points, vertices)
1126        val0 = 2.+2.0/3
1127        val1 = 4.+4.0/3
1128        val2 = 8.+2.0/3
1129        val3 = 2.+8.0/3
1130
1131        zl=zr=-3.75 #Assume constant bed (must be less than stage)       
1132        domain.set_quantity('elevation', zl*ones( (4,3) ))
1133        domain.set_quantity('stage', [[val0, val0-1, val0-2],
1134                                      [val1, val1+1, val1],
1135                                      [val2, val2-2, val2],
1136                                      [val3-0.5, val3, val3]])
1137
1138
1139
1140        domain.order = 1
1141        domain.distribute_to_vertices_and_edges()
1142
1143        #Check that centroid values were distributed to vertices
1144        C = domain.quantities['stage'].centroid_values
1145        for i in range(3): 
1146            assert allclose( domain.quantities['stage'].vertex_values[:,i], C)
1147       
1148
1149    def test_first_order_limiter_variable_z(self):
1150        #Check that first order limiter follows bed_slope
1151        from Numeric import alltrue, greater_equal
1152        from config import epsilon
1153       
1154        a = [0.0, 0.0]
1155        b = [0.0, 2.0]
1156        c = [2.0,0.0]
1157        d = [0.0, 4.0]
1158        e = [2.0, 2.0]
1159        f = [4.0,0.0]
1160
1161        points = [a, b, c, d, e, f]
1162        #bac, bce, ecf, dbe
1163        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
1164       
1165        domain = Domain(points, vertices)
1166        val0 = 2.+2.0/3
1167        val1 = 4.+4.0/3
1168        val2 = 8.+2.0/3
1169        val3 = 2.+8.0/3
1170
1171        domain.set_quantity('elevation', [[0,0,0], [6,0,0],
1172                                          [6,6,6], [6,6,6]])
1173        domain.set_quantity('stage', [[val0, val0, val0],
1174                                      [val1, val1, val1],
1175                                      [val2, val2, val2],
1176                                      [val3, val3, val3]])
1177
1178        E = domain.quantities['elevation'].vertex_values
1179        L = domain.quantities['stage'].vertex_values       
1180
1181       
1182        #Check that some stages are not above elevation (within eps)
1183        #- so that the limiter has something to work with
1184        assert not alltrue(alltrue(greater_equal(L,E-epsilon)))               
1185
1186        domain.order = 1
1187        domain.distribute_to_vertices_and_edges()                           
1188
1189        #Check that all stages are above elevation (within eps)
1190        assert alltrue(alltrue(greater_equal(L,E-epsilon)))
1191
1192
1193    #####################################################   
1194    def test_distribute_basic(self):
1195        #Using test data generated by pyvolution-2
1196        #Assuming no friction and flat bed (0.0)
1197       
1198        a = [0.0, 0.0]
1199        b = [0.0, 2.0]
1200        c = [2.0, 0.0]
1201        d = [0.0, 4.0]
1202        e = [2.0, 2.0]
1203        f = [4.0, 0.0]
1204
1205        points = [a, b, c, d, e, f]
1206        #bac, bce, ecf, dbe
1207        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
1208       
1209        domain = Domain(points, vertices)
1210       
1211        val0 = 2.
1212        val1 = 4.
1213        val2 = 8.
1214        val3 = 2.
1215
1216        domain.set_quantity('stage', [val0, val1, val2, val3], 'centroids')
1217        L = domain.quantities['stage'].vertex_values
1218       
1219        #First order
1220        domain.order = 1       
1221        domain.distribute_to_vertices_and_edges()
1222        assert allclose(L[1], val1)
1223       
1224        #Second order
1225        domain.order = 2       
1226        domain.distribute_to_vertices_and_edges()
1227        assert allclose(L[1], [2.2, 4.9, 4.9])
1228               
1229
1230       
1231    def test_distribute_away_from_bed(self):           
1232        #Using test data generated by pyvolution-2
1233        #Assuming no friction and flat bed (0.0)
1234       
1235        a = [0.0, 0.0]
1236        b = [0.0, 2.0]
1237        c = [2.0, 0.0]
1238        d = [0.0, 4.0]
1239        e = [2.0, 2.0]
1240        f = [4.0, 0.0]
1241
1242        points = [a, b, c, d, e, f]
1243        #bac, bce, ecf, dbe
1244        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
1245       
1246        domain = Domain(points, vertices)
1247        L = domain.quantities['stage'].vertex_values   
1248       
1249        def stage(x,y):
1250            return x**2
1251       
1252        domain.set_quantity('stage', stage, 'centroids')
1253       
1254        a, b = domain.quantities['stage'].compute_gradients()           
1255        assert allclose(a[1], 3.33333334)
1256        assert allclose(b[1], 0.0)
1257       
1258        domain.order = 1
1259        domain.distribute_to_vertices_and_edges()
1260        assert allclose(L[1], 1.77777778)
1261       
1262        domain.order = 2
1263        domain.distribute_to_vertices_and_edges()
1264        assert allclose(L[1], [0.57777777, 2.37777778, 2.37777778])     
1265
1266
1267       
1268    def test_distribute_away_from_bed1(self):           
1269        #Using test data generated by pyvolution-2
1270        #Assuming no friction and flat bed (0.0)
1271       
1272        a = [0.0, 0.0]
1273        b = [0.0, 2.0]
1274        c = [2.0, 0.0]
1275        d = [0.0, 4.0]
1276        e = [2.0, 2.0]
1277        f = [4.0, 0.0]
1278
1279        points = [a, b, c, d, e, f]
1280        #bac, bce, ecf, dbe
1281        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
1282       
1283        domain = Domain(points, vertices)
1284        L = domain.quantities['stage'].vertex_values   
1285       
1286        def stage(x,y):
1287            return x**4+y**2
1288       
1289        domain.set_quantity('stage', stage, 'centroids')
1290        #print domain.quantities['stage'].centroid_values               
1291       
1292        a, b = domain.quantities['stage'].compute_gradients()           
1293        assert allclose(a[1], 25.18518519)
1294        assert allclose(b[1], 3.33333333)
1295       
1296        domain.order = 1
1297        domain.distribute_to_vertices_and_edges()
1298        assert allclose(L[1], 4.9382716)
1299       
1300        domain.order = 2
1301        domain.distribute_to_vertices_and_edges()
1302        assert allclose(L[1], [1.07160494, 6.46058131, 7.28262855])     
1303
1304   
1305   
1306    def test_distribute_near_bed(self):
1307        #Using test data generated by pyvolution-2
1308        #Assuming no friction and flat bed (0.0)
1309
1310        a = [0.0, 0.0]
1311        b = [0.0, 2.0]
1312        c = [2.0, 0.0]
1313        d = [0.0, 4.0]
1314        e = [2.0, 2.0]
1315        f = [4.0, 0.0]
1316
1317        points = [a, b, c, d, e, f]
1318        #bac, bce, ecf, dbe
1319        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
1320       
1321        domain = Domain(points, vertices)
1322
1323
1324        #Set up for a gradient of (3,0) at mid triangle         
1325        def slope(x, y):
1326            return 10*x
1327
1328        h = 0.1
1329        def stage(x,y):
1330            return slope(x,y)+h
1331
1332        domain.set_quantity('elevation', slope)
1333        domain.set_quantity('stage', stage, 'centroids')
1334
1335        #print domain.quantities['elevation'].centroid_values   
1336        #print domain.quantities['stage'].centroid_values
1337       
1338        E = domain.quantities['elevation'].vertex_values
1339        L = domain.quantities['stage'].vertex_values
1340
1341        #print E
1342        domain.order = 1       
1343        domain.distribute_to_vertices_and_edges()
1344        ##assert allclose(L[1], [0.19999999, 20.05, 20.05])
1345        assert allclose(L[1], [0.1, 20.1, 20.1])
1346       
1347        domain.order = 2       
1348        domain.distribute_to_vertices_and_edges()
1349        assert allclose(L[1], [0.1, 20.1, 20.1])
1350
1351    def test_distribute_near_bed1(self):
1352        #Using test data generated by pyvolution-2
1353        #Assuming no friction and flat bed (0.0)
1354
1355        a = [0.0, 0.0]
1356        b = [0.0, 2.0]
1357        c = [2.0, 0.0]
1358        d = [0.0, 4.0]
1359        e = [2.0, 2.0]
1360        f = [4.0, 0.0]
1361
1362        points = [a, b, c, d, e, f]
1363        #bac, bce, ecf, dbe
1364        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
1365       
1366        domain = Domain(points, vertices)
1367
1368
1369        #Set up for a gradient of (3,0) at mid triangle         
1370        def slope(x, y):
1371            return x**4+y**2
1372
1373        h = 0.1
1374        def stage(x,y):
1375            return slope(x,y)+h
1376
1377        domain.set_quantity('elevation', slope)
1378        domain.set_quantity('stage', stage)
1379
1380        #print domain.quantities['elevation'].centroid_values   
1381        #print domain.quantities['stage'].centroid_values
1382       
1383        E = domain.quantities['elevation'].vertex_values
1384        L = domain.quantities['stage'].vertex_values
1385
1386        #print E
1387        domain.order = 1
1388        domain.distribute_to_vertices_and_edges()
1389        ##assert allclose(L[1], [4.19999999, 16.07142857, 20.02857143])
1390        assert allclose(L[1], [4.1, 16.1, 20.1])       
1391
1392        domain.order = 2       
1393        domain.distribute_to_vertices_and_edges()
1394        assert allclose(L[1], [4.1, 16.1, 20.1])
1395
1396 
1397
1398    def test_second_order_distribute_real_data(self):
1399        #Using test data generated by pyvolution-2
1400        #Assuming no friction and flat bed (0.0)
1401
1402        a = [0.0, 0.0]
1403        b = [0.0, 1.0/5]
1404        c = [0.0, 2.0/5]
1405        d = [1.0/5, 0.0]
1406        e = [1.0/5, 1.0/5]
1407        f = [1.0/5, 2.0/5]
1408        g = [2.0/5, 2.0/5]
1409
1410        points = [a, b, c, d, e, f, g]
1411        #bae, efb, cbf, feg
1412        vertices = [ [1,0,4], [4,5,1], [2,1,5], [5,4,6]]
1413       
1414        domain = Domain(points, vertices)
1415
1416        def slope(x, y):
1417            return -x/3
1418
1419        domain.set_quantity('elevation', slope)
1420        domain.set_quantity('stage', 
1421                            [0.01298164, 0.00365611, 0.01440365, -0.0381856437096], 
1422                            'centroids')
1423        domain.set_quantity('xmomentum', 
1424                            [0.00670439, 0.01263789, 0.00647805, 0.0178180740668],
1425                            'centroids')
1426        domain.set_quantity('ymomentum',                           
1427                            [-7.23510980e-004, -6.30413883e-005, 6.30413883e-005, 0.000200907255866],
1428                            'centroids')                           
1429       
1430        E = domain.quantities['elevation'].vertex_values
1431        L = domain.quantities['stage'].vertex_values
1432        X = domain.quantities['xmomentum'].vertex_values       
1433        Y = domain.quantities['ymomentum'].vertex_values               
1434
1435        #print E
1436        domain.order = 2
1437        domain.beta_h = 0.0 #Use first order in h-limiter       
1438        domain.distribute_to_vertices_and_edges()
1439
1440        #print L[1,:]
1441        #print X[1,:]
1442        #print Y[1,:]
1443
1444        assert allclose(L[1,:], [-0.00825735775384, -0.00801881482869, 0.0272445025825])
1445        assert allclose(X[1,:], [0.0143507718962, 0.0142502147066, 0.00931268339717])
1446        assert allclose(Y[1,:], [-0.000117062180693, 7.94434448109e-005, -0.000151505429018])
1447
1448       
1449       
1450    def test_balance_deep_and_shallow(self):
1451        """Test that balanced limiters preserve conserved quantites.
1452        """
1453        import copy
1454       
1455        a = [0.0, 0.0]
1456        b = [0.0, 2.0]
1457        c = [2.0, 0.0]
1458        d = [0.0, 4.0]
1459        e = [2.0, 2.0]
1460        f = [4.0, 0.0]
1461
1462        points = [a, b, c, d, e, f]
1463
1464        #bac, bce, ecf, dbe
1465        elements = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4] ]
1466
1467        mesh = Domain(points, elements)       
1468        mesh.check_integrity()
1469
1470        #Create a deliberate overshoot
1471        mesh.set_quantity('stage', [[3,0,3], [2,2,6], [5,3,8], [8,3,5]])
1472        mesh.set_quantity('elevation', 0) #Flat bed
1473        stage = mesh.quantities['stage']       
1474       
1475        ref_centroid_values = copy.copy(stage.centroid_values[:]) #Copy
1476   
1477        #Limit
1478        mesh.distribute_to_vertices_and_edges()
1479
1480        #Assert that quantities are conserved
1481        from Numeric import sum
1482        for k in range(mesh.number_of_elements):
1483            assert allclose (ref_centroid_values[k],
1484                             sum(stage.vertex_values[k,:])/3)
1485       
1486
1487        #Now try with a non-flat bed - closely hugging initial stage in places
1488        #This will create alphas in the range [0, 0.478260, 1]
1489        mesh.set_quantity('stage', [[3,0,3], [2,2,6], [5,3,8], [8,3,5]])
1490        mesh.set_quantity('elevation', [[0,0,0], 
1491                                        [1.8,1.9,5.9], 
1492                                        [4.6,0,0], 
1493                                        [0,2,4]])
1494        stage = mesh.quantities['stage']       
1495       
1496        ref_centroid_values = copy.copy(stage.centroid_values[:]) #Copy
1497        ref_vertex_values = copy.copy(stage.vertex_values[:]) #Copy     
1498   
1499        #Limit
1500        mesh.distribute_to_vertices_and_edges()
1501
1502       
1503        #Assert that all vertex quantities have changed
1504        for k in range(mesh.number_of_elements):
1505            #print ref_vertex_values[k,:], stage.vertex_values[k,:]     
1506            assert not allclose (ref_vertex_values[k,:], stage.vertex_values[k,:])             
1507        #and assert that quantities are still conserved
1508        from Numeric import sum
1509        for k in range(mesh.number_of_elements):
1510            assert allclose (ref_centroid_values[k],
1511                             sum(stage.vertex_values[k,:])/3)
1512       
1513   
1514        #Also check that Python and C version produce the same
1515        assert allclose (stage.vertex_values,
1516                         [[2,2,2],
1517                          [1.93333333, 2.03333333, 6.03333333],
1518                          [6.93333333, 4.53333333, 4.53333333],
1519                          [5.33333333, 5.33333333, 5.33333333]])
1520                         
1521                         
1522
1523       
1524    def test_conservation_1(self):
1525        """Test that stage is conserved globally
1526       
1527        This one uses a flat bed, reflective bdries and a suitable
1528        initial condition
1529        """
1530        from mesh_factory import rectangular
1531        from shallow_water import Domain, Reflective_boundary,\
1532             Dirichlet_boundary, Constant_height
1533        from Numeric import array
1534
1535        #Create basic mesh
1536        points, vertices, boundary = rectangular(6, 6)
1537
1538        #Create shallow water domain
1539        domain = Domain(points, vertices, boundary)
1540        domain.smooth = False
1541        domain.default_order=1
1542
1543        #IC
1544        def x_slope(x, y):
1545            return x/3
1546
1547        domain.set_quantity('elevation', 0)
1548        domain.set_quantity('friction', 0)     
1549        domain.set_quantity('stage', x_slope)   
1550       
1551        # Boundary conditions (reflective everywhere)
1552        Br = Reflective_boundary(domain)
1553        domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
1554
1555        domain.check_integrity()
1556
1557        #domain.visualise = True #If you want to take a sticky beak 
1558       
1559        initial_volume = domain.quantities['stage'].get_integral()         
1560        initial_xmom = domain.quantities['xmomentum'].get_integral()     
1561       
1562        #print initial_xmom
1563
1564        #Evolution
1565        for t in domain.evolve(yieldstep = 0.05, finaltime = 5.0):
1566            volume =  domain.quantities['stage'].get_integral()         
1567            assert allclose (volume, initial_volume) 
1568         
1569            #I don't believe that the total momentum should be the same
1570            #It starts with zero and ends with zero though
1571            #xmom = domain.quantities['xmomentum'].get_integral()   
1572            #print xmom
1573            #assert allclose (xmom, initial_xmom)
1574       
1575
1576       
1577    def test_conservation_2(self):
1578        """Test that stage is conserved globally
1579       
1580        This one uses a slopy bed, reflective bdries and a suitable
1581        initial condition
1582        """
1583        from mesh_factory import rectangular
1584        from shallow_water import Domain, Reflective_boundary,\
1585             Dirichlet_boundary, Constant_height
1586        from Numeric import array
1587
1588        #Create basic mesh
1589        points, vertices, boundary = rectangular(6, 6)
1590
1591        #Create shallow water domain
1592        domain = Domain(points, vertices, boundary)
1593        domain.smooth = False
1594        domain.default_order=1
1595
1596        #IC
1597        def x_slope(x, y):
1598            return x/3
1599
1600        domain.set_quantity('elevation', x_slope)
1601        domain.set_quantity('friction', 0)     
1602        domain.set_quantity('stage', 0.4) #Steady       
1603       
1604        # Boundary conditions (reflective everywhere)
1605        Br = Reflective_boundary(domain)
1606        domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
1607
1608        domain.check_integrity()
1609
1610        #domain.visualise = True #If you want to take a sticky beak 
1611       
1612        initial_volume = domain.quantities['stage'].get_integral()         
1613        initial_xmom = domain.quantities['xmomentum'].get_integral()     
1614       
1615        #print initial_xmom
1616
1617        #Evolution
1618        for t in domain.evolve(yieldstep = 0.05, finaltime = 5.0):
1619            volume =  domain.quantities['stage'].get_integral()         
1620            assert allclose (volume, initial_volume) 
1621         
1622            #FIXME: What would we expect from momentum
1623            #xmom = domain.quantities['xmomentum'].get_integral()   
1624            #print xmom
1625            #assert allclose (xmom, initial_xmom)
1626       
1627       
1628    def test_conservation_3(self):
1629        """Test that stage is conserved globally
1630       
1631        This one uses a convoluted bed, reflective bdries and a suitable
1632        initial condition
1633        """
1634        from mesh_factory import rectangular
1635        from shallow_water import Domain, Reflective_boundary,\
1636             Dirichlet_boundary, Constant_height
1637        from Numeric import array
1638
1639        #Create basic mesh
1640        points, vertices, boundary = rectangular(6, 6)
1641
1642        #Create shallow water domain
1643        domain = Domain(points, vertices, boundary)
1644        domain.smooth = False
1645        domain.default_order=1
1646
1647        #IC
1648        def x_slope(x, y):
1649            z = 0*x
1650            for i in range(len(x)):
1651                if x[i] < 0.3: 
1652                    z[i] = x[i]/3
1653                if 0.3 <= x[i] < 0.5: 
1654                    z[i] = -0.5
1655                if 0.5 <= x[i] < 0.7: 
1656                    z[i] = 0.4       
1657                if 0.7 <= x[i]: 
1658                    z[i] = x[i]/3                   
1659            return z                               
1660                                   
1661                   
1662
1663        domain.set_quantity('elevation', x_slope)
1664        domain.set_quantity('friction', 0)     
1665        domain.set_quantity('stage', 0.4) #Steady       
1666       
1667        # Boundary conditions (reflective everywhere)
1668        Br = Reflective_boundary(domain)
1669        domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
1670
1671        domain.check_integrity()
1672
1673        #domain.visualise = True #If you want to take a sticky beak 
1674       
1675        initial_volume = domain.quantities['stage'].get_integral()         
1676        initial_xmom = domain.quantities['xmomentum'].get_integral()     
1677       
1678       
1679        #Let surface settle according to limiter
1680        #FIXME: What exactly is the significance of this?
1681        #Is it a problem that an arbitrary initial condition
1682        #needs to respect the law of the limiters?
1683        domain.distribute_to_vertices_and_edges()
1684        initial_volume = domain.quantities['stage'].get_integral()         
1685               
1686        #print initial_volume
1687       
1688        #print initial_xmom
1689
1690        #Evolution
1691        for t in domain.evolve(yieldstep = 0.05, finaltime = 5.0):
1692            volume =  domain.quantities['stage'].get_integral()         
1693            #print t, volume, initial_volume
1694           
1695            assert allclose (volume, initial_volume) 
1696         
1697            #FIXME: What would we expect from momentum
1698            #xmom = domain.quantities['xmomentum'].get_integral()   
1699            #print xmom
1700            #assert allclose (xmom, initial_xmom)
1701       
1702       
1703               
1704    def test_second_order_flat_bed_onestep(self):
1705
1706        from mesh_factory import rectangular
1707        from shallow_water import Domain, Reflective_boundary,\
1708             Dirichlet_boundary, Constant_height
1709        from Numeric import array
1710
1711        #Create basic mesh
1712        points, vertices, boundary = rectangular(6, 6)
1713
1714        #Create shallow water domain
1715        domain = Domain(points, vertices, boundary)
1716        domain.smooth = False
1717        domain.default_order=2
1718
1719        # Boundary conditions
1720        Br = Reflective_boundary(domain)
1721        Bd = Dirichlet_boundary([0.1, 0., 0.])
1722        domain.set_boundary({'left': Bd, 'right': Br, 'top': Br, 'bottom': Br})
1723
1724        domain.check_integrity()
1725
1726        #Evolution
1727        for t in domain.evolve(yieldstep = 0.05, finaltime = 0.05):
1728            pass# domain.write_time()
1729
1730        #Data from earlier version of pyvolution
1731        assert allclose(domain.min_timestep, 0.0396825396825)
1732        assert allclose(domain.max_timestep, 0.0396825396825) 
1733
1734        assert allclose(domain.quantities['stage'].centroid_values[:12],
1735                        [0.00171396, 0.02656103, 0.00241523, 0.02656103,
1736                        0.00241523, 0.02656103, 0.00241523, 0.02656103,
1737                        0.00241523, 0.02656103, 0.00241523, 0.0272623])
1738
1739        domain.distribute_to_vertices_and_edges()
1740        assert allclose(domain.quantities['stage'].vertex_values[:12,0],
1741                        [0.0001714, 0.02656103, 0.00024152,
1742                        0.02656103, 0.00024152, 0.02656103,
1743                        0.00024152, 0.02656103, 0.00024152,
1744                        0.02656103, 0.00024152, 0.0272623])
1745
1746        assert allclose(domain.quantities['stage'].vertex_values[:12,1],
1747                        [0.00315012, 0.02656103, 0.00024152, 0.02656103,
1748                         0.00024152, 0.02656103, 0.00024152, 0.02656103,
1749                         0.00024152, 0.02656103, 0.00040506, 0.0272623])
1750       
1751        assert allclose(domain.quantities['stage'].vertex_values[:12,2],
1752                        [0.00182037, 0.02656103, 0.00676264,
1753                         0.02656103, 0.00676264, 0.02656103,
1754                         0.00676264, 0.02656103, 0.00676264,
1755                         0.02656103, 0.0065991, 0.0272623])
1756
1757        assert allclose(domain.quantities['xmomentum'].centroid_values[:12],
1758                        [0.00113961, 0.01302432, 0.00148672,
1759                         0.01302432, 0.00148672, 0.01302432,
1760                         0.00148672, 0.01302432, 0.00148672 ,
1761                         0.01302432, 0.00148672, 0.01337143])
1762
1763        assert allclose(domain.quantities['ymomentum'].centroid_values[:12],
1764                        [-2.91240050e-004 , 1.22721531e-004,
1765                         -1.22721531e-004, 1.22721531e-004 ,
1766                         -1.22721531e-004, 1.22721531e-004,
1767                         -1.22721531e-004 , 1.22721531e-004,
1768                         -1.22721531e-004, 1.22721531e-004,
1769                         -1.22721531e-004, -4.57969873e-005])
1770
1771
1772       
1773    def test_second_order_flat_bed_moresteps(self):
1774
1775        from mesh_factory import rectangular
1776        from shallow_water import Domain, Reflective_boundary,\
1777             Dirichlet_boundary, Constant_height
1778        from Numeric import array
1779
1780        #Create basic mesh
1781        points, vertices, boundary = rectangular(6, 6)
1782
1783        #Create shallow water domain
1784        domain = Domain(points, vertices, boundary)
1785        domain.smooth = False
1786        domain.default_order=2
1787
1788        # Boundary conditions
1789        Br = Reflective_boundary(domain)
1790        Bd = Dirichlet_boundary([0.1, 0., 0.])
1791        domain.set_boundary({'left': Bd, 'right': Br, 'top': Br, 'bottom': Br})
1792
1793        domain.check_integrity()
1794
1795        #Evolution
1796        for t in domain.evolve(yieldstep = 0.05, finaltime = 0.1):
1797            pass
1798
1799        #Data from earlier version of pyvolution
1800        #assert allclose(domain.min_timestep, 0.0396825396825)
1801        #assert allclose(domain.max_timestep, 0.0396825396825)
1802        #print domain.quantities['stage'].centroid_values
1803
1804
1805    def test_flatbed_first_order(self):
1806        from mesh_factory import rectangular
1807        from shallow_water import Domain,\
1808             Reflective_boundary, Dirichlet_boundary
1809
1810        from Numeric import array
1811
1812        #Create basic mesh
1813        N = 8
1814        points, vertices, boundary = rectangular(N, N)
1815
1816        #Create shallow water domain
1817        domain = Domain(points, vertices, boundary)
1818        domain.smooth = False
1819        domain.visualise = False
1820        domain.default_order=1
1821
1822        # Boundary conditions
1823        Br = Reflective_boundary(domain)
1824        Bd = Dirichlet_boundary([0.2,0.,0.])
1825       
1826        domain.set_boundary({'left': Bd, 'right': Br, 'top': Br, 'bottom': Br})
1827        domain.check_integrity()
1828
1829 
1830        #Evolution
1831        for t in domain.evolve(yieldstep = 0.02, finaltime = 0.5):
1832            pass
1833            #domain.write_time()
1834
1835        #FIXME: These numbers were from version before 25/10       
1836        #assert allclose(domain.min_timestep, 0.0140413643926)
1837        #assert allclose(domain.max_timestep, 0.0140947355753)
1838
1839        for i in range(3):
1840            #assert allclose(domain.quantities['stage'].edge_values[:4,i],
1841            #                [0.10730244,0.12337617,0.11200126,0.12605666])
1842
1843            assert allclose(domain.quantities['xmomentum'].edge_values[:4,i],
1844                            [0.07610894,0.06901572,0.07284461,0.06819712])
1845
1846            #assert allclose(domain.quantities['ymomentum'].edge_values[:4,i],
1847            #                [-0.0060238, -0.00157404, -0.00309633, -0.0001637])                         
1848       
1849    def test_flatbed_second_order(self):
1850        from mesh_factory import rectangular
1851        from shallow_water import Domain,\
1852             Reflective_boundary, Dirichlet_boundary
1853
1854        from Numeric import array
1855
1856        #Create basic mesh
1857        N = 8
1858        points, vertices, boundary = rectangular(N, N)
1859
1860        #Create shallow water domain
1861        domain = Domain(points, vertices, boundary)
1862        domain.smooth = False
1863        domain.visualise = False
1864        domain.default_order=2
1865        #domain.minimum_allowed_height = 0.0 #Makes it like the 'oldstyle' balance
1866
1867        # Boundary conditions
1868        Br = Reflective_boundary(domain)
1869        Bd = Dirichlet_boundary([0.2,0.,0.])
1870
1871        domain.set_boundary({'left': Bd, 'right': Br, 'top': Br, 'bottom': Br})
1872        domain.check_integrity()
1873 
1874        #Evolution
1875        for t in domain.evolve(yieldstep = 0.01, finaltime = 0.03):
1876            pass
1877
1878
1879        assert allclose(domain.min_timestep, 0.0210448446782)
1880        assert allclose(domain.max_timestep, 0.0210448446782) 
1881
1882        #print domain.quantities['stage'].vertex_values[:4,0]
1883        #print domain.quantities['xmomentum'].vertex_values[:4,0]
1884        #print domain.quantities['ymomentum'].vertex_values[:4,0]
1885
1886        #FIXME: These numbers were from version before 25/10
1887        #assert allclose(domain.quantities['stage'].vertex_values[:4,0],
1888        #                [0.00101913,0.05352143,0.00104852,0.05354394])
1889       
1890        assert allclose(domain.quantities['xmomentum'].vertex_values[:4,0],
1891                        [ 0.00064835, 0.03685719, 0.00085073, 0.03687313])
1892
1893        #assert allclose(domain.quantities['xmomentum'].vertex_values[:4,0],
1894        #                [0.00090581,0.03685719,0.00088303,0.03687313])
1895       
1896        assert allclose(domain.quantities['ymomentum'].vertex_values[:4,0],
1897                        [-0.00139463,0.0006156,-0.00060364,0.00061827])
1898
1899
1900       
1901    def test_flatbed_second_order_distribute(self):
1902        #Use real data from pyvolution 2
1903        #painfully setup and extracted.
1904        from mesh_factory import rectangular
1905        from shallow_water import Domain,\
1906             Reflective_boundary, Dirichlet_boundary
1907
1908        from Numeric import array
1909
1910        #Create basic mesh
1911        N = 8
1912        points, vertices, boundary = rectangular(N, N)
1913
1914        #Create shallow water domain
1915        domain = Domain(points, vertices, boundary)
1916        domain.smooth = False
1917        domain.visualise = False
1918        domain.default_order=domain.order=2
1919
1920        # Boundary conditions
1921        Br = Reflective_boundary(domain)
1922        Bd = Dirichlet_boundary([0.2,0.,0.])
1923
1924        domain.set_boundary({'left': Bd, 'right': Br, 'top': Br, 'bottom': Br})
1925        domain.check_integrity()
1926
1927
1928
1929        for V in [False, True]:
1930            if V:
1931                #Set centroids as if system had been evolved           
1932                L = zeros(2*N*N, Float)
1933                L[:32] = [7.21205592e-003, 5.35214298e-002, 1.00910824e-002,
1934                          5.35439433e-002, 1.00910824e-002, 5.35439433e-002,
1935                          1.00910824e-002, 5.35439433e-002, 1.00910824e-002,
1936                          5.35439433e-002, 1.00910824e-002, 5.35439433e-002,
1937                          1.00910824e-002, 5.35393928e-002, 1.02344264e-002,
1938                          5.59605058e-002, 0.00000000e+000, 3.31027800e-004,
1939                          0.00000000e+000, 4.37962142e-005, 0.00000000e+000,
1940                          4.37962142e-005, 0.00000000e+000, 4.37962142e-005,
1941                          0.00000000e+000, 4.37962142e-005, 0.00000000e+000,
1942                          4.37962142e-005, 0.00000000e+000, 4.37962142e-005,
1943                          0.00000000e+000, 5.57305948e-005]
1944       
1945                X = zeros(2*N*N, Float)
1946                X[:32] = [6.48351607e-003, 3.68571894e-002, 8.50733285e-003,
1947                          3.68731327e-002, 8.50733285e-003, 3.68731327e-002,
1948                          8.50733285e-003, 3.68731327e-002, 8.50733285e-003,
1949                          3.68731327e-002, 8.50733285e-003, 3.68731327e-002,
1950                          8.50733285e-003, 3.68693861e-002, 8.65220973e-003,
1951                          3.85055387e-002, 0.00000000e+000, 2.86060840e-004,
1952                          0.00000000e+000, 3.58905503e-005, 0.00000000e+000,
1953                          3.58905503e-005, 0.00000000e+000, 3.58905503e-005,
1954                          0.00000000e+000, 3.58905503e-005, 0.00000000e+000,
1955                          3.58905503e-005, 0.00000000e+000, 3.58905503e-005,
1956                          0.00000000e+000, 4.57662812e-005]
1957
1958                Y = zeros(2*N*N, Float)
1959                Y[:32]=[-1.39463104e-003, 6.15600298e-004, -6.03637382e-004,
1960                        6.18272251e-004, -6.03637382e-004, 6.18272251e-004,
1961                        -6.03637382e-004, 6.18272251e-004, -6.03637382e-004,
1962                        6.18272251e-004, -6.03637382e-004, 6.18272251e-004,
1963                        -6.03637382e-004, 6.18599320e-004, -6.74622797e-004,
1964                        -1.48934756e-004, 0.00000000e+000, -5.35079969e-005,
1965                        0.00000000e+000, -2.57264987e-005, 0.00000000e+000,
1966                        -2.57264987e-005, 0.00000000e+000, -2.57264987e-005,
1967                        0.00000000e+000, -2.57264987e-005, 0.00000000e+000,
1968                        -2.57264987e-005, 0.00000000e+000, -2.57264987e-005,
1969                        0.00000000e+000, -2.57635178e-005]
1970
1971       
1972                domain.set_quantity('stage', L, 'centroids')
1973                domain.set_quantity('xmomentum', X, 'centroids')
1974                domain.set_quantity('ymomentum', Y, 'centroids')
1975       
1976                domain.check_integrity()
1977            else:   
1978                #Evolution
1979                for t in domain.evolve(yieldstep = 0.01, finaltime = 0.03):
1980                    pass
1981                assert allclose(domain.min_timestep, 0.0210448446782)
1982                assert allclose(domain.max_timestep, 0.0210448446782)
1983       
1984
1985            #Centroids were correct but not vertices.
1986            #Hence the check of distribute below.
1987            assert allclose(domain.quantities['stage'].centroid_values[:4],
1988                            [0.00721206,0.05352143,0.01009108,0.05354394])
1989
1990            assert allclose(domain.quantities['xmomentum'].centroid_values[:4],
1991                            [0.00648352,0.03685719,0.00850733,0.03687313])
1992
1993            assert allclose(domain.quantities['ymomentum'].centroid_values[:4],
1994                            [-0.00139463,0.0006156,-0.00060364,0.00061827])
1995
1996            #print 'C17=', domain.quantities['xmomentum'].centroid_values[17]
1997            #print 'C19=', domain.quantities['xmomentum'].centroid_values[19]
1998
1999            #assert allclose(domain.quantities['xmomentum'].centroid_values[17],0.00028606084)
2000            ##print domain.quantities['xmomentum'].centroid_values[17], V
2001            ##print
2002            if not V:           
2003                assert allclose(domain.quantities['xmomentum'].centroid_values[17], 0.0)
2004            else: 
2005                assert allclose(domain.quantities['xmomentum'].centroid_values[17], 0.00028606084)           
2006
2007            import copy
2008            XX = copy.copy(domain.quantities['xmomentum'].centroid_values)
2009            assert allclose(domain.quantities['xmomentum'].centroid_values, XX)
2010       
2011            domain.distribute_to_vertices_and_edges()
2012
2013            #assert allclose(domain.quantities['xmomentum'].centroid_values, XX)
2014
2015            assert allclose(domain.quantities['xmomentum'].centroid_values[17],
2016                            0.0)
2017
2018
2019            #FIXME: These numbers were from version before 25/10
2020            #assert allclose(domain.quantities['stage'].vertex_values[:4,0],
2021            #                [0.00101913,0.05352143,0.00104852,0.05354394])
2022
2023            assert allclose(domain.quantities['ymomentum'].vertex_values[:4,0],
2024                            [-0.00139463,0.0006156,-0.00060364,0.00061827])
2025
2026
2027            assert allclose(domain.quantities['xmomentum'].vertex_values[:4,0],
2028                            [0.00064835,0.03685719,0.00085073,0.03687313])
2029
2030
2031            #NB NO longer relvant:
2032
2033            #This was the culprit. First triangles vertex 0 had an
2034            #x-momentum of 0.0064835 instead of 0.00090581 and
2035            #third triangle had 0.00850733 instead of 0.00088303
2036            #print domain.quantities['xmomentum'].vertex_values[:4,0]
2037
2038            #print domain.quantities['xmomentum'].vertex_values[:4,0]
2039            #assert allclose(domain.quantities['xmomentum'].vertex_values[:4,0],
2040            #                [0.00090581,0.03685719,0.00088303,0.03687313])       
2041
2042
2043
2044       
2045       
2046    def test_bedslope_problem_first_order(self):
2047
2048        from mesh_factory import rectangular
2049        from shallow_water import Domain, Reflective_boundary, Constant_height
2050        from Numeric import array
2051
2052        #Create basic mesh
2053        points, vertices, boundary = rectangular(6, 6)
2054
2055        #Create shallow water domain
2056        domain = Domain(points, vertices, boundary)
2057        domain.smooth = False
2058        domain.default_order=1
2059
2060        #Bed-slope and friction
2061        def x_slope(x, y):
2062            return -x/3
2063
2064        domain.set_quantity('elevation', x_slope)
2065
2066        # Boundary conditions
2067        Br = Reflective_boundary(domain)
2068        domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
2069
2070        #Initial condition
2071        domain.set_quantity('stage', Constant_height(x_slope, 0.05))
2072        domain.check_integrity()
2073
2074        #Evolution
2075        for t in domain.evolve(yieldstep = 0.05, finaltime = 0.05):
2076            pass# domain.write_time()
2077
2078        assert allclose(domain.min_timestep, 0.050010003001)
2079        assert allclose(domain.max_timestep, 0.050010003001)
2080
2081       
2082    def test_bedslope_problem_first_order_moresteps(self):
2083
2084        from mesh_factory import rectangular
2085        from shallow_water import Domain, Reflective_boundary, Constant_height
2086        from Numeric import array
2087
2088        #Create basic mesh
2089        points, vertices, boundary = rectangular(6, 6)
2090
2091        #Create shallow water domain
2092        domain = Domain(points, vertices, boundary)
2093        domain.smooth = False
2094        domain.default_order=1
2095        domain.beta_h = 0.0 #Use first order in h-limiter
2096
2097        #Bed-slope and friction
2098        def x_slope(x, y):
2099            return -x/3
2100
2101        domain.set_quantity('elevation', x_slope)
2102
2103        # Boundary conditions
2104        Br = Reflective_boundary(domain)
2105        domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
2106
2107        #Initial condition
2108        domain.set_quantity('stage', Constant_height(x_slope, 0.05))
2109        domain.check_integrity()
2110
2111        #Evolution
2112        for t in domain.evolve(yieldstep = 0.05, finaltime = 0.5):
2113            pass# domain.write_time()
2114
2115        #Data from earlier version of pyvolution
2116        #print domain.quantities['stage'].centroid_values
2117       
2118        assert allclose(domain.quantities['stage'].centroid_values,
2119                        [-0.02998628, -0.01520652, -0.03043492,
2120                        -0.0149132, -0.03004706, -0.01476251,
2121                        -0.0298215, -0.01467976, -0.02988158,
2122                        -0.01474662, -0.03036161, -0.01442995,
2123                        -0.07624583, -0.06297061, -0.07733792,
2124                        -0.06342237, -0.07695439, -0.06289595,
2125                        -0.07635559, -0.0626065, -0.07633628,
2126                        -0.06280072, -0.07739632, -0.06386738,
2127                        -0.12161738, -0.11028239, -0.1223796,
2128                        -0.11095953, -0.12189744, -0.11048616,
2129                        -0.12074535, -0.10987605, -0.12014311,
2130                        -0.10976691, -0.12096859, -0.11087692,
2131                        -0.16868259, -0.15868061, -0.16801135,
2132                        -0.1588003, -0.16674343, -0.15813323,
2133                        -0.16457595, -0.15693826, -0.16281096,
2134                        -0.15585154, -0.16283873, -0.15540068,
2135                        -0.17450362, -0.19919913, -0.18062882,
2136                        -0.19764131, -0.17783111, -0.19407213,
2137                        -0.1736915, -0.19053624, -0.17228678,
2138                        -0.19105634, -0.17920133, -0.1968828,
2139                        -0.14244395, -0.14604641, -0.14473537,
2140                        -0.1506107, -0.14510055, -0.14919522,
2141                        -0.14175896, -0.14560798, -0.13911658,
2142                        -0.14439383, -0.13924047, -0.14829043])
2143                       
2144       
2145    def test_bedslope_problem_second_order_one_step(self):
2146
2147        from mesh_factory import rectangular
2148        from shallow_water import Domain, Reflective_boundary, Constant_height
2149        from Numeric import array
2150
2151        #Create basic mesh
2152        points, vertices, boundary = rectangular(6, 6)
2153
2154        #Create shallow water domain
2155        domain = Domain(points, vertices, boundary)
2156        domain.smooth = False
2157        domain.default_order=2
2158
2159        #Bed-slope and friction at vertices (and interpolated elsewhere)
2160        def x_slope(x, y):
2161            return -x/3
2162
2163        domain.set_quantity('elevation', x_slope)
2164
2165        # Boundary conditions
2166        Br = Reflective_boundary(domain)
2167        domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
2168
2169        #Initial condition
2170        domain.set_quantity('stage', Constant_height(x_slope, 0.05))
2171        domain.check_integrity()
2172
2173        assert allclose(domain.quantities['stage'].centroid_values,
2174                        [0.01296296, 0.03148148, 0.01296296,
2175                        0.03148148, 0.01296296, 0.03148148,
2176                        0.01296296, 0.03148148, 0.01296296,
2177                        0.03148148, 0.01296296, 0.03148148,
2178                        -0.04259259, -0.02407407, -0.04259259,
2179                        -0.02407407, -0.04259259, -0.02407407,
2180                        -0.04259259, -0.02407407, -0.04259259,
2181                        -0.02407407, -0.04259259, -0.02407407,
2182                        -0.09814815, -0.07962963, -0.09814815,
2183                        -0.07962963, -0.09814815, -0.07962963,
2184                        -0.09814815, -0.07962963, -0.09814815,
2185                        -0.07962963, -0.09814815, -0.07962963,
2186                        -0.1537037 , -0.13518519, -0.1537037,
2187                        -0.13518519, -0.1537037, -0.13518519,
2188                        -0.1537037 , -0.13518519, -0.1537037,
2189                        -0.13518519, -0.1537037, -0.13518519,
2190                        -0.20925926, -0.19074074, -0.20925926,
2191                        -0.19074074, -0.20925926, -0.19074074,
2192                        -0.20925926, -0.19074074, -0.20925926,
2193                        -0.19074074, -0.20925926, -0.19074074,
2194                        -0.26481481, -0.2462963, -0.26481481,
2195                        -0.2462963, -0.26481481, -0.2462963,
2196                        -0.26481481, -0.2462963, -0.26481481,
2197                        -0.2462963, -0.26481481, -0.2462963])
2198
2199
2200        #print domain.quantities['stage'].extrapolate_second_order()
2201        #domain.distribute_to_vertices_and_edges()
2202        #print domain.quantities['stage'].vertex_values[:,0]       
2203       
2204        #Evolution
2205        for t in domain.evolve(yieldstep = 0.05, finaltime = 0.05):
2206            #domain.write_time()           
2207            pass
2208
2209
2210        #print domain.quantities['stage'].centroid_values
2211        assert allclose(domain.quantities['stage'].centroid_values,
2212                 [0.01290985, 0.02356019, 0.01619096, 0.02356019, 0.01619096,
2213                  0.02356019, 0.01619096, 0.02356019, 0.01619096, 0.02356019,
2214                  0.01619096, 0.0268413, -0.04411074, -0.0248011, -0.04186556,
2215                  -0.0248011, -0.04186556, -0.0248011, -0.04186556, -0.0248011,
2216                  -0.04186556, -0.0248011, -0.04186556, -0.02255593,
2217                  -0.09966629, -0.08035666, -0.09742112, -0.08035666,
2218                  -0.09742112, -0.08035666, -0.09742112, -0.08035666,
2219                  -0.09742112, -0.08035666, -0.09742112, -0.07811149,
2220                  -0.15522185, -0.13591222, -0.15297667, -0.13591222,
2221                  -0.15297667, -0.13591222, -0.15297667, -0.13591222,
2222                  -0.15297667, -0.13591222, -0.15297667, -0.13366704,
2223                  -0.2107774, -0.19146777, -0.20853223, -0.19146777,
2224                  -0.20853223, -0.19146777, -0.20853223, -0.19146777,
2225                  -0.20853223, -0.19146777, -0.20853223, -0.1892226,
2226                  -0.26120669, -0.24776246, -0.25865535, -0.24776246,
2227                  -0.25865535, -0.24776246, -0.25865535, -0.24776246,
2228                  -0.25865535, -0.24776246, -0.25865535, -0.24521113])
2229
2230       
2231
2232    def test_bedslope_problem_second_order_two_steps(self):
2233
2234        from mesh_factory import rectangular
2235        from shallow_water import Domain, Reflective_boundary, Constant_height
2236        from Numeric import array
2237
2238        #Create basic mesh
2239        points, vertices, boundary = rectangular(6, 6)
2240
2241        #Create shallow water domain
2242        domain = Domain(points, vertices, boundary)
2243        domain.smooth = False
2244        domain.default_order=2
2245        domain.beta_h = 0.0 #Use first order in h-limiter       
2246
2247        #Bed-slope and friction at vertices (and interpolated elsewhere)
2248        def x_slope(x, y):
2249            return -x/3
2250
2251        domain.set_quantity('elevation', x_slope)
2252
2253        # Boundary conditions
2254        Br = Reflective_boundary(domain)
2255        domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
2256
2257        #Initial condition
2258        domain.set_quantity('stage', Constant_height(x_slope, 0.05))
2259        domain.check_integrity()
2260
2261        assert allclose(domain.quantities['stage'].centroid_values,
2262                        [0.01296296, 0.03148148, 0.01296296,
2263                        0.03148148, 0.01296296, 0.03148148,
2264                        0.01296296, 0.03148148, 0.01296296,
2265                        0.03148148, 0.01296296, 0.03148148,
2266                        -0.04259259, -0.02407407, -0.04259259,
2267                        -0.02407407, -0.04259259, -0.02407407,
2268                        -0.04259259, -0.02407407, -0.04259259,
2269                        -0.02407407, -0.04259259, -0.02407407,
2270                        -0.09814815, -0.07962963, -0.09814815,
2271                        -0.07962963, -0.09814815, -0.07962963,
2272                        -0.09814815, -0.07962963, -0.09814815,
2273                        -0.07962963, -0.09814815, -0.07962963,
2274                        -0.1537037 , -0.13518519, -0.1537037,
2275                        -0.13518519, -0.1537037, -0.13518519,
2276                        -0.1537037 , -0.13518519, -0.1537037,
2277                        -0.13518519, -0.1537037, -0.13518519,
2278                        -0.20925926, -0.19074074, -0.20925926,
2279                        -0.19074074, -0.20925926, -0.19074074,
2280                        -0.20925926, -0.19074074, -0.20925926,
2281                        -0.19074074, -0.20925926, -0.19074074,
2282                        -0.26481481, -0.2462963, -0.26481481,
2283                        -0.2462963, -0.26481481, -0.2462963,
2284                        -0.26481481, -0.2462963, -0.26481481,
2285                        -0.2462963, -0.26481481, -0.2462963])
2286
2287
2288        #print domain.quantities['stage'].extrapolate_second_order()
2289        #domain.distribute_to_vertices_and_edges()
2290        #print domain.quantities['stage'].vertex_values[:,0]       
2291       
2292        #Evolution
2293        for t in domain.evolve(yieldstep = 0.05, finaltime = 0.1):
2294            pass
2295
2296       
2297        #Data from earlier version of pyvolution ft=0.1
2298        assert allclose(domain.min_timestep, 0.0376895634803) 
2299        assert allclose(domain.max_timestep, 0.0415635655309)
2300
2301
2302        assert allclose(domain.quantities['stage'].centroid_values,
2303                        [0.00855788, 0.01575204, 0.00994606, 0.01717072,
2304                         0.01005985, 0.01716362, 0.01005985, 0.01716299,
2305                         0.01007098, 0.01736248, 0.01216452, 0.02026776,
2306                         -0.04462374, -0.02479045, -0.04199789, -0.0229465,
2307                         -0.04184033, -0.02295693, -0.04184013, -0.02295675,
2308                         -0.04184486, -0.0228168, -0.04028876, -0.02036486,
2309                         -0.10029444, -0.08170809, -0.09772846, -0.08021704,
2310                         -0.09760006, -0.08022143, -0.09759984, -0.08022124,
2311                         -0.09760261, -0.08008893, -0.09603914, -0.07758209,
2312                         -0.15584152, -0.13723138, -0.15327266, -0.13572906,
2313                         -0.15314427, -0.13573349, -0.15314405, -0.13573331,
2314                         -0.15314679, -0.13560104, -0.15158523, -0.13310701,
2315                         -0.21208605, -0.19283913, -0.20955631, -0.19134189,
2316                         -0.20942821, -0.19134598, -0.20942799, -0.1913458,
2317                         -0.20943005, -0.19120952, -0.20781177, -0.18869401,
2318                         -0.25384082, -0.2463294, -0.25047649, -0.24464654,
2319                         -0.25031159, -0.24464253, -0.25031112, -0.24464253,
2320                         -0.25031463, -0.24454764, -0.24885323, -0.24286438])
2321
2322       
2323       
2324               
2325    def test_bedslope_problem_second_order_two_yieldsteps(self):
2326
2327        from mesh_factory import rectangular
2328        from shallow_water import Domain, Reflective_boundary, Constant_height
2329        from Numeric import array
2330
2331        #Create basic mesh
2332        points, vertices, boundary = rectangular(6, 6)
2333
2334        #Create shallow water domain
2335        domain = Domain(points, vertices, boundary)
2336        domain.smooth = False
2337        domain.default_order=2
2338        domain.beta_h = 0.0 #Use first order in h-limiter       
2339
2340        #Bed-slope and friction at vertices (and interpolated elsewhere)
2341        def x_slope(x, y):
2342            return -x/3
2343
2344        domain.set_quantity('elevation', x_slope)
2345
2346        # Boundary conditions
2347        Br = Reflective_boundary(domain)
2348        domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
2349
2350        #Initial condition
2351        domain.set_quantity('stage', Constant_height(x_slope, 0.05))
2352        domain.check_integrity()
2353
2354        assert allclose(domain.quantities['stage'].centroid_values,
2355                        [0.01296296, 0.03148148, 0.01296296,
2356                        0.03148148, 0.01296296, 0.03148148,
2357                        0.01296296, 0.03148148, 0.01296296,
2358                        0.03148148, 0.01296296, 0.03148148,
2359                        -0.04259259, -0.02407407, -0.04259259,
2360                        -0.02407407, -0.04259259, -0.02407407,
2361                        -0.04259259, -0.02407407, -0.04259259,
2362                        -0.02407407, -0.04259259, -0.02407407,
2363                        -0.09814815, -0.07962963, -0.09814815,
2364                        -0.07962963, -0.09814815, -0.07962963,
2365                        -0.09814815, -0.07962963, -0.09814815,
2366                        -0.07962963, -0.09814815, -0.07962963,
2367                        -0.1537037 , -0.13518519, -0.1537037,
2368                        -0.13518519, -0.1537037, -0.13518519,
2369                        -0.1537037 , -0.13518519, -0.1537037,
2370                        -0.13518519, -0.1537037, -0.13518519,
2371                        -0.20925926, -0.19074074, -0.20925926,
2372                        -0.19074074, -0.20925926, -0.19074074,
2373                        -0.20925926, -0.19074074, -0.20925926,
2374                        -0.19074074, -0.20925926, -0.19074074,
2375                        -0.26481481, -0.2462963, -0.26481481,
2376                        -0.2462963, -0.26481481, -0.2462963,
2377                        -0.26481481, -0.2462963, -0.26481481,
2378                        -0.2462963, -0.26481481, -0.2462963])
2379
2380
2381        #print domain.quantities['stage'].extrapolate_second_order()
2382        #domain.distribute_to_vertices_and_edges()
2383        #print domain.quantities['stage'].vertex_values[:,0]       
2384       
2385        #Evolution
2386        for t in domain.evolve(yieldstep = 0.05, finaltime = 0.1):   #0.05??
2387            #domain.write_time()           
2388            pass
2389
2390
2391       
2392        assert allclose(domain.quantities['stage'].centroid_values, 
2393                 [0.00855788, 0.01575204, 0.00994606, 0.01717072, 0.01005985,
2394                  0.01716362, 0.01005985, 0.01716299, 0.01007098, 0.01736248,
2395                  0.01216452, 0.02026776, -0.04462374, -0.02479045, -0.04199789,
2396                  -0.0229465, -0.04184033, -0.02295693, -0.04184013,
2397                  -0.02295675, -0.04184486, -0.0228168, -0.04028876,
2398                  -0.02036486, -0.10029444, -0.08170809, -0.09772846,
2399                  -0.08021704, -0.09760006, -0.08022143, -0.09759984,
2400                  -0.08022124, -0.09760261, -0.08008893, -0.09603914,
2401                  -0.07758209, -0.15584152, -0.13723138, -0.15327266,
2402                  -0.13572906, -0.15314427, -0.13573349, -0.15314405,
2403                  -0.13573331, -0.15314679, -0.13560104, -0.15158523,
2404                  -0.13310701, -0.21208605, -0.19283913, -0.20955631,
2405                  -0.19134189, -0.20942821, -0.19134598, -0.20942799,
2406                  -0.1913458, -0.20943005, -0.19120952, -0.20781177,
2407                  -0.18869401, -0.25384082, -0.2463294, -0.25047649,
2408                  -0.24464654, -0.25031159, -0.24464253, -0.25031112,
2409                  -0.24464253, -0.25031463, -0.24454764, -0.24885323,
2410                  -0.24286438])
2411
2412
2413
2414    def test_bedslope_problem_second_order_more_steps(self):
2415
2416        from mesh_factory import rectangular
2417        from shallow_water import Domain, Reflective_boundary, Constant_height
2418        from Numeric import array
2419
2420        #Create basic mesh
2421        points, vertices, boundary = rectangular(6, 6)
2422
2423        #Create shallow water domain
2424        domain = Domain(points, vertices, boundary)
2425        domain.smooth = False
2426        domain.default_order=2
2427        domain.beta_h = 0.0 #Use first order in h-limiter       
2428
2429        #Bed-slope and friction at vertices (and interpolated elsewhere)
2430        def x_slope(x, y):
2431            return -x/3
2432
2433        domain.set_quantity('elevation', x_slope)
2434
2435        # Boundary conditions
2436        Br = Reflective_boundary(domain)
2437        domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
2438
2439        #Initial condition
2440        domain.set_quantity('stage', Constant_height(x_slope, 0.05))
2441        domain.check_integrity()
2442
2443        assert allclose(domain.quantities['stage'].centroid_values,
2444                        [0.01296296, 0.03148148, 0.01296296,
2445                        0.03148148, 0.01296296, 0.03148148,
2446                        0.01296296, 0.03148148, 0.01296296,
2447                        0.03148148, 0.01296296, 0.03148148,
2448                        -0.04259259, -0.02407407, -0.04259259,
2449                        -0.02407407, -0.04259259, -0.02407407,
2450                        -0.04259259, -0.02407407, -0.04259259,
2451                        -0.02407407, -0.04259259, -0.02407407,
2452                        -0.09814815, -0.07962963, -0.09814815,
2453                        -0.07962963, -0.09814815, -0.07962963,
2454                        -0.09814815, -0.07962963, -0.09814815,
2455                        -0.07962963, -0.09814815, -0.07962963,
2456                        -0.1537037 , -0.13518519, -0.1537037,
2457                        -0.13518519, -0.1537037, -0.13518519,
2458                        -0.1537037 , -0.13518519, -0.1537037,
2459                        -0.13518519, -0.1537037, -0.13518519,
2460                        -0.20925926, -0.19074074, -0.20925926,
2461                        -0.19074074, -0.20925926, -0.19074074,
2462                        -0.20925926, -0.19074074, -0.20925926,
2463                        -0.19074074, -0.20925926, -0.19074074,
2464                        -0.26481481, -0.2462963, -0.26481481,
2465                        -0.2462963, -0.26481481, -0.2462963,
2466                        -0.26481481, -0.2462963, -0.26481481,
2467                        -0.2462963, -0.26481481, -0.2462963])
2468
2469
2470        #print domain.quantities['stage'].extrapolate_second_order()
2471        #domain.distribute_to_vertices_and_edges()
2472        #print domain.quantities['stage'].vertex_values[:,0]       
2473       
2474        #Evolution
2475        for t in domain.evolve(yieldstep = 0.05, finaltime = 0.5):
2476            pass
2477
2478       
2479        assert allclose(domain.quantities['stage'].centroid_values,
2480     [-0.02907028, -0.01475478, -0.02973417, -0.01447186, -0.02932665, -0.01428285,
2481      -0.02901975, -0.0141361,  -0.02898816, -0.01418135, -0.02961409, -0.01403487,
2482      -0.07597998, -0.06252591, -0.07664854, -0.06312532, -0.07638287, -0.06265139,
2483      -0.07571145, -0.06235231, -0.0756817,  -0.06245309, -0.07652292, -0.06289946,
2484      -0.12367464, -0.11088981, -0.12237277, -0.11115338, -0.1218934,  -0.1107174, 
2485      -0.12081485, -0.11000491, -0.12038451, -0.11010335, -0.12102113, -0.11012105,
2486      -0.16909116, -0.15831543, -0.16730214, -0.15786249, -0.1665493,  -0.15697919,
2487      -0.16496618, -0.15559852, -0.16338679, -0.15509088, -0.16364092, -0.15424423,
2488      -0.18771107, -0.19903904, -0.18903759, -0.19858437, -0.18701552, -0.19697797,
2489      -0.1833593,  -0.19505871, -0.1818806,  -0.19418042, -0.18586159, -0.19576946,
2490      -0.13986873, -0.14170053, -0.14132188, -0.14560674, -0.14095617, -0.14373292,
2491      -0.13785933, -0.14033364, -0.13592955, -0.13936356, -0.13596008, -0.14216296])
2492
2493        assert allclose(domain.quantities['xmomentum'].centroid_values,
2494     [ 0.00831121,  0.00317948,  0.00731797,  0.00334939,  0.00764717,  0.00348053,
2495       0.00788729,  0.00356522,  0.00780649,  0.00341919,  0.00693525,  0.00310375,
2496       0.02166196,  0.01421475,  0.02017737,  0.01316839,  0.02037015,  0.01368659,
2497       0.02106,     0.01399161,  0.02074514,  0.01354935,  0.01887407,  0.0123113, 
2498       0.03775083,  0.02855197,  0.03689337,  0.02759782,  0.03732848,  0.02812072,
2499       0.03872545,  0.02913348,  0.03880939,  0.02803804,  0.03546499,  0.0260039, 
2500       0.0632131,   0.04730634,  0.0576324,   0.04592336,  0.05790921,  0.04690514,
2501       0.05986467,  0.04871165,  0.06170068,  0.04811572,  0.05657041,  0.04416292,
2502       0.08489642,  0.07188097,  0.07835261,  0.06843406,  0.07986412,  0.0698247, 
2503       0.08201071,  0.07216756,  0.08378418,  0.07273624,  0.080399,    0.06645841,
2504       0.01631548,  0.04691608,  0.0206632,   0.044409,    0.02115518,  0.04560305,
2505       0.02160608,  0.04663725,  0.02174734,  0.04795559,  0.02281427,  0.05667111])
2506
2507       
2508        assert allclose(domain.quantities['ymomentum'].centroid_values, 
2509     [ 1.45876601e-004, -3.24627393e-004, -1.57572719e-004, -2.92790187e-004,
2510      -9.90988382e-005, -3.06677335e-004, -1.62493106e-004, -3.71310004e-004,
2511      -1.99445058e-004, -3.28493467e-004,  6.68217349e-005, -8.42042805e-006,
2512       5.05093371e-004, -1.42842214e-004, -6.81454718e-005, -5.02084057e-004,
2513      -8.50583861e-005, -4.65443981e-004, -1.96406564e-004, -5.88889562e-004,
2514      -2.70160173e-004, -5.35485454e-004,  2.60780997e-004,  3.12145471e-005,
2515       5.16189608e-004,  1.07069062e-004,  9.29989252e-005, -3.71211119e-004,
2516       1.16350246e-004, -3.82407830e-004, -1.62077969e-004, -6.30906636e-004,
2517      -4.74025708e-004, -6.94463009e-004,  6.15092843e-005,  2.22106820e-004,
2518      -6.29589294e-004,  2.43611937e-004, -5.88125094e-004, -6.94293192e-005,
2519      -4.17914641e-004,  6.64609019e-005, -7.68334577e-004, -3.40232101e-004,
2520      -1.67424308e-003, -7.39485066e-004, -1.59966988e-003,  5.68262838e-005,
2521      -1.48470633e-003, -1.84554882e-003, -2.27200099e-003, -1.67506848e-003,
2522      -1.95610258e-003, -1.47638801e-003, -1.73779477e-003, -1.85498791e-003,
2523      -2.01357843e-003, -2.17675471e-003, -1.65783870e-003, -1.15818681e-003,
2524      -1.18663036e-003, -2.94229849e-003, -3.59309018e-003, -5.13496584e-003,
2525      -6.17359400e-003, -5.98761937e-003, -6.00540116e-003, -5.01121966e-003,
2526      -4.50964850e-003, -3.06319963e-003,  6.08950810e-004, -4.79537921e-004])
2527       
2528
2529       
2530
2531    def test_temp_play(self):
2532
2533        from mesh_factory import rectangular
2534        from shallow_water import Domain, Reflective_boundary, Constant_height
2535        from Numeric import array
2536
2537        #Create basic mesh
2538        points, vertices, boundary = rectangular(5, 5)
2539
2540        #Create shallow water domain
2541        domain = Domain(points, vertices, boundary)
2542        domain.smooth = False
2543        domain.default_order=2
2544        domain.beta_h = 0.0 #Use first order in h-limiter       
2545
2546        #Bed-slope and friction at vertices (and interpolated elsewhere)
2547        def x_slope(x, y):
2548            return -x/3
2549
2550        domain.set_quantity('elevation', x_slope)
2551
2552        # Boundary conditions
2553        Br = Reflective_boundary(domain)
2554        domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
2555
2556        #Initial condition
2557        domain.set_quantity('stage', Constant_height(x_slope, 0.05))
2558        domain.check_integrity()
2559
2560        #Evolution
2561        for t in domain.evolve(yieldstep = 0.05, finaltime = 0.1):
2562            pass
2563
2564        assert allclose(domain.quantities['stage'].centroid_values[:4],
2565                        [0.00206836, 0.01296714, 0.00363415, 0.01438924])
2566        assert allclose(domain.quantities['xmomentum'].centroid_values[:4],                     
2567                        [0.01360154, 0.00671133, 0.01264578, 0.00648503])
2568        assert allclose(domain.quantities['ymomentum'].centroid_values[:4],
2569                        [-1.19201077e-003, -7.23647546e-004, 
2570                        -6.39083123e-005, 6.29815168e-005])
2571       
2572       
2573    def test_complex_bed(self):
2574        #No friction is tested here
2575       
2576        from shallow_water import Domain, Reflective_boundary, Dirichlet_boundary,\
2577             Transmissive_boundary, Time_boundary,\
2578             Weir_simple as Weir, Constant_height
2579
2580        from mesh_factory import rectangular
2581        from Numeric import array
2582   
2583        N = 12
2584        points, vertices, boundary = rectangular(N, N/2, len1=1.2,len2=0.6,
2585                                                 origin=(-0.07, 0))
2586
2587
2588        domain = Domain(points, vertices, boundary)
2589        domain.smooth = False
2590        domain.visualise = False
2591        domain.default_order=2
2592
2593       
2594        inflow_stage = 0.1
2595        Z = Weir(inflow_stage)
2596        domain.set_quantity('elevation', Z)
2597
2598        Br = Reflective_boundary(domain)
2599        Bd = Dirichlet_boundary([inflow_stage, 0.0, 0.0])
2600        domain.set_boundary({'left': Bd, 'right': Br, 'bottom': Br, 'top': Br})
2601
2602        domain.set_quantity('stage', Constant_height(Z, 0.))
2603
2604        for t in domain.evolve(yieldstep = 0.02, finaltime = 0.2):
2605            pass
2606
2607
2608        #print domain.quantities['stage'].centroid_values
2609
2610        #FIXME: These numbers were from version before 25/10       
2611        #assert allclose(domain.quantities['stage'].centroid_values,
2612# [3.95822638e-002,  5.61022588e-002,  4.66437868e-002,  5.73081011e-002,
2613#  4.72394613e-002,  5.74684939e-002,  4.74309483e-002,  5.77458084e-002,
2614#  4.80628177e-002,  5.85656225e-002,  4.90498542e-002,  6.02609831e-002,
2615#  1.18470315e-002,  1.75136443e-002,  1.18035266e-002,  2.15565695e-002,
2616#  1.31620268e-002,  2.14351640e-002,  1.32351076e-002,  2.15450687e-002,
2617#  1.36414028e-002,  2.24274619e-002,  1.51689511e-002,  2.21789655e-002,
2618# -7.54337535e-003, -6.86362021e-004, -7.74146760e-003, -1.83756530e-003,
2619# -8.16773628e-003, -4.49916813e-004, -8.08202599e-003, -3.91118720e-004,
2620# -8.10292716e-003, -3.88584984e-004, -7.35226124e-003,  2.73985295e-004,
2621#  1.86166683e-001,  8.74070369e-002,  1.86166712e-001,  8.74035875e-002,
2622#  6.11666935e-002, -3.76173225e-002, -6.38333276e-002, -3.76147365e-002,
2623#  6.11666725e-002,  8.73846774e-002,  1.86166697e-001,  8.74171550e-002,
2624# -4.83333333e-002,  1.18333333e-001, -4.83333333e-002,  1.18333333e-001,
2625# -4.83333333e-002, -6.66666667e-003, -1.73333333e-001, -1.31666667e-001,
2626# -1.73333333e-001, -6.66666667e-003, -4.83333333e-002,  1.18333333e-001,
2627# -2.48333333e-001, -2.31666667e-001, -2.48333333e-001, -2.31666667e-001,
2628# -2.48333333e-001, -2.31666667e-001, -2.48333333e-001, -2.31666667e-001,
2629# -2.48333333e-001, -2.31666667e-001, -2.48333333e-001, -2.31666667e-001,
2630# -4.65000000e-001, -3.65000000e-001, -4.65000000e-001, -3.65000000e-001,
2631# -4.65000000e-001, -3.65000000e-001, -4.65000000e-001, -3.65000000e-001,
2632# -4.65000000e-001, -3.65000000e-001, -4.65000000e-001, -3.65000000e-001,
2633# -5.98333333e-001, -5.81666667e-001, -5.98333333e-001, -5.81666667e-001,
2634# -5.98333333e-001, -5.81666667e-001, -5.98333333e-001, -5.81666667e-001,
2635# -5.98333333e-001, -5.81666667e-001, -5.98333333e-001, -5.81666667e-001,
2636# -6.48333333e-001, -6.31666667e-001, -6.48333333e-001, -6.31666667e-001,
2637# -6.48333333e-001, -6.31666667e-001, -6.48333333e-001, -6.31666667e-001,
2638# -6.48333333e-001, -6.31666667e-001, -6.48333333e-001, -6.31666667e-001,
2639# -5.31666667e-001, -5.98333333e-001, -5.31666667e-001, -5.98333333e-001,
2640# -5.31666667e-001, -5.98333333e-001, -5.31666667e-001, -5.98333333e-001,
2641# -5.31666667e-001, -5.98333333e-001, -5.31666667e-001, -5.98333333e-001,
2642# -4.98333333e-001, -4.81666667e-001, -4.98333333e-001, -4.81666667e-001,
2643# -4.98333333e-001, -4.81666667e-001, -4.98333333e-001, -4.81666667e-001,
2644# -4.98333333e-001, -4.81666667e-001, -4.98333333e-001, -4.81666667e-001,
2645# -5.48333333e-001, -5.31666667e-001, -5.48333333e-001, -5.31666667e-001,
2646# -5.48333333e-001, -5.31666667e-001, -5.48333333e-001, -5.31666667e-001,
2647# -5.48333333e-001, -5.31666667e-001, -5.48333333e-001, -5.31666667e-001])
2648
2649
2650
2651    def test_spatio_temporal_boundary(self):
2652        """Test that boundary values can be read from file and interpolated
2653        in both time and space.
2654
2655        Verify that the same steady state solution is arrived at and that
2656        time interpolation works.
2657
2658        The full solution history is not exactly the same as
2659        file boundary must read and interpolate from *smoothed* version
2660        as stored in sww.
2661        """
2662        import time
2663       
2664        #Create sww file of simple propagation from left to right
2665        #through rectangular domain
2666       
2667        from mesh_factory import rectangular
2668
2669        #Create basic mesh
2670        points, vertices, boundary = rectangular(3, 3)
2671
2672        #Create shallow water domain
2673        domain1 = Domain(points, vertices, boundary)
2674       
2675        from util import mean
2676        domain1.reduction = mean
2677        domain1.smooth = False  #Exact result
2678
2679        domain1.default_order = 2
2680        domain1.store = True   
2681        domain1.set_datadir('.')
2682        domain1.set_name('spatio_temporal_boundary_source' + str(time.time()))
2683       
2684        #FIXME: This is extremely important!
2685        #How can we test if they weren't stored?
2686        domain1.quantities_to_be_stored = ['stage', 'xmomentum', 'ymomentum']
2687
2688               
2689        #Bed-slope and friction at vertices (and interpolated elsewhere)
2690        domain1.set_quantity('elevation', 0)
2691        domain1.set_quantity('friction', 0)     
2692
2693        # Boundary conditions
2694        Br = Reflective_boundary(domain1)
2695        Bd = Dirichlet_boundary([0.3,0,0])
2696        domain1.set_boundary({'left': Bd, 'top': Bd, 'right': Br, 'bottom': Br})
2697        #Initial condition
2698        domain1.set_quantity('stage', 0)
2699        domain1.check_integrity()
2700
2701        finaltime = 5
2702        #Evolution  (full domain - large steps)
2703        for t in domain1.evolve(yieldstep = 0.671, finaltime = finaltime):
2704            pass
2705            #domain1.write_time()
2706
2707        cv1 = domain1.quantities['stage'].centroid_values
2708       
2709       
2710        #Create an triangle shaped domain (reusing coordinates from domain 1),
2711        #formed from the lower and right hand  boundaries and
2712        #the sw-ne diagonal
2713        #from domain 1. Call it domain2
2714       
2715        points = [ [0,0], [1.0/3,0], [1.0/3,1.0/3],
2716                   [2.0/3,0], [2.0/3,1.0/3], [2.0/3,2.0/3],
2717                   [1,0], [1,1.0/3], [1,2.0/3], [1,1]]
2718       
2719        vertices = [ [1,2,0], 
2720                     [3,4,1], [2,1,4], [4,5,2], 
2721                     [6,7,3], [4,3,7], [7,8,4], [5,4,8], [8,9,5]]       
2722
2723        boundary = { (0,1):'bottom', (1,1):'bottom', (4,1): 'bottom',
2724                     (4,2):'right', (6,2):'right', (8,2):'right',
2725                     (0,0):'diagonal', (3,0):'diagonal', (8,0):'diagonal'}   
2726                     
2727        domain2 = Domain(points, vertices, boundary)
2728       
2729        domain2.reduction = domain1.reduction
2730        domain2.smooth = False
2731        domain2.default_order = 2
2732       
2733        #Bed-slope and friction at vertices (and interpolated elsewhere)
2734        domain2.set_quantity('elevation', 0)
2735        domain2.set_quantity('friction', 0)
2736        domain2.set_quantity('stage', 0)       
2737
2738        # Boundary conditions
2739        Br = Reflective_boundary(domain2)       
2740        Bf = Spatio_temporal_boundary(domain1.filename + '.' + domain1.format,
2741                                      domain2)
2742        domain2.set_boundary({'right':Br, 'bottom':Br, 'diagonal':Bf})
2743        domain2.check_integrity()       
2744       
2745
2746
2747        #Evolution (small steps)
2748        for t in domain2.evolve(yieldstep = 0.0711, finaltime = finaltime):
2749            pass
2750
2751       
2752        #Use output from domain1 as spatio-temporal boundary for domain2
2753        #and verify that results at right hand side are close. 
2754
2755        cv2 = domain2.quantities['stage'].centroid_values
2756
2757        #print take(cv1, (12,14,16))  #Right
2758        #print take(cv2, (4,6,8))
2759        #print take(cv1, (0,6,12))  #Bottom     
2760        #print take(cv2, (0,1,4))
2761        #print take(cv1, (0,8,16))  #Diag       
2762        #print take(cv2, (0,3,8))
2763
2764        assert allclose( take(cv1, (0,8,16)), take(cv2, (0,3,8))) #Diag
2765        assert allclose( take(cv1, (0,6,12)), take(cv2, (0,1,4))) #Bottom
2766        assert allclose( take(cv1, (12,14,16)), take(cv2, (4,6,8))) #RHS
2767
2768        #Cleanup
2769        os.remove(domain1.filename + '.' + domain1.format)       
2770       
2771       
2772
2773    def test_spatio_temporal_boundary_2(self):
2774        """Test that boundary values can be read from file and interpolated
2775        in both time and space.
2776        This is a more basic test, verifying that boundary object
2777        produces the expected results
2778
2779
2780        """
2781        import time
2782       
2783        #Create sww file of simple propagation from left to right
2784        #through rectangular domain
2785       
2786        from mesh_factory import rectangular
2787
2788        #Create basic mesh
2789        points, vertices, boundary = rectangular(3, 3)
2790
2791        #Create shallow water domain
2792        domain1 = Domain(points, vertices, boundary)
2793       
2794        from util import mean
2795        domain1.reduction = mean
2796        domain1.smooth = True #To mimic MOST output
2797
2798        domain1.default_order = 2
2799        domain1.store = True   
2800        domain1.set_datadir('.')
2801        domain1.set_name('spatio_temporal_boundary_source' + str(time.time()))
2802       
2803        #FIXME: This is extremely important!
2804        #How can we test if they weren't stored?
2805        domain1.quantities_to_be_stored = ['stage', 'xmomentum', 'ymomentum']
2806
2807               
2808        #Bed-slope and friction at vertices (and interpolated elsewhere)
2809        domain1.set_quantity('elevation', 0)
2810        domain1.set_quantity('friction', 0)     
2811
2812        # Boundary conditions
2813        Br = Reflective_boundary(domain1)
2814        Bd = Dirichlet_boundary([0.3,0,0])     
2815        domain1.set_boundary({'left': Bd, 'top': Bd, 'right': Br, 'bottom': Br})
2816        #Initial condition
2817        domain1.set_quantity('stage', 0)
2818        domain1.check_integrity()
2819
2820        finaltime = 5
2821        #Evolution  (full domain - large steps)
2822        for t in domain1.evolve(yieldstep = 1, finaltime = finaltime):
2823            pass
2824            #domain1.write_time()
2825
2826       
2827        #Create an triangle shaped domain (coinciding with some
2828        #coordinates from domain 1),
2829        #formed from the lower and right hand  boundaries and
2830        #the sw-ne diagonal
2831        #from domain 1. Call it domain2
2832       
2833        points = [ [0,0], [1.0/3,0], [1.0/3,1.0/3],
2834                   [2.0/3,0], [2.0/3,1.0/3], [2.0/3,2.0/3],
2835                   [1,0], [1,1.0/3], [1,2.0/3], [1,1]]
2836       
2837        vertices = [ [1,2,0], 
2838                     [3,4,1], [2,1,4], [4,5,2], 
2839                     [6,7,3], [4,3,7], [7,8,4], [5,4,8], [8,9,5]]       
2840
2841        boundary = { (0,1):'bottom', (1,1):'bottom', (4,1): 'bottom',
2842                     (4,2):'right', (6,2):'right', (8,2):'right',
2843                     (0,0):'diagonal', (3,0):'diagonal', (8,0):'diagonal'}   
2844                     
2845        domain2 = Domain(points, vertices, boundary)
2846       
2847        domain2.reduction = domain1.reduction
2848        domain2.smooth = False
2849        domain2.default_order = 2
2850       
2851        #Bed-slope and friction at vertices (and interpolated elsewhere)
2852        domain2.set_quantity('elevation', 0)
2853        domain2.set_quantity('friction', 0)
2854        domain2.set_quantity('stage', 0)       
2855
2856       
2857        #Read results for specific timesteps t=1 and t=2
2858        from Scientific.IO.NetCDF import NetCDFFile     
2859        fid = NetCDFFile(domain1.filename + '.' + domain1.format)
2860       
2861        x = fid.variables['x'][:]
2862        y = fid.variables['y'][:]       
2863        s1 = fid.variables['stage'][1,:]
2864        s2 = fid.variables['stage'][2,:]       
2865        fid.close()
2866
2867        from Numeric import take, reshape, concatenate 
2868        shp = (len(x), 1)
2869        points = concatenate( (reshape(x, shp), reshape(y, shp)), axis=1)       
2870        #The diagonal points of domain 1 are 0, 5, 10, 15
2871       
2872        #print points[0], points[5], points[10], points[15]
2873        assert allclose( take(points, [0,5,10,15]), 
2874                         [[0,0], [1.0/3, 1.0/3], [2.0/3, 2.0/3], [1,1]]) 
2875
2876       
2877        # Boundary conditions
2878        Br = Reflective_boundary(domain2)       
2879        Bf = Spatio_temporal_boundary(domain1.filename + '.' + domain1.format,
2880                                      domain2)
2881        domain2.set_boundary({'right':Br, 'bottom':Br, 'diagonal':Bf})
2882        domain2.check_integrity()       
2883
2884        #Test that interpolation points are the mid points of the all boundary
2885        #segments       
2886
2887        boundary_midpoints = [[1.0/6, 0], [1.0/2, 0], [5.0/6,0],
2888                              [1.0, 1.0/6], [1.0, 1.0/2], [1.0, 5.0/6],
2889                              [1.0/6, 1.0/6], [0.5, 0.5], [5.0/6, 5.0/6]]
2890                             
2891        boundary_midpoints.sort()
2892        R = Bf.F.interpolation_points.tolist() 
2893        R.sort()
2894        assert allclose(boundary_midpoints, R) 
2895       
2896        #Check spatially interpolated output at time == 1
2897        domain2.time = 1
2898       
2899        #First diagonal midpoint
2900        R0 = Bf.evaluate(0,0)
2901        assert allclose(R0[0], (s1[0] + s1[5])/2)
2902       
2903        #Second diagonal midpoint
2904        R0 = Bf.evaluate(3,0)
2905        assert allclose(R0[0], (s1[5] + s1[10])/2)
2906       
2907        #First diagonal midpoint
2908        R0 = Bf.evaluate(8,0)
2909        assert allclose(R0[0], (s1[10] + s1[15])/2)             
2910               
2911        #Check spatially interpolated output at time == 2
2912        domain2.time = 2
2913       
2914        #First diagonal midpoint
2915        R0 = Bf.evaluate(0,0)
2916        assert allclose(R0[0], (s2[0] + s2[5])/2)
2917       
2918        #Second diagonal midpoint
2919        R0 = Bf.evaluate(3,0)
2920        assert allclose(R0[0], (s2[5] + s2[10])/2)
2921       
2922        #First diagonal midpoint
2923        R0 = Bf.evaluate(8,0)
2924        assert allclose(R0[0], (s2[10] + s2[15])/2)             
2925       
2926       
2927        #Now check temporal interpolation
2928
2929        domain2.time = 1 + 2.0/3
2930       
2931        #First diagonal midpoint
2932        R0 = Bf.evaluate(0,0)
2933        assert allclose(R0[0], ((s1[0] + s1[5])/2 + 2.0*(s2[0] + s2[5])/2)/3)
2934       
2935        #Second diagonal midpoint
2936        R0 = Bf.evaluate(3,0)
2937        assert allclose(R0[0], ((s1[5] + s1[10])/2 + 2.0*(s2[5] + s2[10])/2)/3)
2938       
2939        #First diagonal midpoint
2940        R0 = Bf.evaluate(8,0)
2941        assert allclose(R0[0], ((s1[10] + s1[15])/2 + 2.0*(s2[10] + s2[15])/2)/3)               
2942               
2943       
2944               
2945        #Cleanup
2946        os.remove(domain1.filename + '.' + domain1.format)                 
2947#-------------------------------------------------------------
2948if __name__ == "__main__":
2949    suite = unittest.makeSuite(TestCase,'test')
2950    #suite = unittest.makeSuite(TestCase,'test_balance_deep_and_shallow')
2951    runner = unittest.TextTestRunner()
2952    runner.run(suite)
2953
Note: See TracBrowser for help on using the repository browser.