source: inundation/pyvolution/test_shallow_water.py @ 2544

Last change on this file since 2544 was 2526, checked in by ole, 19 years ago

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