source: inundation/pyvolution/test_shallow_water.py @ 2709

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

Added algorithm for boundary_polygon in the presence of
multiple vertex coordinates. Also tested that new fit_interpolate
can use this new version.

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