source: inundation/pyvolution/test_shallow_water.py @ 2565

Last change on this file since 2565 was 2565, checked in by ole, 18 years ago

Fixed unit tests

File size: 111.0 KB
Line 
1#!/usr/bin/env python
2
3import unittest, os
4from math import sqrt, pi
5
6from config import g, epsilon
7from Numeric import allclose, array, zeros, ones, Float, take
8from 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*#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*#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
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        #FIXME: These numbers were from version before 21/3/6 -
2158        #could be recreated by setting maximum_allowed_speed to 0 maybe
2159        #assert allclose(domain.quantities['xmomentum'].vertex_values[:4,0],
2160        #                [ 0.00064835, 0.03685719, 0.00085073, 0.03687313])
2161       
2162        assert allclose(domain.quantities['xmomentum'].vertex_values[:4,0],
2163                        [ 0.00090581, 0.03685719, 0.00088303, 0.03687313])
2164                       
2165                       
2166
2167        #assert allclose(domain.quantities['xmomentum'].vertex_values[:4,0],
2168        #                [0.00090581,0.03685719,0.00088303,0.03687313])
2169
2170        assert allclose(domain.quantities['ymomentum'].vertex_values[:4,0],
2171                        [-0.00139463,0.0006156,-0.00060364,0.00061827])
2172
2173
2174        os.remove(domain.filename + '.sww')
2175
2176
2177    def test_flatbed_second_order_distribute(self):
2178        #Use real data from pyvolution 2
2179        #painfully setup and extracted.
2180        from mesh_factory import rectangular
2181        from shallow_water import Domain,\
2182             Reflective_boundary, Dirichlet_boundary
2183
2184        from Numeric import array
2185
2186        #Create basic mesh
2187        N = 8
2188        points, vertices, boundary = rectangular(N, N)
2189
2190        #Create shallow water domain
2191        domain = Domain(points, vertices, boundary)
2192        domain.smooth = False
2193        domain.visualise = False
2194        domain.default_order=domain.order=2
2195
2196        # Boundary conditions
2197        Br = Reflective_boundary(domain)
2198        Bd = Dirichlet_boundary([0.2,0.,0.])
2199
2200        domain.set_boundary({'left': Bd, 'right': Br, 'top': Br, 'bottom': Br})
2201        domain.check_integrity()
2202
2203
2204
2205        for V in [False, True]:
2206            if V:
2207                #Set centroids as if system had been evolved
2208                L = zeros(2*N*N, Float)
2209                L[:32] = [7.21205592e-003, 5.35214298e-002, 1.00910824e-002,
2210                          5.35439433e-002, 1.00910824e-002, 5.35439433e-002,
2211                          1.00910824e-002, 5.35439433e-002, 1.00910824e-002,
2212                          5.35439433e-002, 1.00910824e-002, 5.35439433e-002,
2213                          1.00910824e-002, 5.35393928e-002, 1.02344264e-002,
2214                          5.59605058e-002, 0.00000000e+000, 3.31027800e-004,
2215                          0.00000000e+000, 4.37962142e-005, 0.00000000e+000,
2216                          4.37962142e-005, 0.00000000e+000, 4.37962142e-005,
2217                          0.00000000e+000, 4.37962142e-005, 0.00000000e+000,
2218                          4.37962142e-005, 0.00000000e+000, 4.37962142e-005,
2219                          0.00000000e+000, 5.57305948e-005]
2220
2221                X = zeros(2*N*N, Float)
2222                X[:32] = [6.48351607e-003, 3.68571894e-002, 8.50733285e-003,
2223                          3.68731327e-002, 8.50733285e-003, 3.68731327e-002,
2224                          8.50733285e-003, 3.68731327e-002, 8.50733285e-003,
2225                          3.68731327e-002, 8.50733285e-003, 3.68731327e-002,
2226                          8.50733285e-003, 3.68693861e-002, 8.65220973e-003,
2227                          3.85055387e-002, 0.00000000e+000, 2.86060840e-004,
2228                          0.00000000e+000, 3.58905503e-005, 0.00000000e+000,
2229                          3.58905503e-005, 0.00000000e+000, 3.58905503e-005,
2230                          0.00000000e+000, 3.58905503e-005, 0.00000000e+000,
2231                          3.58905503e-005, 0.00000000e+000, 3.58905503e-005,
2232                          0.00000000e+000, 4.57662812e-005]
2233
2234                Y = zeros(2*N*N, Float)
2235                Y[:32]=[-1.39463104e-003, 6.15600298e-004, -6.03637382e-004,
2236                        6.18272251e-004, -6.03637382e-004, 6.18272251e-004,
2237                        -6.03637382e-004, 6.18272251e-004, -6.03637382e-004,
2238                        6.18272251e-004, -6.03637382e-004, 6.18272251e-004,
2239                        -6.03637382e-004, 6.18599320e-004, -6.74622797e-004,
2240                        -1.48934756e-004, 0.00000000e+000, -5.35079969e-005,
2241                        0.00000000e+000, -2.57264987e-005, 0.00000000e+000,
2242                        -2.57264987e-005, 0.00000000e+000, -2.57264987e-005,
2243                        0.00000000e+000, -2.57264987e-005, 0.00000000e+000,
2244                        -2.57264987e-005, 0.00000000e+000, -2.57264987e-005,
2245                        0.00000000e+000, -2.57635178e-005]
2246
2247
2248                domain.set_quantity('stage', L, location='centroids')
2249                domain.set_quantity('xmomentum', X, location='centroids')
2250                domain.set_quantity('ymomentum', Y, location='centroids')
2251
2252                domain.check_integrity()
2253            else:
2254                #Evolution
2255                for t in domain.evolve(yieldstep = 0.01, finaltime = 0.03):
2256                    pass
2257                assert allclose(domain.min_timestep, 0.0210448446782)
2258                assert allclose(domain.max_timestep, 0.0210448446782)
2259
2260
2261            #Centroids were correct but not vertices.
2262            #Hence the check of distribute below.
2263            assert allclose(domain.quantities['stage'].centroid_values[:4],
2264                            [0.00721206,0.05352143,0.01009108,0.05354394])
2265
2266            assert allclose(domain.quantities['xmomentum'].centroid_values[:4],
2267                            [0.00648352,0.03685719,0.00850733,0.03687313])
2268
2269            assert allclose(domain.quantities['ymomentum'].centroid_values[:4],
2270                            [-0.00139463,0.0006156,-0.00060364,0.00061827])
2271
2272            #print 'C17=', domain.quantities['xmomentum'].centroid_values[17]
2273            #print 'C19=', domain.quantities['xmomentum'].centroid_values[19]
2274
2275            #assert allclose(domain.quantities['xmomentum'].centroid_values[17],0.00028606084)
2276            ##print domain.quantities['xmomentum'].centroid_values[17], V
2277            ##print
2278            if not V:
2279                #FIXME: These numbers were from version before 21/3/6 -
2280                #could be recreated by setting maximum_allowed_speed to 0 maybe           
2281               
2282                #assert allclose(domain.quantities['xmomentum'].centroid_values[17], 0.0)               
2283                assert allclose(domain.quantities['xmomentum'].centroid_values[17], 0.000286060839592)                         
2284
2285            else:
2286                assert allclose(domain.quantities['xmomentum'].centroid_values[17], 0.00028606084)
2287
2288            import copy
2289            XX = copy.copy(domain.quantities['xmomentum'].centroid_values)
2290            assert allclose(domain.quantities['xmomentum'].centroid_values, XX)
2291
2292            domain.distribute_to_vertices_and_edges()
2293
2294            #assert allclose(domain.quantities['xmomentum'].centroid_values, XX)
2295
2296            #assert allclose(domain.quantities['xmomentum'].centroid_values[17],
2297            #                0.0)
2298            assert allclose(domain.quantities['xmomentum'].centroid_values[17], 0.000286060839592)                                                 
2299
2300
2301            #FIXME: These numbers were from version before 25/10
2302            #assert allclose(domain.quantities['stage'].vertex_values[:4,0],
2303            #                [0.00101913,0.05352143,0.00104852,0.05354394])
2304
2305            assert allclose(domain.quantities['ymomentum'].vertex_values[:4,0],
2306                            [-0.00139463,0.0006156,-0.00060364,0.00061827])
2307
2308
2309            assert allclose(domain.quantities['xmomentum'].vertex_values[:4,0],
2310                            [0.00090581,0.03685719,0.00088303,0.03687313])
2311
2312
2313            #NB NO longer relvant:
2314
2315            #This was the culprit. First triangles vertex 0 had an
2316            #x-momentum of 0.0064835 instead of 0.00090581 and
2317            #third triangle had 0.00850733 instead of 0.00088303
2318            #print domain.quantities['xmomentum'].vertex_values[:4,0]
2319
2320            #print domain.quantities['xmomentum'].vertex_values[:4,0]
2321            #assert allclose(domain.quantities['xmomentum'].vertex_values[:4,0],
2322            #                [0.00090581,0.03685719,0.00088303,0.03687313])
2323
2324        os.remove(domain.filename + '.sww')
2325
2326
2327
2328    def test_bedslope_problem_first_order(self):
2329
2330        from mesh_factory import rectangular
2331        from shallow_water import Domain, Reflective_boundary, Constant_height
2332        from Numeric import array
2333
2334        #Create basic mesh
2335        points, vertices, boundary = rectangular(6, 6)
2336
2337        #Create shallow water domain
2338        domain = Domain(points, vertices, boundary)
2339        domain.smooth = False
2340        domain.default_order=1
2341
2342        #Bed-slope and friction
2343        def x_slope(x, y):
2344            return -x/3
2345
2346        domain.set_quantity('elevation', x_slope)
2347
2348        # Boundary conditions
2349        Br = Reflective_boundary(domain)
2350        domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
2351
2352        #Initial condition
2353        domain.set_quantity('stage', Constant_height(x_slope, 0.05))
2354        domain.check_integrity()
2355
2356        #Evolution
2357        for t in domain.evolve(yieldstep = 0.05, finaltime = 0.05):
2358            pass# domain.write_time()
2359
2360        assert allclose(domain.min_timestep, 0.050010003001)
2361        assert allclose(domain.max_timestep, 0.050010003001)
2362
2363
2364        os.remove(domain.filename + '.sww')
2365       
2366    def test_bedslope_problem_first_order_moresteps(self):
2367
2368        from mesh_factory import rectangular
2369        from shallow_water import Domain, Reflective_boundary, Constant_height
2370        from Numeric import array
2371
2372        #Create basic mesh
2373        points, vertices, boundary = rectangular(6, 6)
2374
2375        #Create shallow water domain
2376        domain = Domain(points, vertices, boundary)
2377        domain.smooth = False
2378        domain.default_order=1
2379        domain.beta_h = 0.0 #Use first order in h-limiter
2380
2381        #Bed-slope and friction
2382        def x_slope(x, y):
2383            return -x/3
2384
2385        domain.set_quantity('elevation', x_slope)
2386
2387        # Boundary conditions
2388        Br = Reflective_boundary(domain)
2389        domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
2390
2391        #Initial condition
2392        domain.set_quantity('stage', Constant_height(x_slope, 0.05))
2393        domain.check_integrity()
2394
2395        #Evolution
2396        for t in domain.evolve(yieldstep = 0.05, finaltime = 0.5):
2397            pass# domain.write_time()
2398
2399        #Data from earlier version of pyvolution
2400        #print domain.quantities['stage'].centroid_values
2401
2402        assert allclose(domain.quantities['stage'].centroid_values,
2403                        [-0.02998628, -0.01520652, -0.03043492,
2404                        -0.0149132, -0.03004706, -0.01476251,
2405                        -0.0298215, -0.01467976, -0.02988158,
2406                        -0.01474662, -0.03036161, -0.01442995,
2407                        -0.07624583, -0.06297061, -0.07733792,
2408                        -0.06342237, -0.07695439, -0.06289595,
2409                        -0.07635559, -0.0626065, -0.07633628,
2410                        -0.06280072, -0.07739632, -0.06386738,
2411                        -0.12161738, -0.11028239, -0.1223796,
2412                        -0.11095953, -0.12189744, -0.11048616,
2413                        -0.12074535, -0.10987605, -0.12014311,
2414                        -0.10976691, -0.12096859, -0.11087692,
2415                        -0.16868259, -0.15868061, -0.16801135,
2416                        -0.1588003, -0.16674343, -0.15813323,
2417                        -0.16457595, -0.15693826, -0.16281096,
2418                        -0.15585154, -0.16283873, -0.15540068,
2419                        -0.17450362, -0.19919913, -0.18062882,
2420                        -0.19764131, -0.17783111, -0.19407213,
2421                        -0.1736915, -0.19053624, -0.17228678,
2422                        -0.19105634, -0.17920133, -0.1968828,
2423                        -0.14244395, -0.14604641, -0.14473537,
2424                        -0.1506107, -0.14510055, -0.14919522,
2425                        -0.14175896, -0.14560798, -0.13911658,
2426                        -0.14439383, -0.13924047, -0.14829043])
2427
2428        os.remove(domain.filename + '.sww')
2429       
2430    def test_bedslope_problem_second_order_one_step(self):
2431
2432        from mesh_factory import rectangular
2433        from shallow_water import Domain, Reflective_boundary, Constant_height
2434        from Numeric import array
2435
2436        #Create basic mesh
2437        points, vertices, boundary = rectangular(6, 6)
2438
2439        #Create shallow water domain
2440        domain = Domain(points, vertices, boundary)
2441        domain.smooth = False
2442        domain.default_order=2
2443
2444        #Bed-slope and friction at vertices (and interpolated elsewhere)
2445        def x_slope(x, y):
2446            return -x/3
2447
2448        domain.set_quantity('elevation', x_slope)
2449
2450        # Boundary conditions
2451        Br = Reflective_boundary(domain)
2452        domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
2453
2454        #Initial condition
2455        domain.set_quantity('stage', Constant_height(x_slope, 0.05))
2456        domain.check_integrity()
2457
2458        assert allclose(domain.quantities['stage'].centroid_values,
2459                        [0.01296296, 0.03148148, 0.01296296,
2460                        0.03148148, 0.01296296, 0.03148148,
2461                        0.01296296, 0.03148148, 0.01296296,
2462                        0.03148148, 0.01296296, 0.03148148,
2463                        -0.04259259, -0.02407407, -0.04259259,
2464                        -0.02407407, -0.04259259, -0.02407407,
2465                        -0.04259259, -0.02407407, -0.04259259,
2466                        -0.02407407, -0.04259259, -0.02407407,
2467                        -0.09814815, -0.07962963, -0.09814815,
2468                        -0.07962963, -0.09814815, -0.07962963,
2469                        -0.09814815, -0.07962963, -0.09814815,
2470                        -0.07962963, -0.09814815, -0.07962963,
2471                        -0.1537037 , -0.13518519, -0.1537037,
2472                        -0.13518519, -0.1537037, -0.13518519,
2473                        -0.1537037 , -0.13518519, -0.1537037,
2474                        -0.13518519, -0.1537037, -0.13518519,
2475                        -0.20925926, -0.19074074, -0.20925926,
2476                        -0.19074074, -0.20925926, -0.19074074,
2477                        -0.20925926, -0.19074074, -0.20925926,
2478                        -0.19074074, -0.20925926, -0.19074074,
2479                        -0.26481481, -0.2462963, -0.26481481,
2480                        -0.2462963, -0.26481481, -0.2462963,
2481                        -0.26481481, -0.2462963, -0.26481481,
2482                        -0.2462963, -0.26481481, -0.2462963])
2483
2484
2485        #print domain.quantities['stage'].extrapolate_second_order()
2486        #domain.distribute_to_vertices_and_edges()
2487        #print domain.quantities['stage'].vertex_values[:,0]
2488
2489        #Evolution
2490        for t in domain.evolve(yieldstep = 0.05, finaltime = 0.05):
2491            #domain.write_time()
2492            pass
2493
2494
2495        #print domain.quantities['stage'].centroid_values
2496        assert allclose(domain.quantities['stage'].centroid_values,
2497                 [0.01290985, 0.02356019, 0.01619096, 0.02356019, 0.01619096,
2498                  0.02356019, 0.01619096, 0.02356019, 0.01619096, 0.02356019,
2499                  0.01619096, 0.0268413, -0.04411074, -0.0248011, -0.04186556,
2500                  -0.0248011, -0.04186556, -0.0248011, -0.04186556, -0.0248011,
2501                  -0.04186556, -0.0248011, -0.04186556, -0.02255593,
2502                  -0.09966629, -0.08035666, -0.09742112, -0.08035666,
2503                  -0.09742112, -0.08035666, -0.09742112, -0.08035666,
2504                  -0.09742112, -0.08035666, -0.09742112, -0.07811149,
2505                  -0.15522185, -0.13591222, -0.15297667, -0.13591222,
2506                  -0.15297667, -0.13591222, -0.15297667, -0.13591222,
2507                  -0.15297667, -0.13591222, -0.15297667, -0.13366704,
2508                  -0.2107774, -0.19146777, -0.20853223, -0.19146777,
2509                  -0.20853223, -0.19146777, -0.20853223, -0.19146777,
2510                  -0.20853223, -0.19146777, -0.20853223, -0.1892226,
2511                  -0.26120669, -0.24776246, -0.25865535, -0.24776246,
2512                  -0.25865535, -0.24776246, -0.25865535, -0.24776246,
2513                  -0.25865535, -0.24776246, -0.25865535, -0.24521113])
2514
2515        os.remove(domain.filename + '.sww')
2516
2517    def test_bedslope_problem_second_order_two_steps(self):
2518
2519        from mesh_factory import rectangular
2520        from shallow_water import Domain, Reflective_boundary, Constant_height
2521        from Numeric import array
2522
2523        #Create basic mesh
2524        points, vertices, boundary = rectangular(6, 6)
2525
2526        #Create shallow water domain
2527        domain = Domain(points, vertices, boundary)
2528        domain.smooth = False
2529        domain.default_order=2
2530        domain.beta_h = 0.0 #Use first order in h-limiter
2531
2532        #Bed-slope and friction at vertices (and interpolated elsewhere)
2533        def x_slope(x, y):
2534            return -x/3
2535
2536        domain.set_quantity('elevation', x_slope)
2537
2538        # Boundary conditions
2539        Br = Reflective_boundary(domain)
2540        domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
2541
2542        #Initial condition
2543        domain.set_quantity('stage', Constant_height(x_slope, 0.05))
2544        domain.check_integrity()
2545
2546        assert allclose(domain.quantities['stage'].centroid_values,
2547                        [0.01296296, 0.03148148, 0.01296296,
2548                        0.03148148, 0.01296296, 0.03148148,
2549                        0.01296296, 0.03148148, 0.01296296,
2550                        0.03148148, 0.01296296, 0.03148148,
2551                        -0.04259259, -0.02407407, -0.04259259,
2552                        -0.02407407, -0.04259259, -0.02407407,
2553                        -0.04259259, -0.02407407, -0.04259259,
2554                        -0.02407407, -0.04259259, -0.02407407,
2555                        -0.09814815, -0.07962963, -0.09814815,
2556                        -0.07962963, -0.09814815, -0.07962963,
2557                        -0.09814815, -0.07962963, -0.09814815,
2558                        -0.07962963, -0.09814815, -0.07962963,
2559                        -0.1537037 , -0.13518519, -0.1537037,
2560                        -0.13518519, -0.1537037, -0.13518519,
2561                        -0.1537037 , -0.13518519, -0.1537037,
2562                        -0.13518519, -0.1537037, -0.13518519,
2563                        -0.20925926, -0.19074074, -0.20925926,
2564                        -0.19074074, -0.20925926, -0.19074074,
2565                        -0.20925926, -0.19074074, -0.20925926,
2566                        -0.19074074, -0.20925926, -0.19074074,
2567                        -0.26481481, -0.2462963, -0.26481481,
2568                        -0.2462963, -0.26481481, -0.2462963,
2569                        -0.26481481, -0.2462963, -0.26481481,
2570                        -0.2462963, -0.26481481, -0.2462963])
2571
2572
2573        #print domain.quantities['stage'].extrapolate_second_order()
2574        #domain.distribute_to_vertices_and_edges()
2575        #print domain.quantities['stage'].vertex_values[:,0]
2576
2577        #Evolution
2578        for t in domain.evolve(yieldstep = 0.05, finaltime = 0.1):
2579            pass
2580
2581
2582        #Data from earlier version of pyvolution ft=0.1
2583        assert allclose(domain.min_timestep, 0.0376895634803)
2584        assert allclose(domain.max_timestep, 0.0415635655309)
2585
2586
2587        assert allclose(domain.quantities['stage'].centroid_values,
2588                        [0.00855788, 0.01575204, 0.00994606, 0.01717072,
2589                         0.01005985, 0.01716362, 0.01005985, 0.01716299,
2590                         0.01007098, 0.01736248, 0.01216452, 0.02026776,
2591                         -0.04462374, -0.02479045, -0.04199789, -0.0229465,
2592                         -0.04184033, -0.02295693, -0.04184013, -0.02295675,
2593                         -0.04184486, -0.0228168, -0.04028876, -0.02036486,
2594                         -0.10029444, -0.08170809, -0.09772846, -0.08021704,
2595                         -0.09760006, -0.08022143, -0.09759984, -0.08022124,
2596                         -0.09760261, -0.08008893, -0.09603914, -0.07758209,
2597                         -0.15584152, -0.13723138, -0.15327266, -0.13572906,
2598                         -0.15314427, -0.13573349, -0.15314405, -0.13573331,
2599                         -0.15314679, -0.13560104, -0.15158523, -0.13310701,
2600                         -0.21208605, -0.19283913, -0.20955631, -0.19134189,
2601                         -0.20942821, -0.19134598, -0.20942799, -0.1913458,
2602                         -0.20943005, -0.19120952, -0.20781177, -0.18869401,
2603                         -0.25384082, -0.2463294, -0.25047649, -0.24464654,
2604                         -0.25031159, -0.24464253, -0.25031112, -0.24464253,
2605                         -0.25031463, -0.24454764, -0.24885323, -0.24286438])
2606
2607
2608        os.remove(domain.filename + '.sww')
2609
2610    def test_bedslope_problem_second_order_two_yieldsteps(self):
2611
2612        from mesh_factory import rectangular
2613        from shallow_water import Domain, Reflective_boundary, Constant_height
2614        from Numeric import array
2615
2616        #Create basic mesh
2617        points, vertices, boundary = rectangular(6, 6)
2618
2619        #Create shallow water domain
2620        domain = Domain(points, vertices, boundary)
2621        domain.smooth = False
2622        domain.default_order=2
2623        domain.beta_h = 0.0 #Use first order in h-limiter
2624
2625        #Bed-slope and friction at vertices (and interpolated elsewhere)
2626        def x_slope(x, y):
2627            return -x/3
2628
2629        domain.set_quantity('elevation', x_slope)
2630
2631        # Boundary conditions
2632        Br = Reflective_boundary(domain)
2633        domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
2634
2635        #Initial condition
2636        domain.set_quantity('stage', Constant_height(x_slope, 0.05))
2637        domain.check_integrity()
2638
2639        assert allclose(domain.quantities['stage'].centroid_values,
2640                        [0.01296296, 0.03148148, 0.01296296,
2641                        0.03148148, 0.01296296, 0.03148148,
2642                        0.01296296, 0.03148148, 0.01296296,
2643                        0.03148148, 0.01296296, 0.03148148,
2644                        -0.04259259, -0.02407407, -0.04259259,
2645                        -0.02407407, -0.04259259, -0.02407407,
2646                        -0.04259259, -0.02407407, -0.04259259,
2647                        -0.02407407, -0.04259259, -0.02407407,
2648                        -0.09814815, -0.07962963, -0.09814815,
2649                        -0.07962963, -0.09814815, -0.07962963,
2650                        -0.09814815, -0.07962963, -0.09814815,
2651                        -0.07962963, -0.09814815, -0.07962963,
2652                        -0.1537037 , -0.13518519, -0.1537037,
2653                        -0.13518519, -0.1537037, -0.13518519,
2654                        -0.1537037 , -0.13518519, -0.1537037,
2655                        -0.13518519, -0.1537037, -0.13518519,
2656                        -0.20925926, -0.19074074, -0.20925926,
2657                        -0.19074074, -0.20925926, -0.19074074,
2658                        -0.20925926, -0.19074074, -0.20925926,
2659                        -0.19074074, -0.20925926, -0.19074074,
2660                        -0.26481481, -0.2462963, -0.26481481,
2661                        -0.2462963, -0.26481481, -0.2462963,
2662                        -0.26481481, -0.2462963, -0.26481481,
2663                        -0.2462963, -0.26481481, -0.2462963])
2664
2665
2666        #print domain.quantities['stage'].extrapolate_second_order()
2667        #domain.distribute_to_vertices_and_edges()
2668        #print domain.quantities['stage'].vertex_values[:,0]
2669
2670        #Evolution
2671        for t in domain.evolve(yieldstep = 0.05, finaltime = 0.1):   #0.05??
2672            #domain.write_time()
2673            pass
2674
2675
2676
2677        assert allclose(domain.quantities['stage'].centroid_values,
2678                 [0.00855788, 0.01575204, 0.00994606, 0.01717072, 0.01005985,
2679                  0.01716362, 0.01005985, 0.01716299, 0.01007098, 0.01736248,
2680                  0.01216452, 0.02026776, -0.04462374, -0.02479045, -0.04199789,
2681                  -0.0229465, -0.04184033, -0.02295693, -0.04184013,
2682                  -0.02295675, -0.04184486, -0.0228168, -0.04028876,
2683                  -0.02036486, -0.10029444, -0.08170809, -0.09772846,
2684                  -0.08021704, -0.09760006, -0.08022143, -0.09759984,
2685                  -0.08022124, -0.09760261, -0.08008893, -0.09603914,
2686                  -0.07758209, -0.15584152, -0.13723138, -0.15327266,
2687                  -0.13572906, -0.15314427, -0.13573349, -0.15314405,
2688                  -0.13573331, -0.15314679, -0.13560104, -0.15158523,
2689                  -0.13310701, -0.21208605, -0.19283913, -0.20955631,
2690                  -0.19134189, -0.20942821, -0.19134598, -0.20942799,
2691                  -0.1913458, -0.20943005, -0.19120952, -0.20781177,
2692                  -0.18869401, -0.25384082, -0.2463294, -0.25047649,
2693                  -0.24464654, -0.25031159, -0.24464253, -0.25031112,
2694                  -0.24464253, -0.25031463, -0.24454764, -0.24885323,
2695                  -0.24286438])
2696
2697        os.remove(domain.filename + '.sww')
2698
2699    def test_bedslope_problem_second_order_more_steps(self):
2700
2701        from mesh_factory import rectangular
2702        from shallow_water import Domain, Reflective_boundary, Constant_height
2703        from Numeric import array
2704
2705        #Create basic mesh
2706        points, vertices, boundary = rectangular(6, 6)
2707
2708        #Create shallow water domain
2709        domain = Domain(points, vertices, boundary)
2710        domain.smooth = False
2711        domain.default_order=2
2712        domain.beta_h = 0.0 #Use first order in h-limiter
2713
2714        #Bed-slope and friction at vertices (and interpolated elsewhere)
2715        def x_slope(x, y):
2716            return -x/3
2717
2718        domain.set_quantity('elevation', x_slope)
2719
2720        # Boundary conditions
2721        Br = Reflective_boundary(domain)
2722        domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
2723
2724        #Initial condition
2725        domain.set_quantity('stage', Constant_height(x_slope, 0.05))
2726        domain.check_integrity()
2727
2728        assert allclose(domain.quantities['stage'].centroid_values,
2729                        [0.01296296, 0.03148148, 0.01296296,
2730                        0.03148148, 0.01296296, 0.03148148,
2731                        0.01296296, 0.03148148, 0.01296296,
2732                        0.03148148, 0.01296296, 0.03148148,
2733                        -0.04259259, -0.02407407, -0.04259259,
2734                        -0.02407407, -0.04259259, -0.02407407,
2735                        -0.04259259, -0.02407407, -0.04259259,
2736                        -0.02407407, -0.04259259, -0.02407407,
2737                        -0.09814815, -0.07962963, -0.09814815,
2738                        -0.07962963, -0.09814815, -0.07962963,
2739                        -0.09814815, -0.07962963, -0.09814815,
2740                        -0.07962963, -0.09814815, -0.07962963,
2741                        -0.1537037 , -0.13518519, -0.1537037,
2742                        -0.13518519, -0.1537037, -0.13518519,
2743                        -0.1537037 , -0.13518519, -0.1537037,
2744                        -0.13518519, -0.1537037, -0.13518519,
2745                        -0.20925926, -0.19074074, -0.20925926,
2746                        -0.19074074, -0.20925926, -0.19074074,
2747                        -0.20925926, -0.19074074, -0.20925926,
2748                        -0.19074074, -0.20925926, -0.19074074,
2749                        -0.26481481, -0.2462963, -0.26481481,
2750                        -0.2462963, -0.26481481, -0.2462963,
2751                        -0.26481481, -0.2462963, -0.26481481,
2752                        -0.2462963, -0.26481481, -0.2462963])
2753
2754
2755        #print domain.quantities['stage'].extrapolate_second_order()
2756        #domain.distribute_to_vertices_and_edges()
2757        #print domain.quantities['stage'].vertex_values[:,0]
2758
2759        #Evolution
2760        for t in domain.evolve(yieldstep = 0.05, finaltime = 0.5):
2761            pass
2762
2763
2764        assert allclose(domain.quantities['stage'].centroid_values,
2765     [-0.02907028, -0.01475478, -0.02973417, -0.01447186, -0.02932665, -0.01428285,
2766      -0.02901975, -0.0141361,  -0.02898816, -0.01418135, -0.02961409, -0.01403487,
2767      -0.07597998, -0.06252591, -0.07664854, -0.06312532, -0.07638287, -0.06265139,
2768      -0.07571145, -0.06235231, -0.0756817,  -0.06245309, -0.07652292, -0.06289946,
2769      -0.12367464, -0.11088981, -0.12237277, -0.11115338, -0.1218934,  -0.1107174,
2770      -0.12081485, -0.11000491, -0.12038451, -0.11010335, -0.12102113, -0.11012105,
2771      -0.16909116, -0.15831543, -0.16730214, -0.15786249, -0.1665493,  -0.15697919,
2772      -0.16496618, -0.15559852, -0.16338679, -0.15509088, -0.16364092, -0.15424423,
2773      -0.18771107, -0.19903904, -0.18903759, -0.19858437, -0.18701552, -0.19697797,
2774      -0.1833593,  -0.19505871, -0.1818806,  -0.19418042, -0.18586159, -0.19576946,
2775      -0.13986873, -0.14170053, -0.14132188, -0.14560674, -0.14095617, -0.14373292,
2776      -0.13785933, -0.14033364, -0.13592955, -0.13936356, -0.13596008, -0.14216296])
2777
2778        assert allclose(domain.quantities['xmomentum'].centroid_values,
2779     [ 0.00831121,  0.00317948,  0.00731797,  0.00334939,  0.00764717,  0.00348053,
2780       0.00788729,  0.00356522,  0.00780649,  0.00341919,  0.00693525,  0.00310375,
2781       0.02166196,  0.01421475,  0.02017737,  0.01316839,  0.02037015,  0.01368659,
2782       0.02106,     0.01399161,  0.02074514,  0.01354935,  0.01887407,  0.0123113,
2783       0.03775083,  0.02855197,  0.03689337,  0.02759782,  0.03732848,  0.02812072,
2784       0.03872545,  0.02913348,  0.03880939,  0.02803804,  0.03546499,  0.0260039,
2785       0.0632131,   0.04730634,  0.0576324,   0.04592336,  0.05790921,  0.04690514,
2786       0.05986467,  0.04871165,  0.06170068,  0.04811572,  0.05657041,  0.04416292,
2787       0.08489642,  0.07188097,  0.07835261,  0.06843406,  0.07986412,  0.0698247,
2788       0.08201071,  0.07216756,  0.08378418,  0.07273624,  0.080399,    0.06645841,
2789       0.01631548,  0.04691608,  0.0206632,   0.044409,    0.02115518,  0.04560305,
2790       0.02160608,  0.04663725,  0.02174734,  0.04795559,  0.02281427,  0.05667111])
2791
2792
2793        assert allclose(domain.quantities['ymomentum'].centroid_values,
2794     [ 1.45876601e-004, -3.24627393e-004, -1.57572719e-004, -2.92790187e-004,
2795      -9.90988382e-005, -3.06677335e-004, -1.62493106e-004, -3.71310004e-004,
2796      -1.99445058e-004, -3.28493467e-004,  6.68217349e-005, -8.42042805e-006,
2797       5.05093371e-004, -1.42842214e-004, -6.81454718e-005, -5.02084057e-004,
2798      -8.50583861e-005, -4.65443981e-004, -1.96406564e-004, -5.88889562e-004,
2799      -2.70160173e-004, -5.35485454e-004,  2.60780997e-004,  3.12145471e-005,
2800       5.16189608e-004,  1.07069062e-004,  9.29989252e-005, -3.71211119e-004,
2801       1.16350246e-004, -3.82407830e-004, -1.62077969e-004, -6.30906636e-004,
2802      -4.74025708e-004, -6.94463009e-004,  6.15092843e-005,  2.22106820e-004,
2803      -6.29589294e-004,  2.43611937e-004, -5.88125094e-004, -6.94293192e-005,
2804      -4.17914641e-004,  6.64609019e-005, -7.68334577e-004, -3.40232101e-004,
2805      -1.67424308e-003, -7.39485066e-004, -1.59966988e-003,  5.68262838e-005,
2806      -1.48470633e-003, -1.84554882e-003, -2.27200099e-003, -1.67506848e-003,
2807      -1.95610258e-003, -1.47638801e-003, -1.73779477e-003, -1.85498791e-003,
2808      -2.01357843e-003, -2.17675471e-003, -1.65783870e-003, -1.15818681e-003,
2809      -1.18663036e-003, -2.94229849e-003, -3.59309018e-003, -5.13496584e-003,
2810      -6.17359400e-003, -5.98761937e-003, -6.00540116e-003, -5.01121966e-003,
2811      -4.50964850e-003, -3.06319963e-003,  6.08950810e-004, -4.79537921e-004])
2812
2813        os.remove(domain.filename + '.sww')
2814
2815
2816    def test_temp_play(self):
2817
2818        from mesh_factory import rectangular
2819        from shallow_water import Domain, Reflective_boundary, Constant_height
2820        from Numeric import array
2821
2822        #Create basic mesh
2823        points, vertices, boundary = rectangular(5, 5)
2824
2825        #Create shallow water domain
2826        domain = Domain(points, vertices, boundary)
2827        domain.smooth = False
2828        domain.default_order=2
2829        domain.beta_h = 0.0 #Use first order in h-limiter
2830
2831        #Bed-slope and friction at vertices (and interpolated elsewhere)
2832        def x_slope(x, y):
2833            return -x/3
2834
2835        domain.set_quantity('elevation', x_slope)
2836
2837        # Boundary conditions
2838        Br = Reflective_boundary(domain)
2839        domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
2840
2841        #Initial condition
2842        domain.set_quantity('stage', Constant_height(x_slope, 0.05))
2843        domain.check_integrity()
2844
2845        #Evolution
2846        for t in domain.evolve(yieldstep = 0.05, finaltime = 0.1):
2847            pass
2848
2849        assert allclose(domain.quantities['stage'].centroid_values[:4],
2850                        [0.00206836, 0.01296714, 0.00363415, 0.01438924])
2851        assert allclose(domain.quantities['xmomentum'].centroid_values[:4],
2852                        [0.01360154, 0.00671133, 0.01264578, 0.00648503])
2853        assert allclose(domain.quantities['ymomentum'].centroid_values[:4],
2854                        [-1.19201077e-003, -7.23647546e-004,
2855                        -6.39083123e-005, 6.29815168e-005])
2856
2857        os.remove(domain.filename + '.sww')
2858
2859    def test_complex_bed(self):
2860        #No friction is tested here
2861
2862        from shallow_water import Domain, Reflective_boundary, Dirichlet_boundary,\
2863             Transmissive_boundary, Time_boundary,\
2864             Weir_simple as Weir, Constant_height
2865
2866        from mesh_factory import rectangular
2867        from Numeric import array
2868
2869        N = 12
2870        points, vertices, boundary = rectangular(N, N/2, len1=1.2,len2=0.6,
2871                                                 origin=(-0.07, 0))
2872
2873
2874        domain = Domain(points, vertices, boundary)
2875        domain.smooth = False
2876        domain.visualise = False
2877        domain.default_order=2
2878
2879
2880        inflow_stage = 0.1
2881        Z = Weir(inflow_stage)
2882        domain.set_quantity('elevation', Z)
2883
2884        Br = Reflective_boundary(domain)
2885        Bd = Dirichlet_boundary([inflow_stage, 0.0, 0.0])
2886        domain.set_boundary({'left': Bd, 'right': Br, 'bottom': Br, 'top': Br})
2887
2888        domain.set_quantity('stage', Constant_height(Z, 0.))
2889
2890        for t in domain.evolve(yieldstep = 0.02, finaltime = 0.2):
2891            pass
2892
2893
2894        #print domain.quantities['stage'].centroid_values
2895
2896        #FIXME: These numbers were from version before 25/10
2897        #assert allclose(domain.quantities['stage'].centroid_values,
2898# [3.95822638e-002,  5.61022588e-002,  4.66437868e-002,  5.73081011e-002,
2899#  4.72394613e-002,  5.74684939e-002,  4.74309483e-002,  5.77458084e-002,
2900#  4.80628177e-002,  5.85656225e-002,  4.90498542e-002,  6.02609831e-002,
2901#  1.18470315e-002,  1.75136443e-002,  1.18035266e-002,  2.15565695e-002,
2902#  1.31620268e-002,  2.14351640e-002,  1.32351076e-002,  2.15450687e-002,
2903#  1.36414028e-002,  2.24274619e-002,  1.51689511e-002,  2.21789655e-002,
2904# -7.54337535e-003, -6.86362021e-004, -7.74146760e-003, -1.83756530e-003,
2905# -8.16773628e-003, -4.49916813e-004, -8.08202599e-003, -3.91118720e-004,
2906# -8.10292716e-003, -3.88584984e-004, -7.35226124e-003,  2.73985295e-004,
2907#  1.86166683e-001,  8.74070369e-002,  1.86166712e-001,  8.74035875e-002,
2908#  6.11666935e-002, -3.76173225e-002, -6.38333276e-002, -3.76147365e-002,
2909#  6.11666725e-002,  8.73846774e-002,  1.86166697e-001,  8.74171550e-002,
2910# -4.83333333e-002,  1.18333333e-001, -4.83333333e-002,  1.18333333e-001,
2911# -4.83333333e-002, -6.66666667e-003, -1.73333333e-001, -1.31666667e-001,
2912# -1.73333333e-001, -6.66666667e-003, -4.83333333e-002,  1.18333333e-001,
2913# -2.48333333e-001, -2.31666667e-001, -2.48333333e-001, -2.31666667e-001,
2914# -2.48333333e-001, -2.31666667e-001, -2.48333333e-001, -2.31666667e-001,
2915# -2.48333333e-001, -2.31666667e-001, -2.48333333e-001, -2.31666667e-001,
2916# -4.65000000e-001, -3.65000000e-001, -4.65000000e-001, -3.65000000e-001,
2917# -4.65000000e-001, -3.65000000e-001, -4.65000000e-001, -3.65000000e-001,
2918# -4.65000000e-001, -3.65000000e-001, -4.65000000e-001, -3.65000000e-001,
2919# -5.98333333e-001, -5.81666667e-001, -5.98333333e-001, -5.81666667e-001,
2920# -5.98333333e-001, -5.81666667e-001, -5.98333333e-001, -5.81666667e-001,
2921# -5.98333333e-001, -5.81666667e-001, -5.98333333e-001, -5.81666667e-001,
2922# -6.48333333e-001, -6.31666667e-001, -6.48333333e-001, -6.31666667e-001,
2923# -6.48333333e-001, -6.31666667e-001, -6.48333333e-001, -6.31666667e-001,
2924# -6.48333333e-001, -6.31666667e-001, -6.48333333e-001, -6.31666667e-001,
2925# -5.31666667e-001, -5.98333333e-001, -5.31666667e-001, -5.98333333e-001,
2926# -5.31666667e-001, -5.98333333e-001, -5.31666667e-001, -5.98333333e-001,
2927# -5.31666667e-001, -5.98333333e-001, -5.31666667e-001, -5.98333333e-001,
2928# -4.98333333e-001, -4.81666667e-001, -4.98333333e-001, -4.81666667e-001,
2929# -4.98333333e-001, -4.81666667e-001, -4.98333333e-001, -4.81666667e-001,
2930# -4.98333333e-001, -4.81666667e-001, -4.98333333e-001, -4.81666667e-001,
2931# -5.48333333e-001, -5.31666667e-001, -5.48333333e-001, -5.31666667e-001,
2932# -5.48333333e-001, -5.31666667e-001, -5.48333333e-001, -5.31666667e-001,
2933# -5.48333333e-001, -5.31666667e-001, -5.48333333e-001, -5.31666667e-001])
2934
2935        os.remove(domain.filename + '.sww')
2936
2937    def test_spatio_temporal_boundary_1(self):
2938        """Test that boundary values can be read from file and interpolated
2939        in both time and space.
2940
2941        Verify that the same steady state solution is arrived at and that
2942        time interpolation works.
2943
2944        The full solution history is not exactly the same as
2945        file boundary must read and interpolate from *smoothed* version
2946        as stored in sww.
2947        """
2948        import time
2949
2950        #Create sww file of simple propagation from left to right
2951        #through rectangular domain
2952
2953        from mesh_factory import rectangular
2954
2955        #Create basic mesh
2956        points, vertices, boundary = rectangular(3, 3)
2957
2958        #Create shallow water domain
2959        domain1 = Domain(points, vertices, boundary)
2960
2961        domain1.reduction = mean
2962        domain1.smooth = False  #Exact result
2963
2964        domain1.default_order = 2
2965        domain1.store = True
2966        domain1.set_datadir('.')
2967        domain1.set_name('spatio_temporal_boundary_source' + str(time.time()))
2968
2969        #FIXME: This is extremely important!
2970        #How can we test if they weren't stored?
2971        domain1.quantities_to_be_stored = ['stage', 'xmomentum', 'ymomentum']
2972
2973
2974        #Bed-slope and friction at vertices (and interpolated elsewhere)
2975        domain1.set_quantity('elevation', 0)
2976        domain1.set_quantity('friction', 0)
2977
2978        # Boundary conditions
2979        Br = Reflective_boundary(domain1)
2980        Bd = Dirichlet_boundary([0.3,0,0])
2981        domain1.set_boundary({'left': Bd, 'top': Bd, 'right': Br, 'bottom': Br})
2982        #Initial condition
2983        domain1.set_quantity('stage', 0)
2984        domain1.check_integrity()
2985
2986        finaltime = 5
2987        #Evolution  (full domain - large steps)
2988        for t in domain1.evolve(yieldstep = 0.671, finaltime = finaltime):
2989            pass
2990            #domain1.write_time()
2991
2992        cv1 = domain1.quantities['stage'].centroid_values
2993
2994
2995        #Create an triangle shaped domain (reusing coordinates from domain 1),
2996        #formed from the lower and right hand  boundaries and
2997        #the sw-ne diagonal
2998        #from domain 1. Call it domain2
2999
3000        points = [ [0,0], [1.0/3,0], [1.0/3,1.0/3],
3001                   [2.0/3,0], [2.0/3,1.0/3], [2.0/3,2.0/3],
3002                   [1,0], [1,1.0/3], [1,2.0/3], [1,1]]
3003
3004        vertices = [ [1,2,0], [3,4,1], [2,1,4], [4,5,2],
3005                     [6,7,3], [4,3,7], [7,8,4], [5,4,8], [8,9,5]]
3006
3007        boundary = { (0,1):'bottom', (1,1):'bottom', (4,1): 'bottom',
3008                     (4,2):'right', (6,2):'right', (8,2):'right',
3009                     (0,0):'diagonal', (3,0):'diagonal', (8,0):'diagonal'}
3010
3011        domain2 = Domain(points, vertices, boundary)
3012
3013        domain2.reduction = domain1.reduction
3014        domain2.smooth = False
3015        domain2.default_order = 2
3016
3017        #Bed-slope and friction at vertices (and interpolated elsewhere)
3018        domain2.set_quantity('elevation', 0)
3019        domain2.set_quantity('friction', 0)
3020        domain2.set_quantity('stage', 0)
3021
3022        # Boundary conditions
3023        Br = Reflective_boundary(domain2)
3024        Bf = Spatio_temporal_boundary(domain1.filename + '.' + domain1.format,
3025                                      domain2)
3026        domain2.set_boundary({'right':Br, 'bottom':Br, 'diagonal':Bf})
3027        domain2.check_integrity()
3028
3029
3030
3031        #Evolution (small steps)
3032        for t in domain2.evolve(yieldstep = 0.0711, finaltime = finaltime):
3033            pass
3034
3035
3036        #Use output from domain1 as spatio-temporal boundary for domain2
3037        #and verify that results at right hand side are close.
3038
3039        cv2 = domain2.quantities['stage'].centroid_values
3040
3041        #print take(cv1, (12,14,16))  #Right
3042        #print take(cv2, (4,6,8))
3043        #print take(cv1, (0,6,12))  #Bottom
3044        #print take(cv2, (0,1,4))
3045        #print take(cv1, (0,8,16))  #Diag
3046        #print take(cv2, (0,3,8))
3047
3048        assert allclose( take(cv1, (0,8,16)), take(cv2, (0,3,8))) #Diag
3049        assert allclose( take(cv1, (0,6,12)), take(cv2, (0,1,4))) #Bottom
3050        assert allclose( take(cv1, (12,14,16)), take(cv2, (4,6,8))) #RHS
3051
3052        #Cleanup
3053        os.remove(domain1.filename + '.' + domain1.format)
3054        os.remove(domain2.filename + '.' + domain2.format)       
3055
3056
3057
3058    def test_spatio_temporal_boundary_2(self):
3059        """Test that boundary values can be read from file and interpolated
3060        in both time and space.
3061        This is a more basic test, verifying that boundary object
3062        produces the expected results
3063
3064
3065        """
3066        import time
3067
3068        #Create sww file of simple propagation from left to right
3069        #through rectangular domain
3070
3071        from mesh_factory import rectangular
3072
3073        #Create basic mesh
3074        points, vertices, boundary = rectangular(3, 3)
3075
3076        #Create shallow water domain
3077        domain1 = Domain(points, vertices, boundary)
3078
3079        domain1.reduction = mean
3080        domain1.smooth = True #To mimic MOST output
3081
3082        domain1.default_order = 2
3083        domain1.store = True
3084        domain1.set_datadir('.')
3085        domain1.set_name('spatio_temporal_boundary_source' + str(time.time()))
3086
3087        #FIXME: This is extremely important!
3088        #How can we test if they weren't stored?
3089        domain1.quantities_to_be_stored = ['stage', 'xmomentum', 'ymomentum']
3090
3091
3092        #Bed-slope and friction at vertices (and interpolated elsewhere)
3093        domain1.set_quantity('elevation', 0)
3094        domain1.set_quantity('friction', 0)
3095
3096        # Boundary conditions
3097        Br = Reflective_boundary(domain1)
3098        Bd = Dirichlet_boundary([0.3,0,0])
3099        domain1.set_boundary({'left': Bd, 'top': Bd, 'right': Br, 'bottom': Br})
3100        #Initial condition
3101        domain1.set_quantity('stage', 0)
3102        domain1.check_integrity()
3103
3104        finaltime = 5
3105        #Evolution  (full domain - large steps)
3106        for t in domain1.evolve(yieldstep = 1, finaltime = finaltime):
3107            pass
3108            #domain1.write_time()
3109
3110
3111        #Create an triangle shaped domain (coinciding with some
3112        #coordinates from domain 1),
3113        #formed from the lower and right hand  boundaries and
3114        #the sw-ne diagonal
3115        #from domain 1. Call it domain2
3116
3117        points = [ [0,0], [1.0/3,0], [1.0/3,1.0/3],
3118                   [2.0/3,0], [2.0/3,1.0/3], [2.0/3,2.0/3],
3119                   [1,0], [1,1.0/3], [1,2.0/3], [1,1]]
3120
3121        vertices = [ [1,2,0],
3122                     [3,4,1], [2,1,4], [4,5,2],
3123                     [6,7,3], [4,3,7], [7,8,4], [5,4,8], [8,9,5]]
3124
3125        boundary = { (0,1):'bottom', (1,1):'bottom', (4,1): 'bottom',
3126                     (4,2):'right', (6,2):'right', (8,2):'right',
3127                     (0,0):'diagonal', (3,0):'diagonal', (8,0):'diagonal'}
3128
3129        domain2 = Domain(points, vertices, boundary)
3130
3131        domain2.reduction = domain1.reduction
3132        domain2.smooth = False
3133        domain2.default_order = 2
3134
3135        #Bed-slope and friction at vertices (and interpolated elsewhere)
3136        domain2.set_quantity('elevation', 0)
3137        domain2.set_quantity('friction', 0)
3138        domain2.set_quantity('stage', 0)
3139
3140
3141        #Read results for specific timesteps t=1 and t=2
3142        from Scientific.IO.NetCDF import NetCDFFile
3143        fid = NetCDFFile(domain1.filename + '.' + domain1.format)
3144
3145        x = fid.variables['x'][:]
3146        y = fid.variables['y'][:]
3147        s1 = fid.variables['stage'][1,:]
3148        s2 = fid.variables['stage'][2,:]
3149        fid.close()
3150
3151        from Numeric import take, reshape, concatenate
3152        shp = (len(x), 1)
3153        points = concatenate( (reshape(x, shp), reshape(y, shp)), axis=1)
3154        #The diagonal points of domain 1 are 0, 5, 10, 15
3155
3156        #print points[0], points[5], points[10], points[15]
3157        assert allclose( take(points, [0,5,10,15]),
3158                         [[0,0], [1.0/3, 1.0/3], [2.0/3, 2.0/3], [1,1]])
3159
3160
3161        # Boundary conditions
3162        Br = Reflective_boundary(domain2)
3163        Bf = Spatio_temporal_boundary(domain1.filename + '.' + domain1.format,
3164                                      domain2)
3165        domain2.set_boundary({'right':Br, 'bottom':Br, 'diagonal':Bf})
3166        domain2.check_integrity()
3167
3168        #Test that interpolation points are the mid points of the all boundary
3169        #segments
3170
3171        boundary_midpoints = [[1.0/6, 0], [1.0/2, 0], [5.0/6,0],
3172                              [1.0, 1.0/6], [1.0, 1.0/2], [1.0, 5.0/6],
3173                              [1.0/6, 1.0/6], [0.5, 0.5], [5.0/6, 5.0/6]]
3174
3175        boundary_midpoints.sort()
3176        R = Bf.F.interpolation_points.tolist()
3177        R.sort()
3178        assert allclose(boundary_midpoints, R)
3179
3180        #Check spatially interpolated output at time == 1
3181        domain2.time = 1
3182
3183        #First diagonal midpoint
3184        R0 = Bf.evaluate(0,0)
3185        assert allclose(R0[0], (s1[0] + s1[5])/2)
3186
3187        #Second diagonal midpoint
3188        R0 = Bf.evaluate(3,0)
3189        assert allclose(R0[0], (s1[5] + s1[10])/2)
3190
3191        #First diagonal midpoint
3192        R0 = Bf.evaluate(8,0)
3193        assert allclose(R0[0], (s1[10] + s1[15])/2)
3194
3195        #Check spatially interpolated output at time == 2
3196        domain2.time = 2
3197
3198        #First diagonal midpoint
3199        R0 = Bf.evaluate(0,0)
3200        assert allclose(R0[0], (s2[0] + s2[5])/2)
3201
3202        #Second diagonal midpoint
3203        R0 = Bf.evaluate(3,0)
3204        assert allclose(R0[0], (s2[5] + s2[10])/2)
3205
3206        #First diagonal midpoint
3207        R0 = Bf.evaluate(8,0)
3208        assert allclose(R0[0], (s2[10] + s2[15])/2)
3209
3210
3211        #Now check temporal interpolation
3212
3213        domain2.time = 1 + 2.0/3
3214
3215        #First diagonal midpoint
3216        R0 = Bf.evaluate(0,0)
3217        assert allclose(R0[0], ((s1[0] + s1[5])/2 + 2.0*(s2[0] + s2[5])/2)/3)
3218
3219        #Second diagonal midpoint
3220        R0 = Bf.evaluate(3,0)
3221        assert allclose(R0[0], ((s1[5] + s1[10])/2 + 2.0*(s2[5] + s2[10])/2)/3)
3222
3223        #First diagonal midpoint
3224        R0 = Bf.evaluate(8,0)
3225        assert allclose(R0[0], ((s1[10] + s1[15])/2 + 2.0*(s2[10] + s2[15])/2)/3)
3226
3227
3228
3229        #Cleanup
3230        os.remove(domain1.filename + '.' + domain1.format)
3231
3232
3233        #-------------------------------------------------------------
3234if __name__ == "__main__":
3235    suite = unittest.makeSuite(Test_Shallow_Water,'test')
3236    runner = unittest.TextTestRunner()
3237    runner.run(suite)
Note: See TracBrowser for help on using the repository browser.