# source:anuga_core/source/anuga/shallow_water/test_shallow_water_domain.py@3632

Last change on this file since 3632 was 3632, checked in by ole, 17 years ago

Changed gravity back to 9.8 as unit tests in shallow_water depend on it.
Change this later.

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