source: anuga_core/source/anuga/shallow_water/test_shallow_water_domain.py @ 3563

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

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