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

Last change on this file since 1014 was 1014, checked in by ole, 19 years ago
File size: 103.9 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=2
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=2
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  #FIXME: 2 will fail
1646        domain.beta_h = 0
1647
1648        #IC
1649        def x_slope(x, y):
1650            z = 0*x
1651            for i in range(len(x)):
1652                if x[i] < 0.3: 
1653                    z[i] = x[i]/3
1654                if 0.3 <= x[i] < 0.5: 
1655                    z[i] = -0.5
1656                if 0.5 <= x[i] < 0.7: 
1657                    z[i] = 0.4       
1658                    #z[i] = 0.7                     
1659                if 0.7 <= x[i]: 
1660                    z[i] = x[i]/3                   
1661            return z                               
1662                                   
1663                   
1664
1665        domain.set_quantity('elevation', x_slope)
1666        domain.set_quantity('friction', 0)     
1667        domain.set_quantity('stage', 0.4) #Steady       
1668       
1669        # Boundary conditions (reflective everywhere)
1670        Br = Reflective_boundary(domain)
1671        domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
1672
1673        domain.check_integrity()
1674
1675        #domain.visualise = True #If you want to take a sticky beak 
1676       
1677        initial_volume = domain.quantities['stage'].get_integral()         
1678        initial_xmom = domain.quantities['xmomentum'].get_integral()     
1679       
1680       
1681        #Let surface settle according to limiter
1682        #FIXME: What exactly is the significance of this?
1683        #Is it a problem that an arbitrary initial condition
1684        #needs to respect the law of the limiters?
1685        domain.distribute_to_vertices_and_edges()
1686        initial_volume = domain.quantities['stage'].get_integral()         
1687               
1688        #print initial_volume
1689       
1690        #print initial_xmom
1691
1692        #Evolution
1693        for t in domain.evolve(yieldstep = 0.05, finaltime = 10.0):
1694            volume =  domain.quantities['stage'].get_integral()         
1695            #print t, volume, initial_volume
1696           
1697            assert allclose (volume, initial_volume) 
1698         
1699            #FIXME: What would we expect from momentum
1700            #xmom = domain.quantities['xmomentum'].get_integral()   
1701            #print xmom
1702            #assert allclose (xmom, initial_xmom)
1703       
1704       
1705               
1706    def test_second_order_flat_bed_onestep(self):
1707
1708        from mesh_factory import rectangular
1709        from shallow_water import Domain, Reflective_boundary,\
1710             Dirichlet_boundary, Constant_height
1711        from Numeric import array
1712
1713        #Create basic mesh
1714        points, vertices, boundary = rectangular(6, 6)
1715
1716        #Create shallow water domain
1717        domain = Domain(points, vertices, boundary)
1718        domain.smooth = False
1719        domain.default_order=2
1720
1721        # Boundary conditions
1722        Br = Reflective_boundary(domain)
1723        Bd = Dirichlet_boundary([0.1, 0., 0.])
1724        domain.set_boundary({'left': Bd, 'right': Br, 'top': Br, 'bottom': Br})
1725
1726        domain.check_integrity()
1727
1728        #Evolution
1729        for t in domain.evolve(yieldstep = 0.05, finaltime = 0.05):
1730            pass# domain.write_time()
1731
1732        #Data from earlier version of pyvolution
1733        assert allclose(domain.min_timestep, 0.0396825396825)
1734        assert allclose(domain.max_timestep, 0.0396825396825) 
1735
1736        assert allclose(domain.quantities['stage'].centroid_values[:12],
1737                        [0.00171396, 0.02656103, 0.00241523, 0.02656103,
1738                        0.00241523, 0.02656103, 0.00241523, 0.02656103,
1739                        0.00241523, 0.02656103, 0.00241523, 0.0272623])
1740
1741        domain.distribute_to_vertices_and_edges()
1742        assert allclose(domain.quantities['stage'].vertex_values[:12,0],
1743                        [0.0001714, 0.02656103, 0.00024152,
1744                        0.02656103, 0.00024152, 0.02656103,
1745                        0.00024152, 0.02656103, 0.00024152,
1746                        0.02656103, 0.00024152, 0.0272623])
1747
1748        assert allclose(domain.quantities['stage'].vertex_values[:12,1],
1749                        [0.00315012, 0.02656103, 0.00024152, 0.02656103,
1750                         0.00024152, 0.02656103, 0.00024152, 0.02656103,
1751                         0.00024152, 0.02656103, 0.00040506, 0.0272623])
1752       
1753        assert allclose(domain.quantities['stage'].vertex_values[:12,2],
1754                        [0.00182037, 0.02656103, 0.00676264,
1755                         0.02656103, 0.00676264, 0.02656103,
1756                         0.00676264, 0.02656103, 0.00676264,
1757                         0.02656103, 0.0065991, 0.0272623])
1758
1759        assert allclose(domain.quantities['xmomentum'].centroid_values[:12],
1760                        [0.00113961, 0.01302432, 0.00148672,
1761                         0.01302432, 0.00148672, 0.01302432,
1762                         0.00148672, 0.01302432, 0.00148672 ,
1763                         0.01302432, 0.00148672, 0.01337143])
1764
1765        assert allclose(domain.quantities['ymomentum'].centroid_values[:12],
1766                        [-2.91240050e-004 , 1.22721531e-004,
1767                         -1.22721531e-004, 1.22721531e-004 ,
1768                         -1.22721531e-004, 1.22721531e-004,
1769                         -1.22721531e-004 , 1.22721531e-004,
1770                         -1.22721531e-004, 1.22721531e-004,
1771                         -1.22721531e-004, -4.57969873e-005])
1772
1773
1774       
1775    def test_second_order_flat_bed_moresteps(self):
1776
1777        from mesh_factory import rectangular
1778        from shallow_water import Domain, Reflective_boundary,\
1779             Dirichlet_boundary, Constant_height
1780        from Numeric import array
1781
1782        #Create basic mesh
1783        points, vertices, boundary = rectangular(6, 6)
1784
1785        #Create shallow water domain
1786        domain = Domain(points, vertices, boundary)
1787        domain.smooth = False
1788        domain.default_order=2
1789
1790        # Boundary conditions
1791        Br = Reflective_boundary(domain)
1792        Bd = Dirichlet_boundary([0.1, 0., 0.])
1793        domain.set_boundary({'left': Bd, 'right': Br, 'top': Br, 'bottom': Br})
1794
1795        domain.check_integrity()
1796
1797        #Evolution
1798        for t in domain.evolve(yieldstep = 0.05, finaltime = 0.1):
1799            pass
1800
1801        #Data from earlier version of pyvolution
1802        #assert allclose(domain.min_timestep, 0.0396825396825)
1803        #assert allclose(domain.max_timestep, 0.0396825396825)
1804        #print domain.quantities['stage'].centroid_values
1805
1806
1807    def test_flatbed_first_order(self):
1808        from mesh_factory import rectangular
1809        from shallow_water import Domain,\
1810             Reflective_boundary, Dirichlet_boundary
1811
1812        from Numeric import array
1813
1814        #Create basic mesh
1815        N = 8
1816        points, vertices, boundary = rectangular(N, N)
1817
1818        #Create shallow water domain
1819        domain = Domain(points, vertices, boundary)
1820        domain.smooth = False
1821        domain.visualise = False
1822        domain.default_order=1
1823
1824        # Boundary conditions
1825        Br = Reflective_boundary(domain)
1826        Bd = Dirichlet_boundary([0.2,0.,0.])
1827       
1828        domain.set_boundary({'left': Bd, 'right': Br, 'top': Br, 'bottom': Br})
1829        domain.check_integrity()
1830
1831 
1832        #Evolution
1833        for t in domain.evolve(yieldstep = 0.02, finaltime = 0.5):
1834            pass
1835            #domain.write_time()
1836
1837        #FIXME: These numbers were from version before 25/10       
1838        #assert allclose(domain.min_timestep, 0.0140413643926)
1839        #assert allclose(domain.max_timestep, 0.0140947355753)
1840
1841        for i in range(3):
1842            #assert allclose(domain.quantities['stage'].edge_values[:4,i],
1843            #                [0.10730244,0.12337617,0.11200126,0.12605666])
1844
1845            assert allclose(domain.quantities['xmomentum'].edge_values[:4,i],
1846                            [0.07610894,0.06901572,0.07284461,0.06819712])
1847
1848            #assert allclose(domain.quantities['ymomentum'].edge_values[:4,i],
1849            #                [-0.0060238, -0.00157404, -0.00309633, -0.0001637])                         
1850       
1851    def test_flatbed_second_order(self):
1852        from mesh_factory import rectangular
1853        from shallow_water import Domain,\
1854             Reflective_boundary, Dirichlet_boundary
1855
1856        from Numeric import array
1857
1858        #Create basic mesh
1859        N = 8
1860        points, vertices, boundary = rectangular(N, N)
1861
1862        #Create shallow water domain
1863        domain = Domain(points, vertices, boundary)
1864        domain.smooth = False
1865        domain.visualise = False
1866        domain.default_order=2
1867        #domain.minimum_allowed_height = 0.0 #Makes it like the 'oldstyle' balance
1868
1869        # Boundary conditions
1870        Br = Reflective_boundary(domain)
1871        Bd = Dirichlet_boundary([0.2,0.,0.])
1872
1873        domain.set_boundary({'left': Bd, 'right': Br, 'top': Br, 'bottom': Br})
1874        domain.check_integrity()
1875 
1876        #Evolution
1877        for t in domain.evolve(yieldstep = 0.01, finaltime = 0.03):
1878            pass
1879
1880
1881        assert allclose(domain.min_timestep, 0.0210448446782)
1882        assert allclose(domain.max_timestep, 0.0210448446782) 
1883
1884        #print domain.quantities['stage'].vertex_values[:4,0]
1885        #print domain.quantities['xmomentum'].vertex_values[:4,0]
1886        #print domain.quantities['ymomentum'].vertex_values[:4,0]
1887
1888        #FIXME: These numbers were from version before 25/10
1889        #assert allclose(domain.quantities['stage'].vertex_values[:4,0],
1890        #                [0.00101913,0.05352143,0.00104852,0.05354394])
1891       
1892        assert allclose(domain.quantities['xmomentum'].vertex_values[:4,0],
1893                        [ 0.00064835, 0.03685719, 0.00085073, 0.03687313])
1894
1895        #assert allclose(domain.quantities['xmomentum'].vertex_values[:4,0],
1896        #                [0.00090581,0.03685719,0.00088303,0.03687313])
1897       
1898        assert allclose(domain.quantities['ymomentum'].vertex_values[:4,0],
1899                        [-0.00139463,0.0006156,-0.00060364,0.00061827])
1900
1901
1902       
1903    def test_flatbed_second_order_distribute(self):
1904        #Use real data from pyvolution 2
1905        #painfully setup and extracted.
1906        from mesh_factory import rectangular
1907        from shallow_water import Domain,\
1908             Reflective_boundary, Dirichlet_boundary
1909
1910        from Numeric import array
1911
1912        #Create basic mesh
1913        N = 8
1914        points, vertices, boundary = rectangular(N, N)
1915
1916        #Create shallow water domain
1917        domain = Domain(points, vertices, boundary)
1918        domain.smooth = False
1919        domain.visualise = False
1920        domain.default_order=domain.order=2
1921
1922        # Boundary conditions
1923        Br = Reflective_boundary(domain)
1924        Bd = Dirichlet_boundary([0.2,0.,0.])
1925
1926        domain.set_boundary({'left': Bd, 'right': Br, 'top': Br, 'bottom': Br})
1927        domain.check_integrity()
1928
1929
1930
1931        for V in [False, True]:
1932            if V:
1933                #Set centroids as if system had been evolved           
1934                L = zeros(2*N*N, Float)
1935                L[:32] = [7.21205592e-003, 5.35214298e-002, 1.00910824e-002,
1936                          5.35439433e-002, 1.00910824e-002, 5.35439433e-002,
1937                          1.00910824e-002, 5.35439433e-002, 1.00910824e-002,
1938                          5.35439433e-002, 1.00910824e-002, 5.35439433e-002,
1939                          1.00910824e-002, 5.35393928e-002, 1.02344264e-002,
1940                          5.59605058e-002, 0.00000000e+000, 3.31027800e-004,
1941                          0.00000000e+000, 4.37962142e-005, 0.00000000e+000,
1942                          4.37962142e-005, 0.00000000e+000, 4.37962142e-005,
1943                          0.00000000e+000, 4.37962142e-005, 0.00000000e+000,
1944                          4.37962142e-005, 0.00000000e+000, 4.37962142e-005,
1945                          0.00000000e+000, 5.57305948e-005]
1946       
1947                X = zeros(2*N*N, Float)
1948                X[:32] = [6.48351607e-003, 3.68571894e-002, 8.50733285e-003,
1949                          3.68731327e-002, 8.50733285e-003, 3.68731327e-002,
1950                          8.50733285e-003, 3.68731327e-002, 8.50733285e-003,
1951                          3.68731327e-002, 8.50733285e-003, 3.68731327e-002,
1952                          8.50733285e-003, 3.68693861e-002, 8.65220973e-003,
1953                          3.85055387e-002, 0.00000000e+000, 2.86060840e-004,
1954                          0.00000000e+000, 3.58905503e-005, 0.00000000e+000,
1955                          3.58905503e-005, 0.00000000e+000, 3.58905503e-005,
1956                          0.00000000e+000, 3.58905503e-005, 0.00000000e+000,
1957                          3.58905503e-005, 0.00000000e+000, 3.58905503e-005,
1958                          0.00000000e+000, 4.57662812e-005]
1959
1960                Y = zeros(2*N*N, Float)
1961                Y[:32]=[-1.39463104e-003, 6.15600298e-004, -6.03637382e-004,
1962                        6.18272251e-004, -6.03637382e-004, 6.18272251e-004,
1963                        -6.03637382e-004, 6.18272251e-004, -6.03637382e-004,
1964                        6.18272251e-004, -6.03637382e-004, 6.18272251e-004,
1965                        -6.03637382e-004, 6.18599320e-004, -6.74622797e-004,
1966                        -1.48934756e-004, 0.00000000e+000, -5.35079969e-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.57264987e-005, 0.00000000e+000,
1970                        -2.57264987e-005, 0.00000000e+000, -2.57264987e-005,
1971                        0.00000000e+000, -2.57635178e-005]
1972
1973       
1974                domain.set_quantity('stage', L, 'centroids')
1975                domain.set_quantity('xmomentum', X, 'centroids')
1976                domain.set_quantity('ymomentum', Y, 'centroids')
1977       
1978                domain.check_integrity()
1979            else:   
1980                #Evolution
1981                for t in domain.evolve(yieldstep = 0.01, finaltime = 0.03):
1982                    pass
1983                assert allclose(domain.min_timestep, 0.0210448446782)
1984                assert allclose(domain.max_timestep, 0.0210448446782)
1985       
1986
1987            #Centroids were correct but not vertices.
1988            #Hence the check of distribute below.
1989            assert allclose(domain.quantities['stage'].centroid_values[:4],
1990                            [0.00721206,0.05352143,0.01009108,0.05354394])
1991
1992            assert allclose(domain.quantities['xmomentum'].centroid_values[:4],
1993                            [0.00648352,0.03685719,0.00850733,0.03687313])
1994
1995            assert allclose(domain.quantities['ymomentum'].centroid_values[:4],
1996                            [-0.00139463,0.0006156,-0.00060364,0.00061827])
1997
1998            #print 'C17=', domain.quantities['xmomentum'].centroid_values[17]
1999            #print 'C19=', domain.quantities['xmomentum'].centroid_values[19]
2000
2001            #assert allclose(domain.quantities['xmomentum'].centroid_values[17],0.00028606084)
2002            ##print domain.quantities['xmomentum'].centroid_values[17], V
2003            ##print
2004            if not V:           
2005                assert allclose(domain.quantities['xmomentum'].centroid_values[17], 0.0)
2006            else: 
2007                assert allclose(domain.quantities['xmomentum'].centroid_values[17], 0.00028606084)           
2008
2009            import copy
2010            XX = copy.copy(domain.quantities['xmomentum'].centroid_values)
2011            assert allclose(domain.quantities['xmomentum'].centroid_values, XX)
2012       
2013            domain.distribute_to_vertices_and_edges()
2014
2015            #assert allclose(domain.quantities['xmomentum'].centroid_values, XX)
2016
2017            assert allclose(domain.quantities['xmomentum'].centroid_values[17],
2018                            0.0)
2019
2020
2021            #FIXME: These numbers were from version before 25/10
2022            #assert allclose(domain.quantities['stage'].vertex_values[:4,0],
2023            #                [0.00101913,0.05352143,0.00104852,0.05354394])
2024
2025            assert allclose(domain.quantities['ymomentum'].vertex_values[:4,0],
2026                            [-0.00139463,0.0006156,-0.00060364,0.00061827])
2027
2028
2029            assert allclose(domain.quantities['xmomentum'].vertex_values[:4,0],
2030                            [0.00064835,0.03685719,0.00085073,0.03687313])
2031
2032
2033            #NB NO longer relvant:
2034
2035            #This was the culprit. First triangles vertex 0 had an
2036            #x-momentum of 0.0064835 instead of 0.00090581 and
2037            #third triangle had 0.00850733 instead of 0.00088303
2038            #print domain.quantities['xmomentum'].vertex_values[:4,0]
2039
2040            #print domain.quantities['xmomentum'].vertex_values[:4,0]
2041            #assert allclose(domain.quantities['xmomentum'].vertex_values[:4,0],
2042            #                [0.00090581,0.03685719,0.00088303,0.03687313])       
2043
2044
2045
2046       
2047       
2048    def test_bedslope_problem_first_order(self):
2049
2050        from mesh_factory import rectangular
2051        from shallow_water import Domain, Reflective_boundary, Constant_height
2052        from Numeric import array
2053
2054        #Create basic mesh
2055        points, vertices, boundary = rectangular(6, 6)
2056
2057        #Create shallow water domain
2058        domain = Domain(points, vertices, boundary)
2059        domain.smooth = False
2060        domain.default_order=1
2061
2062        #Bed-slope and friction
2063        def x_slope(x, y):
2064            return -x/3
2065
2066        domain.set_quantity('elevation', x_slope)
2067
2068        # Boundary conditions
2069        Br = Reflective_boundary(domain)
2070        domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
2071
2072        #Initial condition
2073        domain.set_quantity('stage', Constant_height(x_slope, 0.05))
2074        domain.check_integrity()
2075
2076        #Evolution
2077        for t in domain.evolve(yieldstep = 0.05, finaltime = 0.05):
2078            pass# domain.write_time()
2079
2080        assert allclose(domain.min_timestep, 0.050010003001)
2081        assert allclose(domain.max_timestep, 0.050010003001)
2082
2083       
2084    def test_bedslope_problem_first_order_moresteps(self):
2085
2086        from mesh_factory import rectangular
2087        from shallow_water import Domain, Reflective_boundary, Constant_height
2088        from Numeric import array
2089
2090        #Create basic mesh
2091        points, vertices, boundary = rectangular(6, 6)
2092
2093        #Create shallow water domain
2094        domain = Domain(points, vertices, boundary)
2095        domain.smooth = False
2096        domain.default_order=1
2097        domain.beta_h = 0.0 #Use first order in h-limiter
2098
2099        #Bed-slope and friction
2100        def x_slope(x, y):
2101            return -x/3
2102
2103        domain.set_quantity('elevation', x_slope)
2104
2105        # Boundary conditions
2106        Br = Reflective_boundary(domain)
2107        domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
2108
2109        #Initial condition
2110        domain.set_quantity('stage', Constant_height(x_slope, 0.05))
2111        domain.check_integrity()
2112
2113        #Evolution
2114        for t in domain.evolve(yieldstep = 0.05, finaltime = 0.5):
2115            pass# domain.write_time()
2116
2117        #Data from earlier version of pyvolution
2118        #print domain.quantities['stage'].centroid_values
2119       
2120        assert allclose(domain.quantities['stage'].centroid_values,
2121                        [-0.02998628, -0.01520652, -0.03043492,
2122                        -0.0149132, -0.03004706, -0.01476251,
2123                        -0.0298215, -0.01467976, -0.02988158,
2124                        -0.01474662, -0.03036161, -0.01442995,
2125                        -0.07624583, -0.06297061, -0.07733792,
2126                        -0.06342237, -0.07695439, -0.06289595,
2127                        -0.07635559, -0.0626065, -0.07633628,
2128                        -0.06280072, -0.07739632, -0.06386738,
2129                        -0.12161738, -0.11028239, -0.1223796,
2130                        -0.11095953, -0.12189744, -0.11048616,
2131                        -0.12074535, -0.10987605, -0.12014311,
2132                        -0.10976691, -0.12096859, -0.11087692,
2133                        -0.16868259, -0.15868061, -0.16801135,
2134                        -0.1588003, -0.16674343, -0.15813323,
2135                        -0.16457595, -0.15693826, -0.16281096,
2136                        -0.15585154, -0.16283873, -0.15540068,
2137                        -0.17450362, -0.19919913, -0.18062882,
2138                        -0.19764131, -0.17783111, -0.19407213,
2139                        -0.1736915, -0.19053624, -0.17228678,
2140                        -0.19105634, -0.17920133, -0.1968828,
2141                        -0.14244395, -0.14604641, -0.14473537,
2142                        -0.1506107, -0.14510055, -0.14919522,
2143                        -0.14175896, -0.14560798, -0.13911658,
2144                        -0.14439383, -0.13924047, -0.14829043])
2145                       
2146       
2147    def test_bedslope_problem_second_order_one_step(self):
2148
2149        from mesh_factory import rectangular
2150        from shallow_water import Domain, Reflective_boundary, Constant_height
2151        from Numeric import array
2152
2153        #Create basic mesh
2154        points, vertices, boundary = rectangular(6, 6)
2155
2156        #Create shallow water domain
2157        domain = Domain(points, vertices, boundary)
2158        domain.smooth = False
2159        domain.default_order=2
2160
2161        #Bed-slope and friction at vertices (and interpolated elsewhere)
2162        def x_slope(x, y):
2163            return -x/3
2164
2165        domain.set_quantity('elevation', x_slope)
2166
2167        # Boundary conditions
2168        Br = Reflective_boundary(domain)
2169        domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
2170
2171        #Initial condition
2172        domain.set_quantity('stage', Constant_height(x_slope, 0.05))
2173        domain.check_integrity()
2174
2175        assert allclose(domain.quantities['stage'].centroid_values,
2176                        [0.01296296, 0.03148148, 0.01296296,
2177                        0.03148148, 0.01296296, 0.03148148,
2178                        0.01296296, 0.03148148, 0.01296296,
2179                        0.03148148, 0.01296296, 0.03148148,
2180                        -0.04259259, -0.02407407, -0.04259259,
2181                        -0.02407407, -0.04259259, -0.02407407,
2182                        -0.04259259, -0.02407407, -0.04259259,
2183                        -0.02407407, -0.04259259, -0.02407407,
2184                        -0.09814815, -0.07962963, -0.09814815,
2185                        -0.07962963, -0.09814815, -0.07962963,
2186                        -0.09814815, -0.07962963, -0.09814815,
2187                        -0.07962963, -0.09814815, -0.07962963,
2188                        -0.1537037 , -0.13518519, -0.1537037,
2189                        -0.13518519, -0.1537037, -0.13518519,
2190                        -0.1537037 , -0.13518519, -0.1537037,
2191                        -0.13518519, -0.1537037, -0.13518519,
2192                        -0.20925926, -0.19074074, -0.20925926,
2193                        -0.19074074, -0.20925926, -0.19074074,
2194                        -0.20925926, -0.19074074, -0.20925926,
2195                        -0.19074074, -0.20925926, -0.19074074,
2196                        -0.26481481, -0.2462963, -0.26481481,
2197                        -0.2462963, -0.26481481, -0.2462963,
2198                        -0.26481481, -0.2462963, -0.26481481,
2199                        -0.2462963, -0.26481481, -0.2462963])
2200
2201
2202        #print domain.quantities['stage'].extrapolate_second_order()
2203        #domain.distribute_to_vertices_and_edges()
2204        #print domain.quantities['stage'].vertex_values[:,0]       
2205       
2206        #Evolution
2207        for t in domain.evolve(yieldstep = 0.05, finaltime = 0.05):
2208            #domain.write_time()           
2209            pass
2210
2211
2212        #print domain.quantities['stage'].centroid_values
2213        assert allclose(domain.quantities['stage'].centroid_values,
2214                 [0.01290985, 0.02356019, 0.01619096, 0.02356019, 0.01619096,
2215                  0.02356019, 0.01619096, 0.02356019, 0.01619096, 0.02356019,
2216                  0.01619096, 0.0268413, -0.04411074, -0.0248011, -0.04186556,
2217                  -0.0248011, -0.04186556, -0.0248011, -0.04186556, -0.0248011,
2218                  -0.04186556, -0.0248011, -0.04186556, -0.02255593,
2219                  -0.09966629, -0.08035666, -0.09742112, -0.08035666,
2220                  -0.09742112, -0.08035666, -0.09742112, -0.08035666,
2221                  -0.09742112, -0.08035666, -0.09742112, -0.07811149,
2222                  -0.15522185, -0.13591222, -0.15297667, -0.13591222,
2223                  -0.15297667, -0.13591222, -0.15297667, -0.13591222,
2224                  -0.15297667, -0.13591222, -0.15297667, -0.13366704,
2225                  -0.2107774, -0.19146777, -0.20853223, -0.19146777,
2226                  -0.20853223, -0.19146777, -0.20853223, -0.19146777,
2227                  -0.20853223, -0.19146777, -0.20853223, -0.1892226,
2228                  -0.26120669, -0.24776246, -0.25865535, -0.24776246,
2229                  -0.25865535, -0.24776246, -0.25865535, -0.24776246,
2230                  -0.25865535, -0.24776246, -0.25865535, -0.24521113])
2231
2232       
2233
2234    def test_bedslope_problem_second_order_two_steps(self):
2235
2236        from mesh_factory import rectangular
2237        from shallow_water import Domain, Reflective_boundary, Constant_height
2238        from Numeric import array
2239
2240        #Create basic mesh
2241        points, vertices, boundary = rectangular(6, 6)
2242
2243        #Create shallow water domain
2244        domain = Domain(points, vertices, boundary)
2245        domain.smooth = False
2246        domain.default_order=2
2247        domain.beta_h = 0.0 #Use first order in h-limiter       
2248
2249        #Bed-slope and friction at vertices (and interpolated elsewhere)
2250        def x_slope(x, y):
2251            return -x/3
2252
2253        domain.set_quantity('elevation', x_slope)
2254
2255        # Boundary conditions
2256        Br = Reflective_boundary(domain)
2257        domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
2258
2259        #Initial condition
2260        domain.set_quantity('stage', Constant_height(x_slope, 0.05))
2261        domain.check_integrity()
2262
2263        assert allclose(domain.quantities['stage'].centroid_values,
2264                        [0.01296296, 0.03148148, 0.01296296,
2265                        0.03148148, 0.01296296, 0.03148148,
2266                        0.01296296, 0.03148148, 0.01296296,
2267                        0.03148148, 0.01296296, 0.03148148,
2268                        -0.04259259, -0.02407407, -0.04259259,
2269                        -0.02407407, -0.04259259, -0.02407407,
2270                        -0.04259259, -0.02407407, -0.04259259,
2271                        -0.02407407, -0.04259259, -0.02407407,
2272                        -0.09814815, -0.07962963, -0.09814815,
2273                        -0.07962963, -0.09814815, -0.07962963,
2274                        -0.09814815, -0.07962963, -0.09814815,
2275                        -0.07962963, -0.09814815, -0.07962963,
2276                        -0.1537037 , -0.13518519, -0.1537037,
2277                        -0.13518519, -0.1537037, -0.13518519,
2278                        -0.1537037 , -0.13518519, -0.1537037,
2279                        -0.13518519, -0.1537037, -0.13518519,
2280                        -0.20925926, -0.19074074, -0.20925926,
2281                        -0.19074074, -0.20925926, -0.19074074,
2282                        -0.20925926, -0.19074074, -0.20925926,
2283                        -0.19074074, -0.20925926, -0.19074074,
2284                        -0.26481481, -0.2462963, -0.26481481,
2285                        -0.2462963, -0.26481481, -0.2462963,
2286                        -0.26481481, -0.2462963, -0.26481481,
2287                        -0.2462963, -0.26481481, -0.2462963])
2288
2289
2290        #print domain.quantities['stage'].extrapolate_second_order()
2291        #domain.distribute_to_vertices_and_edges()
2292        #print domain.quantities['stage'].vertex_values[:,0]       
2293       
2294        #Evolution
2295        for t in domain.evolve(yieldstep = 0.05, finaltime = 0.1):
2296            pass
2297
2298       
2299        #Data from earlier version of pyvolution ft=0.1
2300        assert allclose(domain.min_timestep, 0.0376895634803) 
2301        assert allclose(domain.max_timestep, 0.0415635655309)
2302
2303
2304        assert allclose(domain.quantities['stage'].centroid_values,
2305                        [0.00855788, 0.01575204, 0.00994606, 0.01717072,
2306                         0.01005985, 0.01716362, 0.01005985, 0.01716299,
2307                         0.01007098, 0.01736248, 0.01216452, 0.02026776,
2308                         -0.04462374, -0.02479045, -0.04199789, -0.0229465,
2309                         -0.04184033, -0.02295693, -0.04184013, -0.02295675,
2310                         -0.04184486, -0.0228168, -0.04028876, -0.02036486,
2311                         -0.10029444, -0.08170809, -0.09772846, -0.08021704,
2312                         -0.09760006, -0.08022143, -0.09759984, -0.08022124,
2313                         -0.09760261, -0.08008893, -0.09603914, -0.07758209,
2314                         -0.15584152, -0.13723138, -0.15327266, -0.13572906,
2315                         -0.15314427, -0.13573349, -0.15314405, -0.13573331,
2316                         -0.15314679, -0.13560104, -0.15158523, -0.13310701,
2317                         -0.21208605, -0.19283913, -0.20955631, -0.19134189,
2318                         -0.20942821, -0.19134598, -0.20942799, -0.1913458,
2319                         -0.20943005, -0.19120952, -0.20781177, -0.18869401,
2320                         -0.25384082, -0.2463294, -0.25047649, -0.24464654,
2321                         -0.25031159, -0.24464253, -0.25031112, -0.24464253,
2322                         -0.25031463, -0.24454764, -0.24885323, -0.24286438])
2323
2324       
2325       
2326               
2327    def test_bedslope_problem_second_order_two_yieldsteps(self):
2328
2329        from mesh_factory import rectangular
2330        from shallow_water import Domain, Reflective_boundary, Constant_height
2331        from Numeric import array
2332
2333        #Create basic mesh
2334        points, vertices, boundary = rectangular(6, 6)
2335
2336        #Create shallow water domain
2337        domain = Domain(points, vertices, boundary)
2338        domain.smooth = False
2339        domain.default_order=2
2340        domain.beta_h = 0.0 #Use first order in h-limiter       
2341
2342        #Bed-slope and friction at vertices (and interpolated elsewhere)
2343        def x_slope(x, y):
2344            return -x/3
2345
2346        domain.set_quantity('elevation', x_slope)
2347
2348        # Boundary conditions
2349        Br = Reflective_boundary(domain)
2350        domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
2351
2352        #Initial condition
2353        domain.set_quantity('stage', Constant_height(x_slope, 0.05))
2354        domain.check_integrity()
2355
2356        assert allclose(domain.quantities['stage'].centroid_values,
2357                        [0.01296296, 0.03148148, 0.01296296,
2358                        0.03148148, 0.01296296, 0.03148148,
2359                        0.01296296, 0.03148148, 0.01296296,
2360                        0.03148148, 0.01296296, 0.03148148,
2361                        -0.04259259, -0.02407407, -0.04259259,
2362                        -0.02407407, -0.04259259, -0.02407407,
2363                        -0.04259259, -0.02407407, -0.04259259,
2364                        -0.02407407, -0.04259259, -0.02407407,
2365                        -0.09814815, -0.07962963, -0.09814815,
2366                        -0.07962963, -0.09814815, -0.07962963,
2367                        -0.09814815, -0.07962963, -0.09814815,
2368                        -0.07962963, -0.09814815, -0.07962963,
2369                        -0.1537037 , -0.13518519, -0.1537037,
2370                        -0.13518519, -0.1537037, -0.13518519,
2371                        -0.1537037 , -0.13518519, -0.1537037,
2372                        -0.13518519, -0.1537037, -0.13518519,
2373                        -0.20925926, -0.19074074, -0.20925926,
2374                        -0.19074074, -0.20925926, -0.19074074,
2375                        -0.20925926, -0.19074074, -0.20925926,
2376                        -0.19074074, -0.20925926, -0.19074074,
2377                        -0.26481481, -0.2462963, -0.26481481,
2378                        -0.2462963, -0.26481481, -0.2462963,
2379                        -0.26481481, -0.2462963, -0.26481481,
2380                        -0.2462963, -0.26481481, -0.2462963])
2381
2382
2383        #print domain.quantities['stage'].extrapolate_second_order()
2384        #domain.distribute_to_vertices_and_edges()
2385        #print domain.quantities['stage'].vertex_values[:,0]       
2386       
2387        #Evolution
2388        for t in domain.evolve(yieldstep = 0.05, finaltime = 0.1):   #0.05??
2389            #domain.write_time()           
2390            pass
2391
2392
2393       
2394        assert allclose(domain.quantities['stage'].centroid_values, 
2395                 [0.00855788, 0.01575204, 0.00994606, 0.01717072, 0.01005985,
2396                  0.01716362, 0.01005985, 0.01716299, 0.01007098, 0.01736248,
2397                  0.01216452, 0.02026776, -0.04462374, -0.02479045, -0.04199789,
2398                  -0.0229465, -0.04184033, -0.02295693, -0.04184013,
2399                  -0.02295675, -0.04184486, -0.0228168, -0.04028876,
2400                  -0.02036486, -0.10029444, -0.08170809, -0.09772846,
2401                  -0.08021704, -0.09760006, -0.08022143, -0.09759984,
2402                  -0.08022124, -0.09760261, -0.08008893, -0.09603914,
2403                  -0.07758209, -0.15584152, -0.13723138, -0.15327266,
2404                  -0.13572906, -0.15314427, -0.13573349, -0.15314405,
2405                  -0.13573331, -0.15314679, -0.13560104, -0.15158523,
2406                  -0.13310701, -0.21208605, -0.19283913, -0.20955631,
2407                  -0.19134189, -0.20942821, -0.19134598, -0.20942799,
2408                  -0.1913458, -0.20943005, -0.19120952, -0.20781177,
2409                  -0.18869401, -0.25384082, -0.2463294, -0.25047649,
2410                  -0.24464654, -0.25031159, -0.24464253, -0.25031112,
2411                  -0.24464253, -0.25031463, -0.24454764, -0.24885323,
2412                  -0.24286438])
2413
2414
2415
2416    def test_bedslope_problem_second_order_more_steps(self):
2417
2418        from mesh_factory import rectangular
2419        from shallow_water import Domain, Reflective_boundary, Constant_height
2420        from Numeric import array
2421
2422        #Create basic mesh
2423        points, vertices, boundary = rectangular(6, 6)
2424
2425        #Create shallow water domain
2426        domain = Domain(points, vertices, boundary)
2427        domain.smooth = False
2428        domain.default_order=2
2429        domain.beta_h = 0.0 #Use first order in h-limiter       
2430
2431        #Bed-slope and friction at vertices (and interpolated elsewhere)
2432        def x_slope(x, y):
2433            return -x/3
2434
2435        domain.set_quantity('elevation', x_slope)
2436
2437        # Boundary conditions
2438        Br = Reflective_boundary(domain)
2439        domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
2440
2441        #Initial condition
2442        domain.set_quantity('stage', Constant_height(x_slope, 0.05))
2443        domain.check_integrity()
2444
2445        assert allclose(domain.quantities['stage'].centroid_values,
2446                        [0.01296296, 0.03148148, 0.01296296,
2447                        0.03148148, 0.01296296, 0.03148148,
2448                        0.01296296, 0.03148148, 0.01296296,
2449                        0.03148148, 0.01296296, 0.03148148,
2450                        -0.04259259, -0.02407407, -0.04259259,
2451                        -0.02407407, -0.04259259, -0.02407407,
2452                        -0.04259259, -0.02407407, -0.04259259,
2453                        -0.02407407, -0.04259259, -0.02407407,
2454                        -0.09814815, -0.07962963, -0.09814815,
2455                        -0.07962963, -0.09814815, -0.07962963,
2456                        -0.09814815, -0.07962963, -0.09814815,
2457                        -0.07962963, -0.09814815, -0.07962963,
2458                        -0.1537037 , -0.13518519, -0.1537037,
2459                        -0.13518519, -0.1537037, -0.13518519,
2460                        -0.1537037 , -0.13518519, -0.1537037,
2461                        -0.13518519, -0.1537037, -0.13518519,
2462                        -0.20925926, -0.19074074, -0.20925926,
2463                        -0.19074074, -0.20925926, -0.19074074,
2464                        -0.20925926, -0.19074074, -0.20925926,
2465                        -0.19074074, -0.20925926, -0.19074074,
2466                        -0.26481481, -0.2462963, -0.26481481,
2467                        -0.2462963, -0.26481481, -0.2462963,
2468                        -0.26481481, -0.2462963, -0.26481481,
2469                        -0.2462963, -0.26481481, -0.2462963])
2470
2471
2472        #print domain.quantities['stage'].extrapolate_second_order()
2473        #domain.distribute_to_vertices_and_edges()
2474        #print domain.quantities['stage'].vertex_values[:,0]       
2475       
2476        #Evolution
2477        for t in domain.evolve(yieldstep = 0.05, finaltime = 0.5):
2478            pass
2479
2480       
2481        assert allclose(domain.quantities['stage'].centroid_values,
2482     [-0.02907028, -0.01475478, -0.02973417, -0.01447186, -0.02932665, -0.01428285,
2483      -0.02901975, -0.0141361,  -0.02898816, -0.01418135, -0.02961409, -0.01403487,
2484      -0.07597998, -0.06252591, -0.07664854, -0.06312532, -0.07638287, -0.06265139,
2485      -0.07571145, -0.06235231, -0.0756817,  -0.06245309, -0.07652292, -0.06289946,
2486      -0.12367464, -0.11088981, -0.12237277, -0.11115338, -0.1218934,  -0.1107174, 
2487      -0.12081485, -0.11000491, -0.12038451, -0.11010335, -0.12102113, -0.11012105,
2488      -0.16909116, -0.15831543, -0.16730214, -0.15786249, -0.1665493,  -0.15697919,
2489      -0.16496618, -0.15559852, -0.16338679, -0.15509088, -0.16364092, -0.15424423,
2490      -0.18771107, -0.19903904, -0.18903759, -0.19858437, -0.18701552, -0.19697797,
2491      -0.1833593,  -0.19505871, -0.1818806,  -0.19418042, -0.18586159, -0.19576946,
2492      -0.13986873, -0.14170053, -0.14132188, -0.14560674, -0.14095617, -0.14373292,
2493      -0.13785933, -0.14033364, -0.13592955, -0.13936356, -0.13596008, -0.14216296])
2494
2495        assert allclose(domain.quantities['xmomentum'].centroid_values,
2496     [ 0.00831121,  0.00317948,  0.00731797,  0.00334939,  0.00764717,  0.00348053,
2497       0.00788729,  0.00356522,  0.00780649,  0.00341919,  0.00693525,  0.00310375,
2498       0.02166196,  0.01421475,  0.02017737,  0.01316839,  0.02037015,  0.01368659,
2499       0.02106,     0.01399161,  0.02074514,  0.01354935,  0.01887407,  0.0123113, 
2500       0.03775083,  0.02855197,  0.03689337,  0.02759782,  0.03732848,  0.02812072,
2501       0.03872545,  0.02913348,  0.03880939,  0.02803804,  0.03546499,  0.0260039, 
2502       0.0632131,   0.04730634,  0.0576324,   0.04592336,  0.05790921,  0.04690514,
2503       0.05986467,  0.04871165,  0.06170068,  0.04811572,  0.05657041,  0.04416292,
2504       0.08489642,  0.07188097,  0.07835261,  0.06843406,  0.07986412,  0.0698247, 
2505       0.08201071,  0.07216756,  0.08378418,  0.07273624,  0.080399,    0.06645841,
2506       0.01631548,  0.04691608,  0.0206632,   0.044409,    0.02115518,  0.04560305,
2507       0.02160608,  0.04663725,  0.02174734,  0.04795559,  0.02281427,  0.05667111])
2508
2509       
2510        assert allclose(domain.quantities['ymomentum'].centroid_values, 
2511     [ 1.45876601e-004, -3.24627393e-004, -1.57572719e-004, -2.92790187e-004,
2512      -9.90988382e-005, -3.06677335e-004, -1.62493106e-004, -3.71310004e-004,
2513      -1.99445058e-004, -3.28493467e-004,  6.68217349e-005, -8.42042805e-006,
2514       5.05093371e-004, -1.42842214e-004, -6.81454718e-005, -5.02084057e-004,
2515      -8.50583861e-005, -4.65443981e-004, -1.96406564e-004, -5.88889562e-004,
2516      -2.70160173e-004, -5.35485454e-004,  2.60780997e-004,  3.12145471e-005,
2517       5.16189608e-004,  1.07069062e-004,  9.29989252e-005, -3.71211119e-004,
2518       1.16350246e-004, -3.82407830e-004, -1.62077969e-004, -6.30906636e-004,
2519      -4.74025708e-004, -6.94463009e-004,  6.15092843e-005,  2.22106820e-004,
2520      -6.29589294e-004,  2.43611937e-004, -5.88125094e-004, -6.94293192e-005,
2521      -4.17914641e-004,  6.64609019e-005, -7.68334577e-004, -3.40232101e-004,
2522      -1.67424308e-003, -7.39485066e-004, -1.59966988e-003,  5.68262838e-005,
2523      -1.48470633e-003, -1.84554882e-003, -2.27200099e-003, -1.67506848e-003,
2524      -1.95610258e-003, -1.47638801e-003, -1.73779477e-003, -1.85498791e-003,
2525      -2.01357843e-003, -2.17675471e-003, -1.65783870e-003, -1.15818681e-003,
2526      -1.18663036e-003, -2.94229849e-003, -3.59309018e-003, -5.13496584e-003,
2527      -6.17359400e-003, -5.98761937e-003, -6.00540116e-003, -5.01121966e-003,
2528      -4.50964850e-003, -3.06319963e-003,  6.08950810e-004, -4.79537921e-004])
2529       
2530
2531       
2532
2533    def test_temp_play(self):
2534
2535        from mesh_factory import rectangular
2536        from shallow_water import Domain, Reflective_boundary, Constant_height
2537        from Numeric import array
2538
2539        #Create basic mesh
2540        points, vertices, boundary = rectangular(5, 5)
2541
2542        #Create shallow water domain
2543        domain = Domain(points, vertices, boundary)
2544        domain.smooth = False
2545        domain.default_order=2
2546        domain.beta_h = 0.0 #Use first order in h-limiter       
2547
2548        #Bed-slope and friction at vertices (and interpolated elsewhere)
2549        def x_slope(x, y):
2550            return -x/3
2551
2552        domain.set_quantity('elevation', x_slope)
2553
2554        # Boundary conditions
2555        Br = Reflective_boundary(domain)
2556        domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
2557
2558        #Initial condition
2559        domain.set_quantity('stage', Constant_height(x_slope, 0.05))
2560        domain.check_integrity()
2561
2562        #Evolution
2563        for t in domain.evolve(yieldstep = 0.05, finaltime = 0.1):
2564            pass
2565
2566        assert allclose(domain.quantities['stage'].centroid_values[:4],
2567                        [0.00206836, 0.01296714, 0.00363415, 0.01438924])
2568        assert allclose(domain.quantities['xmomentum'].centroid_values[:4],                     
2569                        [0.01360154, 0.00671133, 0.01264578, 0.00648503])
2570        assert allclose(domain.quantities['ymomentum'].centroid_values[:4],
2571                        [-1.19201077e-003, -7.23647546e-004, 
2572                        -6.39083123e-005, 6.29815168e-005])
2573       
2574       
2575    def test_complex_bed(self):
2576        #No friction is tested here
2577       
2578        from shallow_water import Domain, Reflective_boundary, Dirichlet_boundary,\
2579             Transmissive_boundary, Time_boundary,\
2580             Weir_simple as Weir, Constant_height
2581
2582        from mesh_factory import rectangular
2583        from Numeric import array
2584   
2585        N = 12
2586        points, vertices, boundary = rectangular(N, N/2, len1=1.2,len2=0.6,
2587                                                 origin=(-0.07, 0))
2588
2589
2590        domain = Domain(points, vertices, boundary)
2591        domain.smooth = False
2592        domain.visualise = False
2593        domain.default_order=2
2594
2595       
2596        inflow_stage = 0.1
2597        Z = Weir(inflow_stage)
2598        domain.set_quantity('elevation', Z)
2599
2600        Br = Reflective_boundary(domain)
2601        Bd = Dirichlet_boundary([inflow_stage, 0.0, 0.0])
2602        domain.set_boundary({'left': Bd, 'right': Br, 'bottom': Br, 'top': Br})
2603
2604        domain.set_quantity('stage', Constant_height(Z, 0.))
2605
2606        for t in domain.evolve(yieldstep = 0.02, finaltime = 0.2):
2607            pass
2608
2609
2610        #print domain.quantities['stage'].centroid_values
2611
2612        #FIXME: These numbers were from version before 25/10       
2613        #assert allclose(domain.quantities['stage'].centroid_values,
2614# [3.95822638e-002,  5.61022588e-002,  4.66437868e-002,  5.73081011e-002,
2615#  4.72394613e-002,  5.74684939e-002,  4.74309483e-002,  5.77458084e-002,
2616#  4.80628177e-002,  5.85656225e-002,  4.90498542e-002,  6.02609831e-002,
2617#  1.18470315e-002,  1.75136443e-002,  1.18035266e-002,  2.15565695e-002,
2618#  1.31620268e-002,  2.14351640e-002,  1.32351076e-002,  2.15450687e-002,
2619#  1.36414028e-002,  2.24274619e-002,  1.51689511e-002,  2.21789655e-002,
2620# -7.54337535e-003, -6.86362021e-004, -7.74146760e-003, -1.83756530e-003,
2621# -8.16773628e-003, -4.49916813e-004, -8.08202599e-003, -3.91118720e-004,
2622# -8.10292716e-003, -3.88584984e-004, -7.35226124e-003,  2.73985295e-004,
2623#  1.86166683e-001,  8.74070369e-002,  1.86166712e-001,  8.74035875e-002,
2624#  6.11666935e-002, -3.76173225e-002, -6.38333276e-002, -3.76147365e-002,
2625#  6.11666725e-002,  8.73846774e-002,  1.86166697e-001,  8.74171550e-002,
2626# -4.83333333e-002,  1.18333333e-001, -4.83333333e-002,  1.18333333e-001,
2627# -4.83333333e-002, -6.66666667e-003, -1.73333333e-001, -1.31666667e-001,
2628# -1.73333333e-001, -6.66666667e-003, -4.83333333e-002,  1.18333333e-001,
2629# -2.48333333e-001, -2.31666667e-001, -2.48333333e-001, -2.31666667e-001,
2630# -2.48333333e-001, -2.31666667e-001, -2.48333333e-001, -2.31666667e-001,
2631# -2.48333333e-001, -2.31666667e-001, -2.48333333e-001, -2.31666667e-001,
2632# -4.65000000e-001, -3.65000000e-001, -4.65000000e-001, -3.65000000e-001,
2633# -4.65000000e-001, -3.65000000e-001, -4.65000000e-001, -3.65000000e-001,
2634# -4.65000000e-001, -3.65000000e-001, -4.65000000e-001, -3.65000000e-001,
2635# -5.98333333e-001, -5.81666667e-001, -5.98333333e-001, -5.81666667e-001,
2636# -5.98333333e-001, -5.81666667e-001, -5.98333333e-001, -5.81666667e-001,
2637# -5.98333333e-001, -5.81666667e-001, -5.98333333e-001, -5.81666667e-001,
2638# -6.48333333e-001, -6.31666667e-001, -6.48333333e-001, -6.31666667e-001,
2639# -6.48333333e-001, -6.31666667e-001, -6.48333333e-001, -6.31666667e-001,
2640# -6.48333333e-001, -6.31666667e-001, -6.48333333e-001, -6.31666667e-001,
2641# -5.31666667e-001, -5.98333333e-001, -5.31666667e-001, -5.98333333e-001,
2642# -5.31666667e-001, -5.98333333e-001, -5.31666667e-001, -5.98333333e-001,
2643# -5.31666667e-001, -5.98333333e-001, -5.31666667e-001, -5.98333333e-001,
2644# -4.98333333e-001, -4.81666667e-001, -4.98333333e-001, -4.81666667e-001,
2645# -4.98333333e-001, -4.81666667e-001, -4.98333333e-001, -4.81666667e-001,
2646# -4.98333333e-001, -4.81666667e-001, -4.98333333e-001, -4.81666667e-001,
2647# -5.48333333e-001, -5.31666667e-001, -5.48333333e-001, -5.31666667e-001,
2648# -5.48333333e-001, -5.31666667e-001, -5.48333333e-001, -5.31666667e-001,
2649# -5.48333333e-001, -5.31666667e-001, -5.48333333e-001, -5.31666667e-001])
2650
2651
2652
2653    def test_spatio_temporal_boundary(self):
2654        """Test that boundary values can be read from file and interpolated
2655        in both time and space.
2656
2657        Verify that the same steady state solution is arrived at and that
2658        time interpolation works.
2659
2660        The full solution history is not exactly the same as
2661        file boundary must read and interpolate from *smoothed* version
2662        as stored in sww.
2663        """
2664        import time
2665       
2666        #Create sww file of simple propagation from left to right
2667        #through rectangular domain
2668       
2669        from mesh_factory import rectangular
2670
2671        #Create basic mesh
2672        points, vertices, boundary = rectangular(3, 3)
2673
2674        #Create shallow water domain
2675        domain1 = Domain(points, vertices, boundary)
2676       
2677        from util import mean
2678        domain1.reduction = mean
2679        domain1.smooth = False  #Exact result
2680
2681        domain1.default_order = 2
2682        domain1.store = True   
2683        domain1.set_datadir('.')
2684        domain1.set_name('spatio_temporal_boundary_source' + str(time.time()))
2685       
2686        #FIXME: This is extremely important!
2687        #How can we test if they weren't stored?
2688        domain1.quantities_to_be_stored = ['stage', 'xmomentum', 'ymomentum']
2689
2690               
2691        #Bed-slope and friction at vertices (and interpolated elsewhere)
2692        domain1.set_quantity('elevation', 0)
2693        domain1.set_quantity('friction', 0)     
2694
2695        # Boundary conditions
2696        Br = Reflective_boundary(domain1)
2697        Bd = Dirichlet_boundary([0.3,0,0])
2698        domain1.set_boundary({'left': Bd, 'top': Bd, 'right': Br, 'bottom': Br})
2699        #Initial condition
2700        domain1.set_quantity('stage', 0)
2701        domain1.check_integrity()
2702
2703        finaltime = 5
2704        #Evolution  (full domain - large steps)
2705        for t in domain1.evolve(yieldstep = 0.671, finaltime = finaltime):
2706            pass
2707            #domain1.write_time()
2708
2709        cv1 = domain1.quantities['stage'].centroid_values
2710       
2711       
2712        #Create an triangle shaped domain (reusing coordinates from domain 1),
2713        #formed from the lower and right hand  boundaries and
2714        #the sw-ne diagonal
2715        #from domain 1. Call it domain2
2716       
2717        points = [ [0,0], [1.0/3,0], [1.0/3,1.0/3],
2718                   [2.0/3,0], [2.0/3,1.0/3], [2.0/3,2.0/3],
2719                   [1,0], [1,1.0/3], [1,2.0/3], [1,1]]
2720       
2721        vertices = [ [1,2,0], 
2722                     [3,4,1], [2,1,4], [4,5,2], 
2723                     [6,7,3], [4,3,7], [7,8,4], [5,4,8], [8,9,5]]       
2724
2725        boundary = { (0,1):'bottom', (1,1):'bottom', (4,1): 'bottom',
2726                     (4,2):'right', (6,2):'right', (8,2):'right',
2727                     (0,0):'diagonal', (3,0):'diagonal', (8,0):'diagonal'}   
2728                     
2729        domain2 = Domain(points, vertices, boundary)
2730       
2731        domain2.reduction = domain1.reduction
2732        domain2.smooth = False
2733        domain2.default_order = 2
2734       
2735        #Bed-slope and friction at vertices (and interpolated elsewhere)
2736        domain2.set_quantity('elevation', 0)
2737        domain2.set_quantity('friction', 0)
2738        domain2.set_quantity('stage', 0)       
2739
2740        # Boundary conditions
2741        Br = Reflective_boundary(domain2)       
2742        Bf = Spatio_temporal_boundary(domain1.filename + '.' + domain1.format,
2743                                      domain2)
2744        domain2.set_boundary({'right':Br, 'bottom':Br, 'diagonal':Bf})
2745        domain2.check_integrity()       
2746       
2747
2748
2749        #Evolution (small steps)
2750        for t in domain2.evolve(yieldstep = 0.0711, finaltime = finaltime):
2751            pass
2752
2753       
2754        #Use output from domain1 as spatio-temporal boundary for domain2
2755        #and verify that results at right hand side are close. 
2756
2757        cv2 = domain2.quantities['stage'].centroid_values
2758
2759        #print take(cv1, (12,14,16))  #Right
2760        #print take(cv2, (4,6,8))
2761        #print take(cv1, (0,6,12))  #Bottom     
2762        #print take(cv2, (0,1,4))
2763        #print take(cv1, (0,8,16))  #Diag       
2764        #print take(cv2, (0,3,8))
2765
2766        assert allclose( take(cv1, (0,8,16)), take(cv2, (0,3,8))) #Diag
2767        assert allclose( take(cv1, (0,6,12)), take(cv2, (0,1,4))) #Bottom
2768        assert allclose( take(cv1, (12,14,16)), take(cv2, (4,6,8))) #RHS
2769
2770        #Cleanup
2771        os.remove(domain1.filename + '.' + domain1.format)       
2772       
2773       
2774
2775    def test_spatio_temporal_boundary_2(self):
2776        """Test that boundary values can be read from file and interpolated
2777        in both time and space.
2778        This is a more basic test, verifying that boundary object
2779        produces the expected results
2780
2781
2782        """
2783        import time
2784       
2785        #Create sww file of simple propagation from left to right
2786        #through rectangular domain
2787       
2788        from mesh_factory import rectangular
2789
2790        #Create basic mesh
2791        points, vertices, boundary = rectangular(3, 3)
2792
2793        #Create shallow water domain
2794        domain1 = Domain(points, vertices, boundary)
2795       
2796        from util import mean
2797        domain1.reduction = mean
2798        domain1.smooth = True #To mimic MOST output
2799
2800        domain1.default_order = 2
2801        domain1.store = True   
2802        domain1.set_datadir('.')
2803        domain1.set_name('spatio_temporal_boundary_source' + str(time.time()))
2804       
2805        #FIXME: This is extremely important!
2806        #How can we test if they weren't stored?
2807        domain1.quantities_to_be_stored = ['stage', 'xmomentum', 'ymomentum']
2808
2809               
2810        #Bed-slope and friction at vertices (and interpolated elsewhere)
2811        domain1.set_quantity('elevation', 0)
2812        domain1.set_quantity('friction', 0)     
2813
2814        # Boundary conditions
2815        Br = Reflective_boundary(domain1)
2816        Bd = Dirichlet_boundary([0.3,0,0])     
2817        domain1.set_boundary({'left': Bd, 'top': Bd, 'right': Br, 'bottom': Br})
2818        #Initial condition
2819        domain1.set_quantity('stage', 0)
2820        domain1.check_integrity()
2821
2822        finaltime = 5
2823        #Evolution  (full domain - large steps)
2824        for t in domain1.evolve(yieldstep = 1, finaltime = finaltime):
2825            pass
2826            #domain1.write_time()
2827
2828       
2829        #Create an triangle shaped domain (coinciding with some
2830        #coordinates from domain 1),
2831        #formed from the lower and right hand  boundaries and
2832        #the sw-ne diagonal
2833        #from domain 1. Call it domain2
2834       
2835        points = [ [0,0], [1.0/3,0], [1.0/3,1.0/3],
2836                   [2.0/3,0], [2.0/3,1.0/3], [2.0/3,2.0/3],
2837                   [1,0], [1,1.0/3], [1,2.0/3], [1,1]]
2838       
2839        vertices = [ [1,2,0], 
2840                     [3,4,1], [2,1,4], [4,5,2], 
2841                     [6,7,3], [4,3,7], [7,8,4], [5,4,8], [8,9,5]]       
2842
2843        boundary = { (0,1):'bottom', (1,1):'bottom', (4,1): 'bottom',
2844                     (4,2):'right', (6,2):'right', (8,2):'right',
2845                     (0,0):'diagonal', (3,0):'diagonal', (8,0):'diagonal'}   
2846                     
2847        domain2 = Domain(points, vertices, boundary)
2848       
2849        domain2.reduction = domain1.reduction
2850        domain2.smooth = False
2851        domain2.default_order = 2
2852       
2853        #Bed-slope and friction at vertices (and interpolated elsewhere)
2854        domain2.set_quantity('elevation', 0)
2855        domain2.set_quantity('friction', 0)
2856        domain2.set_quantity('stage', 0)       
2857
2858       
2859        #Read results for specific timesteps t=1 and t=2
2860        from Scientific.IO.NetCDF import NetCDFFile     
2861        fid = NetCDFFile(domain1.filename + '.' + domain1.format)
2862       
2863        x = fid.variables['x'][:]
2864        y = fid.variables['y'][:]       
2865        s1 = fid.variables['stage'][1,:]
2866        s2 = fid.variables['stage'][2,:]       
2867        fid.close()
2868
2869        from Numeric import take, reshape, concatenate 
2870        shp = (len(x), 1)
2871        points = concatenate( (reshape(x, shp), reshape(y, shp)), axis=1)       
2872        #The diagonal points of domain 1 are 0, 5, 10, 15
2873       
2874        #print points[0], points[5], points[10], points[15]
2875        assert allclose( take(points, [0,5,10,15]), 
2876                         [[0,0], [1.0/3, 1.0/3], [2.0/3, 2.0/3], [1,1]]) 
2877
2878       
2879        # Boundary conditions
2880        Br = Reflective_boundary(domain2)       
2881        Bf = Spatio_temporal_boundary(domain1.filename + '.' + domain1.format,
2882                                      domain2)
2883        domain2.set_boundary({'right':Br, 'bottom':Br, 'diagonal':Bf})
2884        domain2.check_integrity()       
2885
2886        #Test that interpolation points are the mid points of the all boundary
2887        #segments       
2888
2889        boundary_midpoints = [[1.0/6, 0], [1.0/2, 0], [5.0/6,0],
2890                              [1.0, 1.0/6], [1.0, 1.0/2], [1.0, 5.0/6],
2891                              [1.0/6, 1.0/6], [0.5, 0.5], [5.0/6, 5.0/6]]
2892                             
2893        boundary_midpoints.sort()
2894        R = Bf.F.interpolation_points.tolist() 
2895        R.sort()
2896        assert allclose(boundary_midpoints, R) 
2897       
2898        #Check spatially interpolated output at time == 1
2899        domain2.time = 1
2900       
2901        #First diagonal midpoint
2902        R0 = Bf.evaluate(0,0)
2903        assert allclose(R0[0], (s1[0] + s1[5])/2)
2904       
2905        #Second diagonal midpoint
2906        R0 = Bf.evaluate(3,0)
2907        assert allclose(R0[0], (s1[5] + s1[10])/2)
2908       
2909        #First diagonal midpoint
2910        R0 = Bf.evaluate(8,0)
2911        assert allclose(R0[0], (s1[10] + s1[15])/2)             
2912               
2913        #Check spatially interpolated output at time == 2
2914        domain2.time = 2
2915       
2916        #First diagonal midpoint
2917        R0 = Bf.evaluate(0,0)
2918        assert allclose(R0[0], (s2[0] + s2[5])/2)
2919       
2920        #Second diagonal midpoint
2921        R0 = Bf.evaluate(3,0)
2922        assert allclose(R0[0], (s2[5] + s2[10])/2)
2923       
2924        #First diagonal midpoint
2925        R0 = Bf.evaluate(8,0)
2926        assert allclose(R0[0], (s2[10] + s2[15])/2)             
2927       
2928       
2929        #Now check temporal interpolation
2930
2931        domain2.time = 1 + 2.0/3
2932       
2933        #First diagonal midpoint
2934        R0 = Bf.evaluate(0,0)
2935        assert allclose(R0[0], ((s1[0] + s1[5])/2 + 2.0*(s2[0] + s2[5])/2)/3)
2936       
2937        #Second diagonal midpoint
2938        R0 = Bf.evaluate(3,0)
2939        assert allclose(R0[0], ((s1[5] + s1[10])/2 + 2.0*(s2[5] + s2[10])/2)/3)
2940       
2941        #First diagonal midpoint
2942        R0 = Bf.evaluate(8,0)
2943        assert allclose(R0[0], ((s1[10] + s1[15])/2 + 2.0*(s2[10] + s2[15])/2)/3)               
2944               
2945       
2946               
2947        #Cleanup
2948        os.remove(domain1.filename + '.' + domain1.format)                 
2949#-------------------------------------------------------------
2950if __name__ == "__main__":
2951    suite = unittest.makeSuite(TestCase,'test')
2952    #suite = unittest.makeSuite(TestCase,'test_balance_deep_and_shallow')
2953    runner = unittest.TextTestRunner()
2954    runner.run(suite)
2955
Note: See TracBrowser for help on using the repository browser.