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

Last change on this file since 1669 was 1636, checked in by ole, 19 years ago

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