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

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