source: anuga_core/source/anuga/pyvolution/test_shallow_water.py @ 3514

Last change on this file since 3514 was 3514, checked in by duncan, 18 years ago

Hi all,
I'm doing a change in the anuga structure, moving the code to

\anuga_core\source\anuga

After you have done an svn update, the PYTHONPATH has to be changed to;
PYTHONPATH = anuga_core/source/

This is part of changes required to make installation of anuga quicker and reducing the size of our sandpits.

If any imports are broken, try fixing them. With adding anuga. to them for example. If this seems to have really broken things, email/phone me.

Cheers
Duncan

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