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

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

Created new boundary class: Field_boundary designed to replace the generic File_boundary. Field_boundary is specific to the shallow water wave equation and allows mean_stage to be specified adjusting the stage derived from the sww. This should address ticket:132.
On test was added, two updated.

File size: 155.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*#New array
26
27    for k in range(N):
28
29        r = sqrt(x[k]**2 + y[k]**2)
30
31        factor = exp( -(r-0.15)**2 )
32
33        s[k] = 4000 * factor * (cos(t*2*pi/150) + 2)
34
35    return s
36
37
38def scalar_func(t,x,y):
39    """Function that returns a scalar.
40    Used to test error message when Numeric array is expected
41    """
42
43    return 17.7
44
45
46def angle(t,x,y):
47    """Rotating field
48    """
49    from math import sqrt, atan, pi
50
51    x = array(x)
52    y = array(y)
53
54    N = len(x)
55    a = 0*#New array
56
57    for k in range(N):
58        r = sqrt(x[k]**2 + y[k]**2)
59
60        angle = atan(y[k]/x[k])
61
62        if x[k] < 0:
63            angle+=pi #atan in ]-pi/2; pi/2[
64
65        #Take normal direction
66        angle -= pi/2
67
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
106    #FIXME (Ole): Individual flux tests do NOT test C implementation directly.   
107    def test_flux_zero_case(self):
108        ql = zeros( 3, Float )
109        qr = zeros( 3, Float )
110        normal = zeros( 2, Float )
111        zl = zr = 0.
112        flux, max_speed = flux_function(normal, ql, qr, zl, zr)
113
114        assert allclose(flux, [0,0,0])
115        assert max_speed == 0.
116
117    def test_flux_constants(self):
118        w = 2.0
119
120        normal = array([1.,0])
121        ql = array([w, 0, 0])
122        qr = array([w, 0, 0])
123        zl = zr = 0.
124        h = w - (zl+zr)/2
125
126        flux, max_speed = flux_function(normal, ql, qr, zl, zr)
127
128        assert allclose(flux, [0., 0.5*g*h**2, 0.])
129        assert max_speed == sqrt(g*h)
130
131    #def test_flux_slope(self):
132    #    #FIXME: TODO
133    #    w = 2.0
134    #
135    #    normal = array([1.,0])
136    #    ql = array([w, 0, 0])
137    #    qr = array([w, 0, 0])
138    #    zl = zr = 0.
139    #    h = w - (zl+zr)/2
140    #
141    #    flux, max_speed = flux_function(normal, ql, qr, zl, zr)
142    #
143    #    assert allclose(flux, [0., 0.5*g*h**2, 0.])
144    #    assert max_speed == sqrt(g*h)
145
146
147    def test_flux1(self):
148        #Use data from previous version of abstract_2d_finite_volumes
149        normal = array([1.,0])
150        ql = array([-0.2, 2, 3])
151        qr = array([-0.2, 2, 3])
152        zl = zr = -0.5
153        flux, max_speed = flux_function(normal, ql, qr, zl, zr)
154
155        assert allclose(flux, [2.,13.77433333, 20.])
156        assert allclose(max_speed, 8.38130948661)
157
158
159    def test_flux2(self):
160        #Use data from previous version of abstract_2d_finite_volumes
161        normal = array([0., -1.])
162        ql = array([-0.075, 2, 3])
163        qr = array([-0.075, 2, 3])
164        zl = zr = -0.375
165        flux, max_speed = flux_function(normal, ql, qr, zl, zr)
166
167        assert allclose(flux, [-3.,-20.0, -30.441])
168        assert allclose(max_speed, 11.7146428199)
169
170    def test_flux3(self):
171        #Use data from previous version of abstract_2d_finite_volumes
172        normal = array([-sqrt(2)/2, sqrt(2)/2])
173        ql = array([-0.075, 2, 3])
174        qr = array([-0.075, 2, 3])
175        zl = zr = -0.375
176        flux, max_speed = flux_function(normal, ql, qr, zl, zr)
177
178        assert allclose(flux, [sqrt(2)/2, 4.40221112, 7.3829019])
179        assert allclose(max_speed, 4.0716654239)
180
181    def test_flux4(self):
182        #Use data from previous version of abstract_2d_finite_volumes
183        normal = array([-sqrt(2)/2, sqrt(2)/2])
184        ql = array([-0.34319278, 0.10254161, 0.07273855])
185        qr = array([-0.30683287, 0.1071986, 0.05930515])
186        zl = zr = -0.375
187        flux, max_speed = flux_function(normal, ql, qr, zl, zr)
188
189        assert allclose(flux, [-0.04072676, -0.07096636, -0.01604364])
190        assert allclose(max_speed, 1.31414103233)
191
192    def test_sw_domain_simple(self):
193        a = [0.0, 0.0]
194        b = [0.0, 2.0]
195        c = [2.0,0.0]
196        d = [0.0, 4.0]
197        e = [2.0, 2.0]
198        f = [4.0,0.0]
199
200        points = [a, b, c, d, e, f]
201        #bac, bce, ecf, dbe, daf, dae
202        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
203
204
205        #from anuga.abstract_2d_finite_volumes.domain import Domain as Generic_domain
206        #msg = 'The class %s is not a subclass of the generic domain class %s'\
207        #      %(DomainClass, Domain)
208        #assert issubclass(DomainClass, Domain), msg
209
210        domain = Domain(points, vertices)
211        domain.check_integrity()
212
213        for name in ['stage', 'xmomentum', 'ymomentum',
214                     'elevation', 'friction']:
215            assert domain.quantities.has_key(name)
216
217
218        assert domain.get_conserved_quantities(0, edge=1) == 0.
219
220
221    def test_boundary_conditions(self):
222
223        a = [0.0, 0.0]
224        b = [0.0, 2.0]
225        c = [2.0,0.0]
226        d = [0.0, 4.0]
227        e = [2.0, 2.0]
228        f = [4.0,0.0]
229
230        points = [a, b, c, d, e, f]
231        #bac, bce, ecf, dbe
232        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4] ]
233        boundary = { (0, 0): 'Third',
234                     (0, 2): 'First',
235                     (2, 0): 'Second',
236                     (2, 1): 'Second',
237                     (3, 1): 'Second',
238                     (3, 2): 'Third'}
239
240
241        domain = Domain(points, vertices, boundary)
242        domain.check_integrity()
243
244
245        domain.set_quantity('stage', [[1,2,3], [5,5,5],
246                                      [0,0,9], [-6, 3, 3]])
247
248        domain.set_quantity('xmomentum', [[1,1,1], [2,2,2],
249                                          [3,3,3], [4, 4, 4]])
250
251        domain.set_quantity('ymomentum', [[10,10,10], [20,20,20],
252                                          [30,30,30], [40, 40, 40]])
253
254
255        D = Dirichlet_boundary([5,2,1])
256        T = Transmissive_boundary(domain)
257        R = Reflective_boundary(domain)
258        domain.set_boundary( {'First': D, 'Second': T, 'Third': R})
259
260        domain.update_boundary()
261
262        #Stage
263        assert domain.quantities['stage'].boundary_values[0] == 2.5
264        assert domain.quantities['stage'].boundary_values[0] ==\
265               domain.get_conserved_quantities(0, edge=0)[0] #Reflective (2.5)
266        assert domain.quantities['stage'].boundary_values[1] == 5. #Dirichlet
267        assert domain.quantities['stage'].boundary_values[2] ==\
268               domain.get_conserved_quantities(2, edge=0)[0] #Transmissive (4.5)
269        assert domain.quantities['stage'].boundary_values[3] ==\
270               domain.get_conserved_quantities(2, edge=1)[0] #Transmissive (4.5)
271        assert domain.quantities['stage'].boundary_values[4] ==\
272               domain.get_conserved_quantities(3, edge=1)[0] #Transmissive (-1.5)
273        assert domain.quantities['stage'].boundary_values[5] ==\
274               domain.get_conserved_quantities(3, edge=2)[0] #Reflective (-1.5)
275
276        #Xmomentum
277        assert domain.quantities['xmomentum'].boundary_values[0] == 1.0 #Reflective
278        assert domain.quantities['xmomentum'].boundary_values[1] == 2. #Dirichlet
279        assert domain.quantities['xmomentum'].boundary_values[2] ==\
280               domain.get_conserved_quantities(2, edge=0)[1] #Transmissive
281        assert domain.quantities['xmomentum'].boundary_values[3] ==\
282               domain.get_conserved_quantities(2, edge=1)[1] #Transmissive
283        assert domain.quantities['xmomentum'].boundary_values[4] ==\
284               domain.get_conserved_quantities(3, edge=1)[1] #Transmissive
285        assert domain.quantities['xmomentum'].boundary_values[5] == -4.0  #Reflective
286
287        #Ymomentum
288        assert domain.quantities['ymomentum'].boundary_values[0] == -10.0 #Reflective
289        assert domain.quantities['ymomentum'].boundary_values[1] == 1.  #Dirichlet
290        assert domain.quantities['ymomentum'].boundary_values[2] == 30. #Transmissive
291        assert domain.quantities['ymomentum'].boundary_values[3] == 30. #Transmissive
292        assert domain.quantities['ymomentum'].boundary_values[4] == 40. #Transmissive
293        assert domain.quantities['ymomentum'].boundary_values[5] == 40. #Reflective
294
295
296    def test_boundary_conditionsII(self):
297
298        a = [0.0, 0.0]
299        b = [0.0, 2.0]
300        c = [2.0,0.0]
301        d = [0.0, 4.0]
302        e = [2.0, 2.0]
303        f = [4.0,0.0]
304
305        points = [a, b, c, d, e, f]
306        #bac, bce, ecf, dbe
307        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4] ]
308        boundary = { (0, 0): 'Third',
309                     (0, 2): 'First',
310                     (2, 0): 'Second',
311                     (2, 1): 'Second',
312                     (3, 1): 'Second',
313                     (3, 2): 'Third',
314                     (0, 1): 'Internal'}
315
316
317        domain = Domain(points, vertices, boundary)
318        domain.check_integrity()
319
320
321        domain.set_quantity('stage', [[1,2,3], [5,5,5],
322                                      [0,0,9], [-6, 3, 3]])
323
324        domain.set_quantity('xmomentum', [[1,1,1], [2,2,2],
325                                          [3,3,3], [4, 4, 4]])
326
327        domain.set_quantity('ymomentum', [[10,10,10], [20,20,20],
328                                          [30,30,30], [40, 40, 40]])
329
330
331        D = Dirichlet_boundary([5,2,1])
332        T = Transmissive_boundary(domain)
333        R = Reflective_boundary(domain)
334        domain.set_boundary( {'First': D, 'Second': T,
335                              'Third': R, 'Internal': None})
336
337        domain.update_boundary()
338        domain.check_integrity()
339
340
341    def test_compute_fluxes0(self):
342        #Do a full triangle and check that fluxes cancel out for
343        #the constant stage case
344
345        a = [0.0, 0.0]
346        b = [0.0, 2.0]
347        c = [2.0,0.0]
348        d = [0.0, 4.0]
349        e = [2.0, 2.0]
350        f = [4.0,0.0]
351
352        points = [a, b, c, d, e, f]
353        #bac, bce, ecf, dbe
354        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
355
356        domain = Domain(points, vertices)
357        domain.set_quantity('stage', [[2,2,2], [2,2,2],
358                                      [2,2,2], [2,2,2]])
359        domain.check_integrity()
360
361        assert allclose(domain.neighbours, [[-1,1,-1], [2,3,0], [-1,-1,1],[1,-1,-1]])
362        assert allclose(domain.neighbour_edges, [[-1,2,-1], [2,0,1], [-1,-1,0],[1,-1,-1]])
363
364        zl=zr=0. #Assume flat bed
365
366        #Flux across right edge of volume 1
367        normal = domain.get_normal(1,0)
368        ql = domain.get_conserved_quantities(vol_id=1, edge=0)
369        qr = domain.get_conserved_quantities(vol_id=2, edge=2)
370        flux0, max_speed = flux_function(normal, ql, qr, zl, zr)
371
372        #Check that flux seen from other triangles is inverse
373        tmp = qr; qr=ql; ql=tmp
374        normal = domain.get_normal(2,2)
375        flux, max_speed = flux_function(normal, ql, qr, zl, zr)
376        assert allclose(flux + flux0, 0.)
377
378        #Flux across upper edge of volume 1
379        normal = domain.get_normal(1,1)
380        ql = domain.get_conserved_quantities(vol_id=1, edge=1)
381        qr = domain.get_conserved_quantities(vol_id=3, edge=0)
382        flux1, max_speed = flux_function(normal, ql, qr, zl, zr)
383
384        #Check that flux seen from other triangles is inverse
385        tmp = qr; qr=ql; ql=tmp
386        normal = domain.get_normal(3,0)
387        flux, max_speed = flux_function(normal, ql, qr, zl, zr)
388        assert allclose(flux + flux1, 0.)
389
390        #Flux across lower left hypotenuse of volume 1
391        normal = domain.get_normal(1,2)
392        ql = domain.get_conserved_quantities(vol_id=1, edge=2)
393        qr = domain.get_conserved_quantities(vol_id=0, edge=1)
394        flux2, max_speed = flux_function(normal, ql, qr, zl, zr)
395
396        #Check that flux seen from other triangles is inverse
397        tmp = qr; qr=ql; ql=tmp
398        normal = domain.get_normal(0,1)
399        flux, max_speed = flux_function(normal, ql, qr, zl, zr)
400        assert allclose(flux + flux2, 0.)
401
402
403        #Scale by edgelengths, add up anc check that total flux is zero
404        e0 = domain.edgelengths[1, 0]
405        e1 = domain.edgelengths[1, 1]
406        e2 = domain.edgelengths[1, 2]
407
408        assert allclose(e0*flux0+e1*flux1+e2*flux2, 0.)
409
410        #Now check that compute_flux yields zeros as well
411        domain.compute_fluxes()
412
413        for name in ['stage', 'xmomentum', 'ymomentum']:
414            #print name, domain.quantities[name].explicit_update
415            assert allclose(domain.quantities[name].explicit_update[1], 0)
416
417
418
419    def test_compute_fluxes1(self):
420        #Use values from previous version
421
422        a = [0.0, 0.0]
423        b = [0.0, 2.0]
424        c = [2.0,0.0]
425        d = [0.0, 4.0]
426        e = [2.0, 2.0]
427        f = [4.0,0.0]
428
429        points = [a, b, c, d, e, f]
430        #bac, bce, ecf, dbe
431        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
432
433        domain = Domain(points, vertices)
434        val0 = 2.+2.0/3
435        val1 = 4.+4.0/3
436        val2 = 8.+2.0/3
437        val3 = 2.+8.0/3
438
439        domain.set_quantity('stage', [[val0, val0, val0], [val1, val1, val1],
440                                      [val2, val2, val2], [val3, val3, val3]])
441        domain.check_integrity()
442
443        zl=zr=0. #Assume flat bed
444
445        #Flux across right edge of volume 1
446        normal = domain.get_normal(1,0) #Get normal 0 of triangle 1
447        assert allclose(normal, [1, 0])
448       
449        ql = domain.get_conserved_quantities(vol_id=1, edge=0)
450        assert allclose(ql, [val1, 0, 0])
451       
452        qr = domain.get_conserved_quantities(vol_id=2, edge=2)
453        assert allclose(qr, [val2, 0, 0])
454       
455        flux0, max_speed = flux_function(normal, ql, qr, zl, zr)
456
457        #Flux across edge in the east direction (as per normal vector)
458        assert allclose(flux0, [-15.3598804, 253.71111111, 0.])
459        assert allclose(max_speed, 9.21592824046)
460
461
462        #Flux across edge in the west direction (opposite sign for xmomentum)
463        normal_opposite = domain.get_normal(2,2) #Get normal 2 of triangle 2
464        assert allclose(normal_opposite, [-1, 0])
465        flux_opposite, max_speed = flux_function([-1, 0], ql, qr, zl, zr)
466        assert allclose(flux_opposite, [-15.3598804, -253.71111111, 0.])
467       
468
469        #Flux across upper edge of volume 1
470        normal = domain.get_normal(1,1)
471        ql = domain.get_conserved_quantities(vol_id=1, edge=1)
472        qr = domain.get_conserved_quantities(vol_id=3, edge=0)
473        flux1, max_speed = flux_function(normal, ql, qr, zl, zr)
474        assert allclose(flux1, [2.4098563, 0., 123.04444444])
475        assert allclose(max_speed, 7.22956891292)
476
477        #Flux across lower left hypotenuse of volume 1
478        normal = domain.get_normal(1,2)
479        ql = domain.get_conserved_quantities(vol_id=1, edge=2)
480        qr = domain.get_conserved_quantities(vol_id=0, edge=1)
481        flux2, max_speed = flux_function(normal, ql, qr, zl, zr)
482
483        assert allclose(flux2, [9.63942522, -61.59685738, -61.59685738])
484        assert allclose(max_speed, 7.22956891292)
485
486        #Scale, add up and check that compute_fluxes is correct for vol 1
487        e0 = domain.edgelengths[1, 0]
488        e1 = domain.edgelengths[1, 1]
489        e2 = domain.edgelengths[1, 2]
490
491        total_flux = -(e0*flux0+e1*flux1+e2*flux2)/domain.areas[1]
492        assert allclose(total_flux, [-0.68218178, -166.6, -35.93333333])
493
494
495        domain.compute_fluxes()
496
497        #assert allclose(total_flux, domain.explicit_update[1,:])
498        for i, name in enumerate(['stage', 'xmomentum', 'ymomentum']):
499            assert allclose(total_flux[i],
500                            domain.quantities[name].explicit_update[1])
501
502        #assert allclose(domain.explicit_update, [
503        #    [0., -69.68888889, -69.68888889],
504        #    [-0.68218178, -166.6, -35.93333333],
505        #    [-111.77316251, 69.68888889, 0.],
506        #    [-35.68522449, 0., 69.68888889]])
507
508        assert allclose(domain.quantities['stage'].explicit_update,
509                        [0., -0.68218178, -111.77316251, -35.68522449])
510        assert allclose(domain.quantities['xmomentum'].explicit_update,
511                        [-69.68888889, -166.6, 69.68888889, 0])
512        assert allclose(domain.quantities['ymomentum'].explicit_update,
513                        [-69.68888889, -35.93333333, 0., 69.68888889])
514
515
516        #assert allclose(domain.quantities[name].explicit_update
517
518
519
520
521
522    def test_compute_fluxes2(self):
523        #Random values, incl momentum
524
525        a = [0.0, 0.0]
526        b = [0.0, 2.0]
527        c = [2.0,0.0]
528        d = [0.0, 4.0]
529        e = [2.0, 2.0]
530        f = [4.0,0.0]
531
532        points = [a, b, c, d, e, f]
533        #bac, bce, ecf, dbe
534        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
535
536        domain = Domain(points, vertices)
537        val0 = 2.+2.0/3
538        val1 = 4.+4.0/3
539        val2 = 8.+2.0/3
540        val3 = 2.+8.0/3
541
542        zl=zr=0 #Assume flat zero bed
543
544        domain.set_quantity('elevation', zl*ones( (4,3) ))
545
546
547        domain.set_quantity('stage', [[val0, val0-1, val0-2],
548                                      [val1, val1+1, val1],
549                                      [val2, val2-2, val2],
550                                      [val3-0.5, val3, val3]])
551
552        domain.set_quantity('xmomentum',
553                            [[1, 2, 3], [3, 4, 5],
554                             [1, -1, 0], [0, -2, 2]])
555
556        domain.set_quantity('ymomentum',
557                            [[1, -1, 0], [0, -3, 2],
558                             [0, 1, 0], [-1, 2, 2]])
559
560
561        domain.check_integrity()
562
563
564
565        #Flux across right edge of volume 1
566        normal = domain.get_normal(1,0)
567        ql = domain.get_conserved_quantities(vol_id=1, edge=0)
568        qr = domain.get_conserved_quantities(vol_id=2, edge=2)
569        flux0, max_speed = flux_function(normal, ql, qr, zl, zr)
570
571        #Flux across upper edge of volume 1
572        normal = domain.get_normal(1,1)
573        ql = domain.get_conserved_quantities(vol_id=1, edge=1)
574        qr = domain.get_conserved_quantities(vol_id=3, edge=0)
575        flux1, max_speed = flux_function(normal, ql, qr, zl, zr)
576
577        #Flux across lower left hypotenuse of volume 1
578        normal = domain.get_normal(1,2)
579        ql = domain.get_conserved_quantities(vol_id=1, edge=2)
580        qr = domain.get_conserved_quantities(vol_id=0, edge=1)
581        flux2, max_speed = flux_function(normal, ql, qr, zl, zr)
582
583        #Scale, add up and check that compute_fluxes is correct for vol 1
584        e0 = domain.edgelengths[1, 0]
585        e1 = domain.edgelengths[1, 1]
586        e2 = domain.edgelengths[1, 2]
587
588        total_flux = -(e0*flux0+e1*flux1+e2*flux2)/domain.areas[1]
589
590
591        domain.compute_fluxes()
592        for i, name in enumerate(['stage', 'xmomentum', 'ymomentum']):
593            assert allclose(total_flux[i],
594                            domain.quantities[name].explicit_update[1])
595        #assert allclose(total_flux, domain.explicit_update[1,:])
596
597
598    # FIXME (Ole): Need test like this for fluxes in very shallow water.   
599    def test_compute_fluxes3(self):
600        #Random values, incl momentum
601
602        a = [0.0, 0.0]
603        b = [0.0, 2.0]
604        c = [2.0,0.0]
605        d = [0.0, 4.0]
606        e = [2.0, 2.0]
607        f = [4.0,0.0]
608
609        points = [a, b, c, d, e, f]
610        #bac, bce, ecf, dbe
611        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
612
613        domain = Domain(points, vertices)
614        val0 = 2.+2.0/3
615        val1 = 4.+4.0/3
616        val2 = 8.+2.0/3
617        val3 = 2.+8.0/3
618
619        zl=zr=-3.75 #Assume constant bed (must be less than stage)
620        domain.set_quantity('elevation', zl*ones( (4,3) ))
621
622
623        domain.set_quantity('stage', [[val0, val0-1, val0-2],
624                                      [val1, val1+1, val1],
625                                      [val2, val2-2, val2],
626                                      [val3-0.5, val3, val3]])
627
628        domain.set_quantity('xmomentum',
629                            [[1, 2, 3], [3, 4, 5],
630                             [1, -1, 0], [0, -2, 2]])
631
632        domain.set_quantity('ymomentum',
633                            [[1, -1, 0], [0, -3, 2],
634                             [0, 1, 0], [-1, 2, 2]])
635
636
637        domain.check_integrity()
638
639
640
641        #Flux across right edge of volume 1
642        normal = domain.get_normal(1,0)
643        ql = domain.get_conserved_quantities(vol_id=1, edge=0)
644        qr = domain.get_conserved_quantities(vol_id=2, edge=2)
645        flux0, max_speed = flux_function(normal, ql, qr, zl, zr)
646
647        #Flux across upper edge of volume 1
648        normal = domain.get_normal(1,1)
649        ql = domain.get_conserved_quantities(vol_id=1, edge=1)
650        qr = domain.get_conserved_quantities(vol_id=3, edge=0)
651        flux1, max_speed = flux_function(normal, ql, qr, zl, zr)
652
653        #Flux across lower left hypotenuse of volume 1
654        normal = domain.get_normal(1,2)
655        ql = domain.get_conserved_quantities(vol_id=1, edge=2)
656        qr = domain.get_conserved_quantities(vol_id=0, edge=1)
657        flux2, max_speed = flux_function(normal, ql, qr, zl, zr)
658
659        #Scale, add up and check that compute_fluxes is correct for vol 1
660        e0 = domain.edgelengths[1, 0]
661        e1 = domain.edgelengths[1, 1]
662        e2 = domain.edgelengths[1, 2]
663
664        total_flux = -(e0*flux0+e1*flux1+e2*flux2)/domain.areas[1]
665
666        domain.compute_fluxes()
667        for i, name in enumerate(['stage', 'xmomentum', 'ymomentum']):
668            assert allclose(total_flux[i],
669                            domain.quantities[name].explicit_update[1])
670
671
672
673    def xtest_catching_negative_heights(self):
674
675        #OBSOLETE
676       
677        a = [0.0, 0.0]
678        b = [0.0, 2.0]
679        c = [2.0,0.0]
680        d = [0.0, 4.0]
681        e = [2.0, 2.0]
682        f = [4.0,0.0]
683
684        points = [a, b, c, d, e, f]
685        #bac, bce, ecf, dbe
686        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
687
688        domain = Domain(points, vertices)
689        val0 = 2.+2.0/3
690        val1 = 4.+4.0/3
691        val2 = 8.+2.0/3
692        val3 = 2.+8.0/3
693
694        zl=zr=4  #Too large
695        domain.set_quantity('elevation', zl*ones( (4,3) ))
696        domain.set_quantity('stage', [[val0, val0-1, val0-2],
697                                      [val1, val1+1, val1],
698                                      [val2, val2-2, val2],
699                                      [val3-0.5, val3, val3]])
700
701        #Should fail
702        try:
703            domain.check_integrity()
704        except:
705            pass
706
707
708
709    def test_get_wet_elements(self):
710
711        a = [0.0, 0.0]
712        b = [0.0, 2.0]
713        c = [2.0,0.0]
714        d = [0.0, 4.0]
715        e = [2.0, 2.0]
716        f = [4.0,0.0]
717
718        points = [a, b, c, d, e, f]
719        #bac, bce, ecf, dbe
720        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
721
722        domain = Domain(points, vertices)
723        val0 = 2.+2.0/3
724        val1 = 4.+4.0/3
725        val2 = 8.+2.0/3
726        val3 = 2.+8.0/3
727
728        zl=zr=5
729        domain.set_quantity('elevation', zl*ones( (4,3) ))
730        domain.set_quantity('stage', [[val0, val0-1, val0-2],
731                                      [val1, val1+1, val1],
732                                      [val2, val2-2, val2],
733                                      [val3-0.5, val3, val3]])
734
735
736
737        #print domain.get_quantity('elevation').get_values(location='centroids')
738        #print domain.get_quantity('stage').get_values(location='centroids')       
739        domain.check_integrity()
740
741        indices = domain.get_wet_elements()
742        assert allclose(indices, [1,2])
743
744        indices = domain.get_wet_elements(indices=[0,1,3])
745        assert allclose(indices, [1])
746       
747
748
749    def test_get_maximum_inundation_1(self):
750
751        a = [0.0, 0.0]
752        b = [0.0, 2.0]
753        c = [2.0,0.0]
754        d = [0.0, 4.0]
755        e = [2.0, 2.0]
756        f = [4.0,0.0]
757
758        points = [a, b, c, d, e, f]
759        #bac, bce, ecf, dbe
760        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
761
762        domain = Domain(points, vertices)
763
764        domain.set_quantity('elevation', lambda x, y: x+2*y) #2 4 4 6
765        domain.set_quantity('stage', 3)
766
767        domain.check_integrity()
768
769        indices = domain.get_wet_elements()
770        assert allclose(indices, [0])
771
772        q = domain.get_maximum_inundation_elevation()
773        assert allclose(q, domain.get_quantity('elevation').get_values(location='centroids')[0])
774
775        x, y = domain.get_maximum_inundation_location()
776        assert allclose([x, y], domain.get_centroid_coordinates()[0])
777
778
779    def test_get_maximum_inundation_2(self):
780        """test_get_maximum_inundation_2(self)
781
782        Test multiple wet cells with same elevation
783        """
784       
785        a = [0.0, 0.0]
786        b = [0.0, 2.0]
787        c = [2.0,0.0]
788        d = [0.0, 4.0]
789        e = [2.0, 2.0]
790        f = [4.0,0.0]
791
792        points = [a, b, c, d, e, f]
793        #bac, bce, ecf, dbe
794        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
795
796        domain = Domain(points, vertices)
797
798        domain.set_quantity('elevation', lambda x, y: x+2*y) #2 4 4 6
799        domain.set_quantity('stage', 4.1)
800
801        domain.check_integrity()
802
803        indices = domain.get_wet_elements()
804        assert allclose(indices, [0,1,2])
805
806        q = domain.get_maximum_inundation_elevation()
807        assert allclose(q, 4)       
808
809        x, y = domain.get_maximum_inundation_location()
810        assert allclose([x, y], domain.get_centroid_coordinates()[1])       
811       
812
813    def test_get_maximum_inundation_3(self):
814        """test_get_maximum_inundation_3(self)
815
816        Test real runup example
817        """
818
819        from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular_cross
820
821        initial_runup_height = -0.4
822        final_runup_height = -0.3
823
824
825        #--------------------------------------------------------------
826        # Setup computational domain
827        #--------------------------------------------------------------
828        N = 5
829        points, vertices, boundary = rectangular_cross(N, N) 
830        domain = Domain(points, vertices, boundary)           
831
832        #--------------------------------------------------------------
833        # Setup initial conditions
834        #--------------------------------------------------------------
835        def topography(x,y):
836            return -x/2                             # linear bed slope
837           
838
839        domain.set_quantity('elevation', topography)       # Use function for elevation
840        domain.set_quantity('friction', 0.)                # Zero friction
841        domain.set_quantity('stage', initial_runup_height) # Constant negative initial stage
842
843
844        #--------------------------------------------------------------
845        # Setup boundary conditions
846        #--------------------------------------------------------------
847        Br = Reflective_boundary(domain)              # Reflective wall
848        Bd = Dirichlet_boundary([final_runup_height,  # Constant inflow
849                                 0,
850                                 0])
851
852        # All reflective to begin with (still water)
853        domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
854
855
856        #--------------------------------------------------------------
857        # Test initial inundation height
858        #--------------------------------------------------------------
859
860        indices = domain.get_wet_elements()
861        z = domain.get_quantity('elevation').\
862            get_values(location='centroids', indices=indices)
863        assert alltrue(z<initial_runup_height)
864
865        q = domain.get_maximum_inundation_elevation()
866        assert allclose(q, initial_runup_height, rtol = 1.0/N) # First order accuracy
867
868        x, y = domain.get_maximum_inundation_location()
869
870        qref = domain.get_quantity('elevation').get_values(interpolation_points = [[x, y]])
871        assert allclose(q, qref)
872
873
874        wet_elements = domain.get_wet_elements()
875        wet_elevations = domain.get_quantity('elevation').get_values(location='centroids',
876                                                                     indices=wet_elements)
877        assert alltrue(wet_elevations<initial_runup_height)
878        assert allclose(wet_elevations, z)       
879
880
881        #print domain.get_quantity('elevation').get_maximum_value(indices=wet_elements)
882        #print domain.get_quantity('elevation').get_maximum_location(indices=wet_elements)
883        #print domain.get_quantity('elevation').get_maximum_index(indices=wet_elements)
884
885       
886        #--------------------------------------------------------------
887        # Let triangles adjust
888        #--------------------------------------------------------------
889        for t in domain.evolve(yieldstep = 0.1, finaltime = 1.0):
890            pass
891
892
893        #--------------------------------------------------------------
894        # Test inundation height again
895        #--------------------------------------------------------------
896
897        indices = domain.get_wet_elements()
898        z = domain.get_quantity('elevation').\
899            get_values(location='centroids', indices=indices)
900
901        assert alltrue(z<initial_runup_height)
902
903        q = domain.get_maximum_inundation_elevation()
904        assert allclose(q, initial_runup_height, rtol = 1.0/N) # First order accuracy
905       
906        x, y = domain.get_maximum_inundation_location()
907        qref = domain.get_quantity('elevation').get_values(interpolation_points = [[x, y]])
908        assert allclose(q, qref)       
909
910
911        #--------------------------------------------------------------
912        # Update boundary to allow inflow
913        #--------------------------------------------------------------
914        domain.set_boundary({'right': Bd})
915
916       
917        #--------------------------------------------------------------
918        # Evolve system through time
919        #--------------------------------------------------------------
920        for t in domain.evolve(yieldstep = 0.1, finaltime = 3.0):
921            #domain.write_time()
922            pass
923   
924        #--------------------------------------------------------------
925        # Test inundation height again
926        #--------------------------------------------------------------
927
928        indices = domain.get_wet_elements()
929        z = domain.get_quantity('elevation').\
930            get_values(location='centroids', indices=indices)
931
932        assert alltrue(z<final_runup_height)
933
934        q = domain.get_maximum_inundation_elevation()
935        assert allclose(q, final_runup_height, rtol = 1.0/N) # First order accuracy
936
937        x, y = domain.get_maximum_inundation_location()
938        qref = domain.get_quantity('elevation').get_values(interpolation_points = [[x, y]])
939        assert allclose(q, qref)       
940
941
942        wet_elements = domain.get_wet_elements()
943        wet_elevations = domain.get_quantity('elevation').get_values(location='centroids',
944                                                                     indices=wet_elements)
945        assert alltrue(wet_elevations<final_runup_height)
946        assert allclose(wet_elevations, z)       
947       
948
949
950    def test_another_runup_example(self):
951        """test_another_runup_example
952
953        Test runup example where actual timeseries at interpolated
954        points are tested.
955        """
956
957        #-----------------------------------------------------------------
958        # Import necessary modules
959        #-----------------------------------------------------------------
960
961        from anuga.pmesh.mesh_interface import create_mesh_from_regions
962        from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular_cross
963        from anuga.shallow_water import Domain
964        from anuga.shallow_water import Reflective_boundary
965        from anuga.shallow_water import Dirichlet_boundary
966
967
968        #-----------------------------------------------------------------
969        # Setup computational domain
970        #-----------------------------------------------------------------
971        points, vertices, boundary = rectangular_cross(10, 10) # Basic mesh
972        domain = Domain(points, vertices, boundary) # Create domain
973        domain.set_quantities_to_be_stored(None)
974        domain.set_maximum_allowed_speed(100) #
975        domain.H0 = 0 # Backwards compatibility (6/2/7)
976        domain.beta_h = 0.2 # Backwards compatibility (14/2/7)
977
978        #-----------------------------------------------------------------
979        # Setup initial conditions
980        #-----------------------------------------------------------------
981
982        def topography(x,y): 
983            return -x/2                              # linear bed slope
984
985        domain.set_quantity('elevation', topography) 
986        domain.set_quantity('friction', 0.0)         
987        domain.set_quantity('stage', expression='elevation')           
988
989
990        #----------------------------------------------------------------
991        # Setup boundary conditions
992        #----------------------------------------------------------------
993
994        Br = Reflective_boundary(domain)      # Solid reflective wall
995        Bd = Dirichlet_boundary([-0.2,0.,0.]) # Constant boundary values
996        domain.set_boundary({'left': Br, 'right': Bd, 'top': Br, 'bottom': Br})
997
998
999        #----------------------------------------------------------------
1000        # Evolve system through time
1001        #----------------------------------------------------------------
1002
1003        interpolation_points = [[0.4,0.5], [0.6,0.5], [0.8,0.5], [0.9,0.5]]
1004        gauge_values = []
1005        for _ in interpolation_points:
1006            gauge_values.append([])
1007
1008        time = []
1009        for t in domain.evolve(yieldstep = 0.1, finaltime = 5.0):
1010            # Record time series at known points
1011            time.append(domain.get_time())
1012           
1013            stage = domain.get_quantity('stage')
1014            w = stage.get_values(interpolation_points=interpolation_points)
1015           
1016            for i, _ in enumerate(interpolation_points):
1017                gauge_values[i].append(w[i])
1018
1019
1020        #print
1021        #print time
1022        #print
1023        #for i, (x,y) in enumerate(interpolation_points):
1024        #    print i, gauge_values[i]
1025        #    print
1026
1027        #Reference (nautilus 13/10/2006)
1028
1029        G0 = [-0.20000000000000001, -0.19999681443389281, -0.1986192343695303, -0.19147413648863046, -0.19132688908678019, -0.17642317476621105, -0.167376262630034, -0.16192452887426961, -0.15609171725778803, -0.15127107084302249, -0.14048864340360018, -0.19296484125327093, -0.19997006390580363, -0.19999999999937063, -0.19999999999937063, -0.19999999999938772, -0.19999999999938772, -0.19999999999938772, -0.19999999999938772, -0.19974288463035494, -0.19951636867991712, -0.19966301435195755, -0.19981082259800226, -0.19978575003960128, -0.19992942471933109, -0.19999999931029933, -0.19999999999906989, -0.19999999999906989, -0.19999999999906989, -0.19999999999906989, -0.19999999999906989, -0.19999999999906989, -0.19999999999906989, -0.19999999999906989, -0.19999999999906989, -0.19999999999906989, -0.19999999999906989, -0.19999999999906989, -0.19999999999906989, -0.19999999999906989, -0.19999999999906989, -0.19999999999906989, -0.19999999999906989, -0.19999999999906989, -0.19999999999906989, -0.19999999999906989, -0.19999999999906989, -0.19999999999906989, -0.19999999999906989, -0.19999999999906989, -0.19999999999906989]
1030
1031
1032        G1 = [-0.29999999999999993, -0.29988962537199199, -0.29293904425532025, -0.28329367722887888, -0.25999146407696289, -0.22613875068011896, -0.21190705052094994, -0.19900707995208217, -0.18876305176191882, -0.18132447501091936, -0.17395459512711151, -0.15562414200985644, -0.16212999953643359, -0.18964422820514618, -0.20871181844346975, -0.21672207791083464, -0.21774940291862779, -0.21482868050219833, -0.21057786776704043, -0.20649663432591045, -0.20294932949211578, -0.19974459897911329, -0.19733648772704043, -0.19641404599824669, -0.19654095699184146, -0.19709942852191994, -0.19780873983410741, -0.19853259125123518, -0.19916495938961168, -0.19965391267799168, -0.19993539587158982, -0.2001383705551133, -0.20029344332295113, -0.20035349748150011, -0.20029886541561631, -0.20015541958920294, -0.19997273066429103, -0.19979879448668514, -0.19966016997024041, -0.19957558009501869, -0.19955725674938532, -0.19958083002853366, -0.19961752462568647, -0.19965296611330258, -0.19968998132634594, -0.19972532942208607, -0.19975372922008239, -0.19977196116929855, -0.19977951443660594, -0.19977792107284789, -0.19976991595502003]
1033
1034        G2 = [-0.40000000000000002, -0.39011996186687281, -0.33359026016903887, -0.29757449757405952, -0.27594124995715791, -0.25970211955309436, -0.24482929492054245, -0.23156757139219822, -0.21956485769139392, -0.20844522129026694, -0.19856327660654355, -0.18962303467030903, -0.17371085465024955, -0.16429840256208336, -0.17793711732368575, -0.19287799702389993, -0.20236271260796762, -0.20700727993623128, -0.20847704371373174, -0.20796895600687262, -0.20653398626186478, -0.20480656169870676, -0.20295863990994492, -0.20100199602968896, -0.19940642689498472, -0.19858371478015749, -0.19838672154605322, -0.19851093923669558, -0.19878191998909323, -0.19910827645394291, -0.19943514333832094, -0.19971231361970535, -0.19992429278849655, -0.20010744405928019, -0.20025927002359642, -0.20034751667523681, -0.20035504591467249, -0.20029401385620157, -0.20019492358237226, -0.20008934249434918, -0.19999808924091636, -0.19993869218976712, -0.19991589568150098, -0.19991815777945968, -0.19993012995477188, -0.19994576118144997, -0.19996497193815974, -0.19998586151236197, -0.20000487253824847, -0.20001903000364174, -0.20002698661385457]
1035
1036        G3 = [-0.45000000000000001, -0.37713945714588398, -0.33029565026933816, -0.30598209033945367, -0.28847101155177313, -0.27211191064563195, -0.25701544058818926, -0.24298945948410997, -0.23010402733784807, -0.21820351802867713, -0.20709938367218383, -0.19719881806182216, -0.18568281604361933, -0.16828653906676322, -0.16977310768235579, -0.1832707289594605, -0.19483524345250974, -0.20233480051649216, -0.20630757214159207, -0.20763927857964531, -0.20724458160595791, -0.20599191745446047, -0.20438329669495012, -0.20256105512496606, -0.20071269486729407, -0.19934403619901719, -0.19866860191898347, -0.19849975056296071, -0.19860870923007437, -0.19885838217851401, -0.19916422433758982, -0.19946861981642039, -0.19972267778871666, -0.19993013816258154, -0.20011063428833351, -0.20024891930311628, -0.20031882555219671, -0.20031326268593497, -0.20024881068472311, -0.20015443214902759, -0.20005669097631221, -0.19997542564643309, -0.19992564006223304, -0.19990746148869892, -0.19990923999172872, -0.19991956416813192, -0.19993484556273733, -0.1999538628054662, -0.19997381636620407, -0.19999130900268777, -0.20000388227457688]
1037       
1038        #FIXME (DSG):This is a hack so the anuga install, not precompiled
1039        # works on DSG's win2000, python 2.3
1040        #The problem is the gauge_values[X] are 52 long, not 51.
1041        #
1042        # This was probably fixed by Stephen in changeset:3804
1043        if len(gauge_values[0]) == 52: gauge_values[0].pop()
1044        if len(gauge_values[1]) == 52: gauge_values[1].pop()
1045        if len(gauge_values[2]) == 52: gauge_values[2].pop()
1046        if len(gauge_values[3]) == 52: gauge_values[3].pop()
1047
1048##         print len(G0), len(gauge_values[0])
1049##         print len(G1), len(gauge_values[1])
1050##         print gauge_values[0]
1051##         print G0
1052
1053       
1054       
1055        assert allclose(gauge_values[0], G0)
1056        assert allclose(gauge_values[1], G1)
1057        assert allclose(gauge_values[2], G2)
1058        assert allclose(gauge_values[3], G3)       
1059
1060
1061
1062
1063
1064
1065
1066    #####################################################
1067    def test_initial_condition(self):
1068        """Test that initial condition is output at time == 0
1069        """
1070
1071        from anuga.config import g
1072        import copy
1073
1074        a = [0.0, 0.0]
1075        b = [0.0, 2.0]
1076        c = [2.0, 0.0]
1077        d = [0.0, 4.0]
1078        e = [2.0, 2.0]
1079        f = [4.0, 0.0]
1080
1081        points = [a, b, c, d, e, f]
1082        #bac, bce, ecf, dbe
1083        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
1084
1085        domain = Domain(points, vertices)
1086
1087        #Set up for a gradient of (3,0) at mid triangle
1088        def slope(x, y):
1089            return 3*x
1090
1091        h = 0.1
1092        def stage(x,y):
1093            return slope(x,y)+h
1094
1095        domain.set_quantity('elevation', slope)
1096        domain.set_quantity('stage', stage)
1097
1098        initial_stage = copy.copy(domain.quantities['stage'].vertex_values)
1099
1100        domain.set_boundary({'exterior': Reflective_boundary(domain)})
1101
1102        #Evolution
1103        for t in domain.evolve(yieldstep = 1.0, finaltime = 2.0):
1104            stage = domain.quantities['stage'].vertex_values
1105
1106            if t == 0.0:
1107                assert allclose(stage, initial_stage)
1108            else:
1109                assert not allclose(stage, initial_stage)
1110
1111        os.remove(domain.get_name() + '.sww')
1112
1113
1114
1115    #####################################################
1116    def test_gravity(self):
1117        #Assuming no friction
1118
1119        from anuga.config import g
1120
1121        a = [0.0, 0.0]
1122        b = [0.0, 2.0]
1123        c = [2.0, 0.0]
1124        d = [0.0, 4.0]
1125        e = [2.0, 2.0]
1126        f = [4.0, 0.0]
1127
1128        points = [a, b, c, d, e, f]
1129        #bac, bce, ecf, dbe
1130        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
1131
1132        domain = Domain(points, vertices)
1133
1134        #Set up for a gradient of (3,0) at mid triangle
1135        def slope(x, y):
1136            return 3*x
1137
1138        h = 0.1
1139        def stage(x,y):
1140            return slope(x,y)+h
1141
1142        domain.set_quantity('elevation', slope)
1143        domain.set_quantity('stage', stage)
1144
1145        for name in domain.conserved_quantities:
1146            assert allclose(domain.quantities[name].explicit_update, 0)
1147            assert allclose(domain.quantities[name].semi_implicit_update, 0)
1148
1149        domain.compute_forcing_terms()
1150
1151        assert allclose(domain.quantities['stage'].explicit_update, 0)
1152        assert allclose(domain.quantities['xmomentum'].explicit_update, -g*h*3)
1153        assert allclose(domain.quantities['ymomentum'].explicit_update, 0)
1154
1155
1156    def test_manning_friction(self):
1157        from anuga.config import g
1158
1159        a = [0.0, 0.0]
1160        b = [0.0, 2.0]
1161        c = [2.0, 0.0]
1162        d = [0.0, 4.0]
1163        e = [2.0, 2.0]
1164        f = [4.0, 0.0]
1165
1166        points = [a, b, c, d, e, f]
1167        #bac, bce, ecf, dbe
1168        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
1169
1170        domain = Domain(points, vertices)
1171
1172        #Set up for a gradient of (3,0) at mid triangle
1173        def slope(x, y):
1174            return 3*x
1175
1176        h = 0.1
1177        def stage(x,y):
1178            return slope(x,y)+h
1179
1180        eta = 0.07
1181        domain.set_quantity('elevation', slope)
1182        domain.set_quantity('stage', stage)
1183        domain.set_quantity('friction', eta)
1184
1185        for name in domain.conserved_quantities:
1186            assert allclose(domain.quantities[name].explicit_update, 0)
1187            assert allclose(domain.quantities[name].semi_implicit_update, 0)
1188
1189        domain.compute_forcing_terms()
1190
1191        assert allclose(domain.quantities['stage'].explicit_update, 0)
1192        assert allclose(domain.quantities['xmomentum'].explicit_update, -g*h*3)
1193        assert allclose(domain.quantities['ymomentum'].explicit_update, 0)
1194
1195        assert allclose(domain.quantities['stage'].semi_implicit_update, 0)
1196        assert allclose(domain.quantities['xmomentum'].semi_implicit_update, 0)
1197        assert allclose(domain.quantities['ymomentum'].semi_implicit_update, 0)
1198
1199        #Create some momentum for friction to work with
1200        domain.set_quantity('xmomentum', 1)
1201        S = -g * eta**2 / h**(7.0/3)
1202
1203        domain.compute_forcing_terms()
1204        assert allclose(domain.quantities['stage'].semi_implicit_update, 0)
1205        assert allclose(domain.quantities['xmomentum'].semi_implicit_update, S)
1206        assert allclose(domain.quantities['ymomentum'].semi_implicit_update, 0)
1207
1208        #A more complex example
1209        domain.quantities['stage'].semi_implicit_update[:] = 0.0
1210        domain.quantities['xmomentum'].semi_implicit_update[:] = 0.0
1211        domain.quantities['ymomentum'].semi_implicit_update[:] = 0.0
1212
1213        domain.set_quantity('xmomentum', 3)
1214        domain.set_quantity('ymomentum', 4)
1215
1216        S = -g * eta**2 * 5 / h**(7.0/3)
1217
1218
1219        domain.compute_forcing_terms()
1220
1221        assert allclose(domain.quantities['stage'].semi_implicit_update, 0)
1222        assert allclose(domain.quantities['xmomentum'].semi_implicit_update, 3*S)
1223        assert allclose(domain.quantities['ymomentum'].semi_implicit_update, 4*S)
1224
1225    def test_constant_wind_stress(self):
1226        from anuga.config import rho_a, rho_w, eta_w
1227        from math import pi, cos, sin, sqrt
1228
1229        a = [0.0, 0.0]
1230        b = [0.0, 2.0]
1231        c = [2.0, 0.0]
1232        d = [0.0, 4.0]
1233        e = [2.0, 2.0]
1234        f = [4.0, 0.0]
1235
1236        points = [a, b, c, d, e, f]
1237        #bac, bce, ecf, dbe
1238        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
1239
1240
1241        domain = Domain(points, vertices)
1242
1243        #Flat surface with 1m of water
1244        domain.set_quantity('elevation', 0)
1245        domain.set_quantity('stage', 1.0)
1246        domain.set_quantity('friction', 0)
1247
1248        Br = Reflective_boundary(domain)
1249        domain.set_boundary({'exterior': Br})
1250
1251        #Setup only one forcing term, constant wind stress
1252        s = 100
1253        phi = 135
1254        domain.forcing_terms = []
1255        domain.forcing_terms.append( Wind_stress(s, phi) )
1256
1257        domain.compute_forcing_terms()
1258
1259
1260        const = eta_w*rho_a/rho_w
1261
1262        #Convert to radians
1263        phi = phi*pi/180
1264
1265        #Compute velocity vector (u, v)
1266        u = s*cos(phi)
1267        v = s*sin(phi)
1268
1269        #Compute wind stress
1270        S = const * sqrt(u**2 + v**2)
1271
1272        assert allclose(domain.quantities['stage'].explicit_update, 0)
1273        assert allclose(domain.quantities['xmomentum'].explicit_update, S*u)
1274        assert allclose(domain.quantities['ymomentum'].explicit_update, S*v)
1275
1276
1277    def test_variable_wind_stress(self):
1278        from anuga.config import rho_a, rho_w, eta_w
1279        from math import pi, cos, sin, sqrt
1280
1281        a = [0.0, 0.0]
1282        b = [0.0, 2.0]
1283        c = [2.0, 0.0]
1284        d = [0.0, 4.0]
1285        e = [2.0, 2.0]
1286        f = [4.0, 0.0]
1287
1288        points = [a, b, c, d, e, f]
1289        #bac, bce, ecf, dbe
1290        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
1291
1292        domain = Domain(points, vertices)
1293
1294        #Flat surface with 1m of water
1295        domain.set_quantity('elevation', 0)
1296        domain.set_quantity('stage', 1.0)
1297        domain.set_quantity('friction', 0)
1298
1299        Br = Reflective_boundary(domain)
1300        domain.set_boundary({'exterior': Br})
1301
1302
1303        domain.time = 5.54 #Take a random time (not zero)
1304
1305        #Setup only one forcing term, constant wind stress
1306        s = 100
1307        phi = 135
1308        domain.forcing_terms = []
1309        domain.forcing_terms.append( Wind_stress(s = speed, phi = angle) )
1310
1311        domain.compute_forcing_terms()
1312
1313        #Compute reference solution
1314        const = eta_w*rho_a/rho_w
1315
1316        N = len(domain) # number_of_triangles
1317
1318        xc = domain.get_centroid_coordinates()
1319        t = domain.time
1320
1321        x = xc[:,0]
1322        y = xc[:,1]
1323        s_vec = speed(t,x,y)
1324        phi_vec = angle(t,x,y)
1325
1326
1327        for k in range(N):
1328            #Convert to radians
1329            phi = phi_vec[k]*pi/180
1330            s = s_vec[k]
1331
1332            #Compute velocity vector (u, v)
1333            u = s*cos(phi)
1334            v = s*sin(phi)
1335
1336            #Compute wind stress
1337            S = const * sqrt(u**2 + v**2)
1338
1339            assert allclose(domain.quantities['stage'].explicit_update[k], 0)
1340            assert allclose(domain.quantities['xmomentum'].explicit_update[k], S*u)
1341            assert allclose(domain.quantities['ymomentum'].explicit_update[k], S*v)
1342
1343
1344
1345
1346    def test_windfield_from_file(self):
1347        from anuga.config import rho_a, rho_w, eta_w
1348        from math import pi, cos, sin, sqrt
1349        from anuga.config import time_format
1350        from anuga.abstract_2d_finite_volumes.util import file_function
1351        import time
1352
1353
1354        a = [0.0, 0.0]
1355        b = [0.0, 2.0]
1356        c = [2.0, 0.0]
1357        d = [0.0, 4.0]
1358        e = [2.0, 2.0]
1359        f = [4.0, 0.0]
1360
1361        points = [a, b, c, d, e, f]
1362        #bac, bce, ecf, dbe
1363        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
1364
1365        domain = Domain(points, vertices)
1366
1367        #Flat surface with 1m of water
1368        domain.set_quantity('elevation', 0)
1369        domain.set_quantity('stage', 1.0)
1370        domain.set_quantity('friction', 0)
1371
1372        Br = Reflective_boundary(domain)
1373        domain.set_boundary({'exterior': Br})
1374
1375
1376        domain.time = 7 #Take a time that is represented in file (not zero)
1377
1378        #Write wind stress file (ensure that domain.time is covered)
1379        #Take x=1 and y=0
1380        filename = 'test_windstress_from_file'
1381        start = time.mktime(time.strptime('2000', '%Y'))
1382        fid = open(filename + '.txt', 'w')
1383        dt = 1  #One second interval
1384        t = 0.0
1385        while t <= 10.0:
1386            t_string = time.strftime(time_format, time.gmtime(t+start))
1387
1388            fid.write('%s, %f %f\n' %(t_string,
1389                                      speed(t,[1],[0])[0],
1390                                      angle(t,[1],[0])[0]))
1391            t += dt
1392
1393        fid.close()
1394
1395
1396        #Convert ASCII file to NetCDF (Which is what we really like!)
1397        from data_manager import timefile2netcdf       
1398        timefile2netcdf(filename)
1399        os.remove(filename + '.txt')
1400
1401       
1402        #Setup wind stress
1403        F = file_function(filename + '.tms', quantities = ['Attribute0',
1404                                                           'Attribute1'])
1405        os.remove(filename + '.tms')
1406       
1407
1408        #print 'F(5)', F(5)
1409
1410        #print 'F(5,x,y)', F(5,x=zeros(3),y=zeros(3))
1411       
1412        #print dir(F)
1413        #print F.T
1414        #print F.precomputed_values
1415        #
1416        #F = file_function(filename + '.txt')       
1417        #
1418        #print dir(F)
1419        #print F.T       
1420        #print F.Q
1421       
1422        W = Wind_stress(F)
1423
1424        domain.forcing_terms = []
1425        domain.forcing_terms.append(W)
1426
1427        domain.compute_forcing_terms()
1428
1429        #Compute reference solution
1430        const = eta_w*rho_a/rho_w
1431
1432        N = len(domain) # number_of_triangles
1433
1434        t = domain.time
1435
1436        s = speed(t,[1],[0])[0]
1437        phi = angle(t,[1],[0])[0]
1438
1439        #Convert to radians
1440        phi = phi*pi/180
1441
1442
1443        #Compute velocity vector (u, v)
1444        u = s*cos(phi)
1445        v = s*sin(phi)
1446
1447        #Compute wind stress
1448        S = const * sqrt(u**2 + v**2)
1449
1450        for k in range(N):
1451            assert allclose(domain.quantities['stage'].explicit_update[k], 0)
1452            assert allclose(domain.quantities['xmomentum'].explicit_update[k], S*u)
1453            assert allclose(domain.quantities['ymomentum'].explicit_update[k], S*v)
1454
1455
1456    def test_windfield_from_file_seconds(self):
1457        from anuga.config import rho_a, rho_w, eta_w
1458        from math import pi, cos, sin, sqrt
1459        from anuga.config import time_format
1460        from anuga.abstract_2d_finite_volumes.util import file_function
1461        import time
1462
1463
1464        a = [0.0, 0.0]
1465        b = [0.0, 2.0]
1466        c = [2.0, 0.0]
1467        d = [0.0, 4.0]
1468        e = [2.0, 2.0]
1469        f = [4.0, 0.0]
1470
1471        points = [a, b, c, d, e, f]
1472        #bac, bce, ecf, dbe
1473        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
1474
1475        domain = Domain(points, vertices)
1476
1477        #Flat surface with 1m of water
1478        domain.set_quantity('elevation', 0)
1479        domain.set_quantity('stage', 1.0)
1480        domain.set_quantity('friction', 0)
1481
1482        Br = Reflective_boundary(domain)
1483        domain.set_boundary({'exterior': Br})
1484
1485
1486        domain.time = 7 #Take a time that is represented in file (not zero)
1487
1488        #Write wind stress file (ensure that domain.time is covered)
1489        #Take x=1 and y=0
1490        filename = 'test_windstress_from_file'
1491        start = time.mktime(time.strptime('2000', '%Y'))
1492        fid = open(filename + '.txt', 'w')
1493        dt = 0.5 #1  #One second interval
1494        t = 0.0
1495        while t <= 10.0:
1496            fid.write('%s, %f %f\n' %(str(t),
1497                                      speed(t,[1],[0])[0],
1498                                      angle(t,[1],[0])[0]))
1499            t += dt
1500
1501        fid.close()
1502
1503
1504        #Convert ASCII file to NetCDF (Which is what we really like!)
1505        from data_manager import timefile2netcdf       
1506        timefile2netcdf(filename, time_as_seconds=True)
1507        os.remove(filename + '.txt')
1508
1509       
1510        #Setup wind stress
1511        F = file_function(filename + '.tms', quantities = ['Attribute0',
1512                                                           'Attribute1'])
1513        os.remove(filename + '.tms')
1514       
1515
1516        #print 'F(5)', F(5)
1517
1518        #print 'F(5,x,y)', F(5,x=zeros(3),y=zeros(3))
1519       
1520        #print dir(F)
1521        #print F.T
1522        #print F.precomputed_values
1523        #
1524        #F = file_function(filename + '.txt')       
1525        #
1526        #print dir(F)
1527        #print F.T       
1528        #print F.Q
1529       
1530        W = Wind_stress(F)
1531
1532        domain.forcing_terms = []
1533        domain.forcing_terms.append(W)
1534
1535        domain.compute_forcing_terms()
1536
1537        #Compute reference solution
1538        const = eta_w*rho_a/rho_w
1539
1540        N = len(domain) # number_of_triangles
1541
1542        t = domain.time
1543
1544        s = speed(t,[1],[0])[0]
1545        phi = angle(t,[1],[0])[0]
1546
1547        #Convert to radians
1548        phi = phi*pi/180
1549
1550
1551        #Compute velocity vector (u, v)
1552        u = s*cos(phi)
1553        v = s*sin(phi)
1554
1555        #Compute wind stress
1556        S = const * sqrt(u**2 + v**2)
1557
1558        for k in range(N):
1559            assert allclose(domain.quantities['stage'].explicit_update[k], 0)
1560            assert allclose(domain.quantities['xmomentum'].explicit_update[k], S*u)
1561            assert allclose(domain.quantities['ymomentum'].explicit_update[k], S*v)
1562
1563
1564       
1565
1566    def test_wind_stress_error_condition(self):
1567        """Test that windstress reacts properly when forcing functions
1568        are wrong - e.g. returns a scalar
1569        """
1570
1571        from anuga.config import rho_a, rho_w, eta_w
1572        from math import pi, cos, sin, sqrt
1573
1574        a = [0.0, 0.0]
1575        b = [0.0, 2.0]
1576        c = [2.0, 0.0]
1577        d = [0.0, 4.0]
1578        e = [2.0, 2.0]
1579        f = [4.0, 0.0]
1580
1581        points = [a, b, c, d, e, f]
1582        #bac, bce, ecf, dbe
1583        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
1584
1585        domain = Domain(points, vertices)
1586
1587        #Flat surface with 1m of water
1588        domain.set_quantity('elevation', 0)
1589        domain.set_quantity('stage', 1.0)
1590        domain.set_quantity('friction', 0)
1591
1592        Br = Reflective_boundary(domain)
1593        domain.set_boundary({'exterior': Br})
1594
1595
1596        domain.time = 5.54 #Take a random time (not zero)
1597
1598        #Setup only one forcing term, bad func
1599        domain.forcing_terms = []
1600
1601        try:
1602            domain.forcing_terms.append(Wind_stress(s = scalar_func,
1603                                                    phi = angle))
1604        except AssertionError:
1605            pass
1606        else:
1607            msg = 'Should have raised exception'
1608            raise msg
1609
1610
1611        try:
1612            domain.forcing_terms.append(Wind_stress(s = speed,
1613                                                    phi = scalar_func))
1614        except AssertionError:
1615            pass
1616        else:
1617            msg = 'Should have raised exception'
1618            raise msg
1619
1620        try:
1621            domain.forcing_terms.append(Wind_stress(s = speed,
1622                                                    phi = 'xx'))
1623        except:
1624            pass
1625        else:
1626            msg = 'Should have raised exception'
1627            raise msg
1628
1629
1630    #####################################################
1631    def test_first_order_extrapolator_const_z(self):
1632
1633        a = [0.0, 0.0]
1634        b = [0.0, 2.0]
1635        c = [2.0, 0.0]
1636        d = [0.0, 4.0]
1637        e = [2.0, 2.0]
1638        f = [4.0, 0.0]
1639
1640        points = [a, b, c, d, e, f]
1641        #bac, bce, ecf, dbe
1642        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
1643
1644        domain = Domain(points, vertices)
1645        val0 = 2.+2.0/3
1646        val1 = 4.+4.0/3
1647        val2 = 8.+2.0/3
1648        val3 = 2.+8.0/3
1649
1650        zl=zr=-3.75 #Assume constant bed (must be less than stage)
1651        domain.set_quantity('elevation', zl*ones( (4,3) ))
1652        domain.set_quantity('stage', [[val0, val0-1, val0-2],
1653                                      [val1, val1+1, val1],
1654                                      [val2, val2-2, val2],
1655                                      [val3-0.5, val3, val3]])
1656
1657
1658
1659        domain._order_ = 1
1660        domain.distribute_to_vertices_and_edges()
1661
1662        #Check that centroid values were distributed to vertices
1663        C = domain.quantities['stage'].centroid_values
1664        for i in range(3):
1665            assert allclose( domain.quantities['stage'].vertex_values[:,i], C)
1666
1667
1668    def test_first_order_limiter_variable_z(self):
1669        #Check that first order limiter follows bed_slope
1670        from Numeric import alltrue, greater_equal
1671        from anuga.config import epsilon
1672
1673        a = [0.0, 0.0]
1674        b = [0.0, 2.0]
1675        c = [2.0,0.0]
1676        d = [0.0, 4.0]
1677        e = [2.0, 2.0]
1678        f = [4.0,0.0]
1679
1680        points = [a, b, c, d, e, f]
1681        #bac, bce, ecf, dbe
1682        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
1683
1684        domain = Domain(points, vertices)
1685        val0 = 2.+2.0/3
1686        val1 = 4.+4.0/3
1687        val2 = 8.+2.0/3
1688        val3 = 2.+8.0/3
1689
1690        domain.set_quantity('elevation', [[0,0,0], [6,0,0],
1691                                          [6,6,6], [6,6,6]])
1692        domain.set_quantity('stage', [[val0, val0, val0],
1693                                      [val1, val1, val1],
1694                                      [val2, val2, val2],
1695                                      [val3, val3, val3]])
1696
1697        E = domain.quantities['elevation'].vertex_values
1698        L = domain.quantities['stage'].vertex_values
1699
1700
1701        #Check that some stages are not above elevation (within eps)
1702        #- so that the limiter has something to work with
1703        assert not alltrue(alltrue(greater_equal(L,E-epsilon)))
1704
1705        domain._order_ = 1
1706        domain.distribute_to_vertices_and_edges()
1707
1708        #Check that all stages are above elevation (within eps)
1709        assert alltrue(alltrue(greater_equal(L,E-epsilon)))
1710
1711
1712    #####################################################
1713    def test_distribute_basic(self):
1714        #Using test data generated by abstract_2d_finite_volumes-2
1715        #Assuming no friction and flat bed (0.0)
1716
1717        a = [0.0, 0.0]
1718        b = [0.0, 2.0]
1719        c = [2.0, 0.0]
1720        d = [0.0, 4.0]
1721        e = [2.0, 2.0]
1722        f = [4.0, 0.0]
1723
1724        points = [a, b, c, d, e, f]
1725        #bac, bce, ecf, dbe
1726        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
1727
1728        domain = Domain(points, vertices)
1729
1730        val0 = 2.
1731        val1 = 4.
1732        val2 = 8.
1733        val3 = 2.
1734
1735        domain.set_quantity('stage', [val0, val1, val2, val3],
1736                            location='centroids')
1737        L = domain.quantities['stage'].vertex_values
1738
1739        #First order
1740        domain._order_ = 1
1741        domain.distribute_to_vertices_and_edges()
1742        assert allclose(L[1], val1)
1743
1744        #Second order
1745        domain._order_ = 2
1746        domain.beta_w      = 0.9
1747        domain.beta_w_dry  = 0.9
1748        domain.beta_uh     = 0.9
1749        domain.beta_uh_dry = 0.9
1750        domain.beta_vh     = 0.9
1751        domain.beta_vh_dry = 0.9
1752        domain.distribute_to_vertices_and_edges()
1753        assert allclose(L[1], [2.2, 4.9, 4.9])
1754
1755
1756
1757    def test_distribute_away_from_bed(self):
1758        #Using test data generated by abstract_2d_finite_volumes-2
1759        #Assuming no friction and flat bed (0.0)
1760
1761        a = [0.0, 0.0]
1762        b = [0.0, 2.0]
1763        c = [2.0, 0.0]
1764        d = [0.0, 4.0]
1765        e = [2.0, 2.0]
1766        f = [4.0, 0.0]
1767
1768        points = [a, b, c, d, e, f]
1769        #bac, bce, ecf, dbe
1770        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
1771
1772        domain = Domain(points, vertices)
1773        L = domain.quantities['stage'].vertex_values
1774
1775        def stage(x,y):
1776            return x**2
1777
1778        domain.set_quantity('stage', stage, location='centroids')
1779
1780        a, b = domain.quantities['stage'].compute_gradients()
1781        assert allclose(a[1], 3.33333334)
1782        assert allclose(b[1], 0.0)
1783
1784        domain._order_ = 1
1785        domain.distribute_to_vertices_and_edges()
1786        assert allclose(L[1], 1.77777778)
1787
1788        domain._order_ = 2
1789        domain.beta_w      = 0.9
1790        domain.beta_w_dry  = 0.9
1791        domain.beta_uh     = 0.9
1792        domain.beta_uh_dry = 0.9
1793        domain.beta_vh     = 0.9
1794        domain.beta_vh_dry = 0.9
1795        domain.distribute_to_vertices_and_edges()
1796        assert allclose(L[1], [0.57777777, 2.37777778, 2.37777778])
1797
1798
1799
1800    def test_distribute_away_from_bed1(self):
1801        #Using test data generated by abstract_2d_finite_volumes-2
1802        #Assuming no friction and flat bed (0.0)
1803
1804        a = [0.0, 0.0]
1805        b = [0.0, 2.0]
1806        c = [2.0, 0.0]
1807        d = [0.0, 4.0]
1808        e = [2.0, 2.0]
1809        f = [4.0, 0.0]
1810
1811        points = [a, b, c, d, e, f]
1812        #bac, bce, ecf, dbe
1813        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
1814
1815        domain = Domain(points, vertices)
1816        L = domain.quantities['stage'].vertex_values
1817
1818        def stage(x,y):
1819            return x**4+y**2
1820
1821        domain.set_quantity('stage', stage, location='centroids')
1822        #print domain.quantities['stage'].centroid_values
1823
1824        a, b = domain.quantities['stage'].compute_gradients()
1825        assert allclose(a[1], 25.18518519)
1826        assert allclose(b[1], 3.33333333)
1827
1828        domain._order_ = 1
1829        domain.distribute_to_vertices_and_edges()
1830        assert allclose(L[1], 4.9382716)
1831
1832        domain._order_ = 2
1833        domain.beta_w      = 0.9
1834        domain.beta_w_dry  = 0.9
1835        domain.beta_uh     = 0.9
1836        domain.beta_uh_dry = 0.9
1837        domain.beta_vh     = 0.9
1838        domain.beta_vh_dry = 0.9
1839        domain.distribute_to_vertices_and_edges()
1840        assert allclose(L[1], [1.07160494, 6.46058131, 7.28262855])
1841
1842
1843
1844    def test_distribute_near_bed(self):
1845        #Using test data generated by abstract_2d_finite_volumes-2
1846        #Assuming no friction and flat bed (0.0)
1847
1848        a = [0.0, 0.0]
1849        b = [0.0, 2.0]
1850        c = [2.0, 0.0]
1851        d = [0.0, 4.0]
1852        e = [2.0, 2.0]
1853        f = [4.0, 0.0]
1854
1855        points = [a, b, c, d, e, f]
1856        #bac, bce, ecf, dbe
1857        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
1858
1859        domain = Domain(points, vertices)
1860
1861
1862        #Set up for a gradient of (3,0) at mid triangle
1863        def slope(x, y):
1864            return 10*x
1865
1866        h = 0.1
1867        def stage(x,y):
1868            return slope(x,y)+h
1869
1870        domain.set_quantity('elevation', slope)
1871        domain.set_quantity('stage', stage, location='centroids')
1872
1873        #print domain.quantities['elevation'].centroid_values
1874        #print domain.quantities['stage'].centroid_values
1875
1876        E = domain.quantities['elevation'].vertex_values
1877        L = domain.quantities['stage'].vertex_values
1878
1879        #print E
1880        domain._order_ = 1
1881        domain.distribute_to_vertices_and_edges()
1882        ##assert allclose(L[1], [0.19999999, 20.05, 20.05])
1883        assert allclose(L[1], [0.1, 20.1, 20.1])
1884
1885        domain._order_ = 2
1886        domain.distribute_to_vertices_and_edges()
1887        assert allclose(L[1], [0.1, 20.1, 20.1])
1888
1889    def test_distribute_near_bed1(self):
1890        #Using test data generated by abstract_2d_finite_volumes-2
1891        #Assuming no friction and flat bed (0.0)
1892
1893        a = [0.0, 0.0]
1894        b = [0.0, 2.0]
1895        c = [2.0, 0.0]
1896        d = [0.0, 4.0]
1897        e = [2.0, 2.0]
1898        f = [4.0, 0.0]
1899
1900        points = [a, b, c, d, e, f]
1901        #bac, bce, ecf, dbe
1902        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
1903
1904        domain = Domain(points, vertices)
1905
1906
1907        #Set up for a gradient of (3,0) at mid triangle
1908        def slope(x, y):
1909            return x**4+y**2
1910
1911        h = 0.1
1912        def stage(x,y):
1913            return slope(x,y)+h
1914
1915        domain.set_quantity('elevation', slope)
1916        domain.set_quantity('stage', stage)
1917
1918        #print domain.quantities['elevation'].centroid_values
1919        #print domain.quantities['stage'].centroid_values
1920
1921        E = domain.quantities['elevation'].vertex_values
1922        L = domain.quantities['stage'].vertex_values
1923
1924        #print E
1925        domain._order_ = 1
1926        domain.distribute_to_vertices_and_edges()
1927        ##assert allclose(L[1], [4.19999999, 16.07142857, 20.02857143])
1928        assert allclose(L[1], [4.1, 16.1, 20.1])
1929
1930        domain._order_ = 2
1931        domain.distribute_to_vertices_and_edges()
1932        assert allclose(L[1], [4.1, 16.1, 20.1])
1933
1934
1935
1936    def test_second_order_distribute_real_data(self):
1937        #Using test data generated by abstract_2d_finite_volumes-2
1938        #Assuming no friction and flat bed (0.0)
1939
1940        a = [0.0, 0.0]
1941        b = [0.0, 1.0/5]
1942        c = [0.0, 2.0/5]
1943        d = [1.0/5, 0.0]
1944        e = [1.0/5, 1.0/5]
1945        f = [1.0/5, 2.0/5]
1946        g = [2.0/5, 2.0/5]
1947
1948        points = [a, b, c, d, e, f, g]
1949        #bae, efb, cbf, feg
1950        vertices = [ [1,0,4], [4,5,1], [2,1,5], [5,4,6]]
1951
1952        domain = Domain(points, vertices)
1953
1954        def slope(x, y):
1955            return -x/3
1956
1957        domain.set_quantity('elevation', slope)
1958        domain.set_quantity('stage',
1959                            [0.01298164, 0.00365611,
1960                             0.01440365, -0.0381856437096],
1961                            location='centroids')
1962        domain.set_quantity('xmomentum',
1963                            [0.00670439, 0.01263789,
1964                             0.00647805, 0.0178180740668],
1965                            location='centroids')
1966        domain.set_quantity('ymomentum',
1967                            [-7.23510980e-004, -6.30413883e-005,
1968                             6.30413883e-005, 0.000200907255866],
1969                            location='centroids')
1970
1971        E = domain.quantities['elevation'].vertex_values
1972        L = domain.quantities['stage'].vertex_values
1973        X = domain.quantities['xmomentum'].vertex_values
1974        Y = domain.quantities['ymomentum'].vertex_values
1975
1976        #print E
1977        domain._order_ = 2
1978        domain.beta_w      = 0.9
1979        domain.beta_w_dry  = 0.9
1980        domain.beta_uh     = 0.9
1981        domain.beta_uh_dry = 0.9
1982        domain.beta_vh     = 0.9
1983        domain.beta_vh_dry = 0.9
1984        domain.beta_h = 0.0 #Use first order in h-limiter
1985        domain.distribute_to_vertices_and_edges()
1986
1987        #print L[1,:]
1988        #print X[1,:]
1989        #print Y[1,:]
1990
1991        assert allclose(L[1,:], [-0.00825735775384,
1992                                 -0.00801881482869,
1993                                 0.0272445025825])
1994        assert allclose(X[1,:], [0.0143507718962,
1995                                 0.0142502147066,
1996                                 0.00931268339717])
1997        assert allclose(Y[1,:], [-0.000117062180693,
1998                                 7.94434448109e-005,
1999                                 -0.000151505429018])
2000
2001
2002
2003    def test_balance_deep_and_shallow(self):
2004        """Test that balanced limiters preserve conserved quantites.
2005        """
2006        import copy
2007
2008        a = [0.0, 0.0]
2009        b = [0.0, 2.0]
2010        c = [2.0, 0.0]
2011        d = [0.0, 4.0]
2012        e = [2.0, 2.0]
2013        f = [4.0, 0.0]
2014
2015        points = [a, b, c, d, e, f]
2016
2017        #bac, bce, ecf, dbe
2018        elements = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4] ]
2019
2020        mesh = Domain(points, elements)
2021        mesh.check_integrity()
2022
2023        #Create a deliberate overshoot
2024        mesh.set_quantity('stage', [[3,0,3], [2,2,6], [5,3,8], [8,3,5]])
2025        mesh.set_quantity('elevation', 0) #Flat bed
2026        stage = mesh.quantities['stage']
2027
2028        ref_centroid_values = copy.copy(stage.centroid_values[:]) #Copy
2029
2030        #Limit
2031        mesh.distribute_to_vertices_and_edges()
2032
2033        #Assert that quantities are conserved
2034        from Numeric import sum
2035        for k in range(len(mesh)):
2036            assert allclose (ref_centroid_values[k],
2037                             sum(stage.vertex_values[k,:])/3)
2038
2039
2040        #Now try with a non-flat bed - closely hugging initial stage in places
2041        #This will create alphas in the range [0, 0.478260, 1]
2042        mesh.set_quantity('stage', [[3,0,3], [2,2,6], [5,3,8], [8,3,5]])
2043        mesh.set_quantity('elevation', [[0,0,0],
2044                                        [1.8,1.9,5.9],
2045                                        [4.6,0,0],
2046                                        [0,2,4]])
2047        stage = mesh.quantities['stage']
2048
2049        ref_centroid_values = copy.copy(stage.centroid_values[:]) #Copy
2050        ref_vertex_values = copy.copy(stage.vertex_values[:]) #Copy
2051
2052        #Limit
2053        mesh.distribute_to_vertices_and_edges()
2054
2055
2056        #Assert that all vertex quantities have changed
2057        for k in range(len(mesh)):
2058            #print ref_vertex_values[k,:], stage.vertex_values[k,:]
2059            assert not allclose (ref_vertex_values[k,:], stage.vertex_values[k,:])
2060        #and assert that quantities are still conserved
2061        from Numeric import sum
2062        for k in range(len(mesh)):
2063            assert allclose (ref_centroid_values[k],
2064                             sum(stage.vertex_values[k,:])/3)
2065
2066
2067        #Also check that Python and C version produce the same
2068        assert allclose (stage.vertex_values,
2069                         [[2,2,2],
2070                          [1.93333333, 2.03333333, 6.03333333],
2071                          [6.93333333, 4.53333333, 4.53333333],
2072                          [5.33333333, 5.33333333, 5.33333333]])
2073
2074
2075
2076
2077    def test_conservation_1(self):
2078        """Test that stage is conserved globally
2079
2080        This one uses a flat bed, reflective bdries and a suitable
2081        initial condition
2082        """
2083        from mesh_factory import rectangular
2084        from Numeric import array
2085
2086        #Create basic mesh
2087        points, vertices, boundary = rectangular(6, 6)
2088
2089        #Create shallow water domain
2090        domain = Domain(points, vertices, boundary)
2091        domain.smooth = False
2092        domain.default_order=2
2093
2094        #IC
2095        def x_slope(x, y):
2096            return x/3
2097
2098        domain.set_quantity('elevation', 0)
2099        domain.set_quantity('friction', 0)
2100        domain.set_quantity('stage', x_slope)
2101
2102        # Boundary conditions (reflective everywhere)
2103        Br = Reflective_boundary(domain)
2104        domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
2105
2106        domain.check_integrity()
2107
2108        initial_volume = domain.quantities['stage'].get_integral()
2109        initial_xmom = domain.quantities['xmomentum'].get_integral()
2110
2111        #print initial_xmom
2112
2113        #Evolution
2114        for t in domain.evolve(yieldstep = 0.05, finaltime = 5.0):
2115            volume =  domain.quantities['stage'].get_integral()
2116            assert allclose (volume, initial_volume)
2117
2118            #I don't believe that the total momentum should be the same
2119            #It starts with zero and ends with zero though
2120            #xmom = domain.quantities['xmomentum'].get_integral()
2121            #print xmom
2122            #assert allclose (xmom, initial_xmom)
2123
2124        os.remove(domain.get_name() + '.sww')
2125
2126
2127    def test_conservation_2(self):
2128        """Test that stage is conserved globally
2129
2130        This one uses a slopy bed, reflective bdries and a suitable
2131        initial condition
2132        """
2133        from mesh_factory import rectangular
2134        from Numeric import array
2135
2136        #Create basic mesh
2137        points, vertices, boundary = rectangular(6, 6)
2138
2139        #Create shallow water domain
2140        domain = Domain(points, vertices, boundary)
2141        domain.smooth = False
2142        domain.default_order=2
2143
2144        #IC
2145        def x_slope(x, y):
2146            return x/3
2147
2148        domain.set_quantity('elevation', x_slope)
2149        domain.set_quantity('friction', 0)
2150        domain.set_quantity('stage', 0.4) #Steady
2151
2152        # Boundary conditions (reflective everywhere)
2153        Br = Reflective_boundary(domain)
2154        domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
2155
2156        domain.check_integrity()
2157
2158        initial_volume = domain.quantities['stage'].get_integral()
2159        initial_xmom = domain.quantities['xmomentum'].get_integral()
2160
2161        #print initial_xmom
2162
2163        #Evolution
2164        for t in domain.evolve(yieldstep = 0.05, finaltime = 5.0):
2165            volume =  domain.quantities['stage'].get_integral()
2166            assert allclose (volume, initial_volume)
2167
2168            #FIXME: What would we expect from momentum
2169            #xmom = domain.quantities['xmomentum'].get_integral()
2170            #print xmom
2171            #assert allclose (xmom, initial_xmom)
2172
2173        os.remove(domain.get_name() + '.sww')
2174
2175    def test_conservation_3(self):
2176        """Test that stage is conserved globally
2177
2178        This one uses a larger grid, convoluted bed, reflective bdries and a suitable
2179        initial condition
2180        """
2181        from mesh_factory import rectangular
2182        from Numeric import array
2183
2184        #Create basic mesh
2185        points, vertices, boundary = rectangular(2, 1)
2186
2187        #Create shallow water domain
2188        domain = Domain(points, vertices, boundary)
2189        domain.smooth = False
2190        domain.default_order = 2
2191        domain.beta_h = 0.2
2192        domain.set_quantities_to_be_stored(['stage'])
2193
2194        #IC
2195        def x_slope(x, y):
2196            z = 0*x
2197            for i in range(len(x)):
2198                if x[i] < 0.3:
2199                    z[i] = x[i]/3
2200                if 0.3 <= x[i] < 0.5:
2201                    z[i] = -0.5
2202                if 0.5 <= x[i] < 0.7:
2203                    z[i] = 0.39
2204                if 0.7 <= x[i]:
2205                    z[i] = x[i]/3
2206            return z
2207
2208
2209
2210        domain.set_quantity('elevation', x_slope)
2211        domain.set_quantity('friction', 0)
2212        domain.set_quantity('stage', 0.4) #Steady
2213
2214        # Boundary conditions (reflective everywhere)
2215        Br = Reflective_boundary(domain)
2216        domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
2217
2218        domain.check_integrity()
2219
2220        initial_volume = domain.quantities['stage'].get_integral()
2221        initial_xmom = domain.quantities['xmomentum'].get_integral()
2222
2223        import copy
2224        ref_centroid_values =\
2225                 copy.copy(domain.quantities['stage'].centroid_values)
2226
2227        #print 'ORG', domain.quantities['stage'].centroid_values
2228        domain.distribute_to_vertices_and_edges()
2229
2230
2231        #print domain.quantities['stage'].centroid_values
2232        assert allclose(domain.quantities['stage'].centroid_values,
2233                        ref_centroid_values)
2234
2235
2236        #Check that initial limiter doesn't violate cons quan
2237        assert allclose (domain.quantities['stage'].get_integral(),
2238                         initial_volume)
2239
2240        #Evolution
2241        for t in domain.evolve(yieldstep = 0.05, finaltime = 10):
2242            volume =  domain.quantities['stage'].get_integral()
2243            #print t, volume, initial_volume
2244            assert allclose (volume, initial_volume)
2245
2246        os.remove(domain.get_name() + '.sww')
2247
2248    def test_conservation_4(self):
2249        """Test that stage is conserved globally
2250
2251        This one uses a larger grid, convoluted bed, reflective bdries and a suitable
2252        initial condition
2253        """
2254        from mesh_factory import rectangular
2255        from Numeric import array
2256
2257        #Create basic mesh
2258        points, vertices, boundary = rectangular(6, 6)
2259
2260        #Create shallow water domain
2261        domain = Domain(points, vertices, boundary)
2262        domain.smooth = False
2263        domain.default_order=2
2264        domain.beta_h = 0.0
2265
2266        #IC
2267        def x_slope(x, y):
2268            z = 0*x
2269            for i in range(len(x)):
2270                if x[i] < 0.3:
2271                    z[i] = x[i]/3
2272                if 0.3 <= x[i] < 0.5:
2273                    z[i] = -0.5
2274                if 0.5 <= x[i] < 0.7:
2275                    #z[i] = 0.3 #OK with beta == 0.2
2276                    z[i] = 0.34 #OK with beta == 0.0
2277                    #z[i] = 0.35#Fails after 80 timesteps with an error
2278                                #of the order 1.0e-5
2279                if 0.7 <= x[i]:
2280                    z[i] = x[i]/3
2281            return z
2282
2283
2284
2285        domain.set_quantity('elevation', x_slope)
2286        domain.set_quantity('friction', 0)
2287        domain.set_quantity('stage', 0.4) #Steady
2288
2289        # Boundary conditions (reflective everywhere)
2290        Br = Reflective_boundary(domain)
2291        domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
2292
2293        domain.check_integrity()
2294
2295        initial_volume = domain.quantities['stage'].get_integral()
2296        initial_xmom = domain.quantities['xmomentum'].get_integral()
2297
2298        import copy
2299        ref_centroid_values =\
2300                 copy.copy(domain.quantities['stage'].centroid_values)
2301
2302        #Test limiter by itself
2303        domain.distribute_to_vertices_and_edges()
2304
2305        #Check that initial limiter doesn't violate cons quan
2306        assert allclose (domain.quantities['stage'].get_integral(),
2307                         initial_volume)
2308        #NOTE: This would fail if any initial stage was less than the
2309        #corresponding bed elevation - but that is reasonable.
2310
2311
2312        #Evolution
2313        for t in domain.evolve(yieldstep = 0.05, finaltime = 10.0):
2314            volume =  domain.quantities['stage'].get_integral()
2315
2316            #print t, volume, initial_volume
2317
2318            assert allclose (volume, initial_volume)
2319
2320
2321        os.remove(domain.get_name() + '.sww')
2322
2323
2324    def test_conservation_5(self):
2325        """Test that momentum is conserved globally in
2326        steady state scenario
2327
2328        This one uses a slopy bed, dirichlet and reflective bdries
2329        """
2330        from mesh_factory import rectangular
2331        from Numeric import array
2332
2333        #Create basic mesh
2334        points, vertices, boundary = rectangular(6, 6)
2335
2336        #Create shallow water domain
2337        domain = Domain(points, vertices, boundary)
2338        domain.smooth = False
2339        domain.default_order=2
2340
2341        #IC
2342        def x_slope(x, y):
2343            return x/3
2344
2345        domain.set_quantity('elevation', x_slope)
2346        domain.set_quantity('friction', 0)
2347        domain.set_quantity('stage', 0.4) #Steady
2348
2349        # Boundary conditions (reflective everywhere)
2350        Br = Reflective_boundary(domain)
2351        Bleft = Dirichlet_boundary([0.5,0,0])
2352        Bright = Dirichlet_boundary([0.1,0,0])
2353        domain.set_boundary({'left': Bleft, 'right': Bright,
2354                             'top': Br, 'bottom': Br})
2355
2356        domain.check_integrity()
2357
2358        initial_volume = domain.quantities['stage'].get_integral()
2359        initial_xmom = domain.quantities['xmomentum'].get_integral()
2360
2361
2362        #Evolution
2363        for t in domain.evolve(yieldstep = 0.05, finaltime = 15.0):
2364            stage =  domain.quantities['stage'].get_integral()
2365            xmom = domain.quantities['xmomentum'].get_integral()
2366            ymom = domain.quantities['ymomentum'].get_integral()
2367
2368            if allclose(t, 6):  #Steady state reached
2369                steady_xmom = domain.quantities['xmomentum'].get_integral()
2370                steady_ymom = domain.quantities['ymomentum'].get_integral()
2371                steady_stage = domain.quantities['stage'].get_integral()
2372
2373            if t > 6:
2374                #print '%.2f %14.8f %14.8f' %(t, ymom, steady_ymom)
2375                assert allclose(xmom, steady_xmom)
2376                assert allclose(ymom, steady_ymom)
2377                assert allclose(stage, steady_stage)
2378
2379
2380        os.remove(domain.get_name() + '.sww')
2381
2382
2383
2384
2385
2386    def test_conservation_real(self):
2387        """Test that momentum is conserved globally
2388       
2389        Stephen finally made a test that revealed the problem.
2390        This test failed with code prior to 25 July 2005
2391        """
2392
2393        yieldstep = 0.01
2394        finaltime = 0.05
2395        min_depth = 1.0e-2
2396
2397
2398        import sys
2399        from os import sep; sys.path.append('..'+sep+'abstract_2d_finite_volumes')
2400        from mesh_factory import rectangular
2401
2402
2403        #Create shallow water domain
2404        points, vertices, boundary = rectangular(10, 10, len1=500, len2=500)
2405        domain = Domain(points, vertices, boundary)
2406        domain.smooth = False
2407        domain.default_order = 1
2408        domain.minimum_allowed_height = min_depth
2409
2410        # Set initial condition
2411        class Set_IC:
2412            """Set an initial condition with a constant value, for x0<x<x1
2413            """
2414
2415            def __init__(self, x0=0.25, x1=0.5, h=1.0):
2416                self.x0 = x0
2417                self.x1 = x1
2418                self.= h
2419
2420            def __call__(self, x, y):
2421                return self.h*((x>self.x0)&(x<self.x1))
2422
2423
2424        domain.set_quantity('stage', Set_IC(200.0,300.0,5.0))
2425
2426
2427        #Boundaries
2428        R = Reflective_boundary(domain)
2429        domain.set_boundary( {'left': R, 'right': R, 'top':R, 'bottom': R})
2430
2431        ref = domain.quantities['stage'].get_integral()
2432
2433        # Evolution
2434        for t in domain.evolve(yieldstep = yieldstep, finaltime = finaltime):
2435            pass
2436            #print 'Integral stage = ',\
2437            #      domain.quantities['stage'].get_integral(),\
2438            #      ' Time = ',domain.time
2439
2440
2441        now = domain.quantities['stage'].get_integral()
2442
2443        msg = 'Stage not conserved: was %f, now %f' %(ref, now)
2444        assert allclose(ref, now), msg
2445
2446        os.remove(domain.get_name() + '.sww')
2447
2448    def test_second_order_flat_bed_onestep(self):
2449
2450        from mesh_factory import rectangular
2451        from Numeric import array
2452
2453        #Create basic mesh
2454        points, vertices, boundary = rectangular(6, 6)
2455
2456        #Create shallow water domain
2457        domain = Domain(points, vertices, boundary)
2458        domain.smooth = False
2459        domain.default_order=2
2460        domain.beta_w      = 0.9
2461        domain.beta_w_dry  = 0.9
2462        domain.beta_uh     = 0.9
2463        domain.beta_uh_dry = 0.9
2464        domain.beta_vh     = 0.9
2465        domain.beta_vh_dry = 0.9
2466        domain.H0 = 0 # Backwards compatibility (6/2/7)       
2467       
2468        # Boundary conditions
2469        Br = Reflective_boundary(domain)
2470        Bd = Dirichlet_boundary([0.1, 0., 0.])
2471        domain.set_boundary({'left': Bd, 'right': Br, 'top': Br, 'bottom': Br})
2472
2473        domain.check_integrity()
2474
2475        #Evolution
2476        for t in domain.evolve(yieldstep = 0.05, finaltime = 0.05):
2477            pass# domain.write_time()
2478
2479        #Data from earlier version of abstract_2d_finite_volumes
2480        assert allclose(domain.min_timestep, 0.0396825396825)
2481        assert allclose(domain.max_timestep, 0.0396825396825)
2482
2483        assert allclose(domain.quantities['stage'].centroid_values[:12],
2484                        [0.00171396, 0.02656103, 0.00241523, 0.02656103,
2485                        0.00241523, 0.02656103, 0.00241523, 0.02656103,
2486                        0.00241523, 0.02656103, 0.00241523, 0.0272623])
2487
2488        domain.distribute_to_vertices_and_edges()
2489        assert allclose(domain.quantities['stage'].vertex_values[:12,0],
2490                        [0.0001714, 0.02656103, 0.00024152,
2491                        0.02656103, 0.00024152, 0.02656103,
2492                        0.00024152, 0.02656103, 0.00024152,
2493                        0.02656103, 0.00024152, 0.0272623])
2494
2495        assert allclose(domain.quantities['stage'].vertex_values[:12,1],
2496                        [0.00315012, 0.02656103, 0.00024152, 0.02656103,
2497                         0.00024152, 0.02656103, 0.00024152, 0.02656103,
2498                         0.00024152, 0.02656103, 0.00040506, 0.0272623])
2499
2500        assert allclose(domain.quantities['stage'].vertex_values[:12,2],
2501                        [0.00182037, 0.02656103, 0.00676264,
2502                         0.02656103, 0.00676264, 0.02656103,
2503                         0.00676264, 0.02656103, 0.00676264,
2504                         0.02656103, 0.0065991, 0.0272623])
2505
2506        assert allclose(domain.quantities['xmomentum'].centroid_values[:12],
2507                        [0.00113961, 0.01302432, 0.00148672,
2508                         0.01302432, 0.00148672, 0.01302432,
2509                         0.00148672, 0.01302432, 0.00148672 ,
2510                         0.01302432, 0.00148672, 0.01337143])
2511
2512        assert allclose(domain.quantities['ymomentum'].centroid_values[:12],
2513                        [-2.91240050e-004 , 1.22721531e-004,
2514                         -1.22721531e-004, 1.22721531e-004 ,
2515                         -1.22721531e-004, 1.22721531e-004,
2516                         -1.22721531e-004 , 1.22721531e-004,
2517                         -1.22721531e-004, 1.22721531e-004,
2518                         -1.22721531e-004, -4.57969873e-005])
2519
2520        os.remove(domain.get_name() + '.sww')
2521
2522
2523    def test_second_order_flat_bed_moresteps(self):
2524
2525        from mesh_factory import rectangular
2526        from Numeric import array
2527
2528        #Create basic mesh
2529        points, vertices, boundary = rectangular(6, 6)
2530
2531        #Create shallow water domain
2532        domain = Domain(points, vertices, boundary)
2533        domain.smooth = False
2534        domain.default_order=2
2535
2536        # Boundary conditions
2537        Br = Reflective_boundary(domain)
2538        Bd = Dirichlet_boundary([0.1, 0., 0.])
2539        domain.set_boundary({'left': Bd, 'right': Br, 'top': Br, 'bottom': Br})
2540
2541        domain.check_integrity()
2542
2543        #Evolution
2544        for t in domain.evolve(yieldstep = 0.05, finaltime = 0.1):
2545            pass
2546
2547        #Data from earlier version of abstract_2d_finite_volumes
2548        #assert allclose(domain.min_timestep, 0.0396825396825)
2549        #assert allclose(domain.max_timestep, 0.0396825396825)
2550        #print domain.quantities['stage'].centroid_values
2551
2552        os.remove(domain.get_name() + '.sww')       
2553
2554
2555    def test_flatbed_first_order(self):
2556        from mesh_factory import rectangular
2557        from Numeric import array
2558
2559        #Create basic mesh
2560        N = 8
2561        points, vertices, boundary = rectangular(N, N)
2562
2563        #Create shallow water domain
2564        domain = Domain(points, vertices, boundary)
2565        domain.smooth = False
2566        domain.default_order=1
2567        domain.H0 = 0 # Backwards compatibility (6/2/7)       
2568
2569        # Boundary conditions
2570        Br = Reflective_boundary(domain)
2571        Bd = Dirichlet_boundary([0.2,0.,0.])
2572
2573        domain.set_boundary({'left': Bd, 'right': Br, 'top': Br, 'bottom': Br})
2574        domain.check_integrity()
2575
2576
2577        #Evolution
2578        for t in domain.evolve(yieldstep = 0.02, finaltime = 0.5):
2579            pass
2580            #domain.write_time()
2581
2582        #FIXME: These numbers were from version before 25/10
2583        #assert allclose(domain.min_timestep, 0.0140413643926)
2584        #assert allclose(domain.max_timestep, 0.0140947355753)
2585
2586        for i in range(3):
2587            #assert allclose(domain.quantities['stage'].edge_values[:4,i],
2588            #                [0.10730244,0.12337617,0.11200126,0.12605666])
2589
2590            assert allclose(domain.quantities['xmomentum'].edge_values[:4,i],
2591                            [0.07610894,0.06901572,0.07284461,0.06819712])
2592
2593            #assert allclose(domain.quantities['ymomentum'].edge_values[:4,i],
2594            #                [-0.0060238, -0.00157404, -0.00309633, -0.0001637])
2595
2596
2597        os.remove(domain.get_name() + '.sww')           
2598
2599    def test_flatbed_second_order(self):
2600        from mesh_factory import rectangular
2601        from Numeric import array
2602
2603        #Create basic mesh
2604        N = 8
2605        points, vertices, boundary = rectangular(N, N)
2606
2607        #Create shallow water domain
2608        domain = Domain(points, vertices, boundary)
2609        domain.smooth = False
2610        domain.default_order=2
2611        domain.beta_w      = 0.9
2612        domain.beta_w_dry  = 0.9
2613        domain.beta_uh     = 0.9
2614        domain.beta_uh_dry = 0.9
2615        domain.beta_vh     = 0.9
2616        domain.beta_vh_dry = 0.9       
2617        #domain.minimum_allowed_height = 0.0 #Makes it like the 'oldstyle' balance
2618        domain.H0 = 0 # Backwards compatibility (6/2/7)       
2619
2620        # Boundary conditions
2621        Br = Reflective_boundary(domain)
2622        Bd = Dirichlet_boundary([0.2,0.,0.])
2623
2624        domain.set_boundary({'left': Bd, 'right': Br, 'top': Br, 'bottom': Br})
2625        domain.check_integrity()
2626
2627        #Evolution
2628        for t in domain.evolve(yieldstep = 0.01, finaltime = 0.03):
2629            pass
2630
2631
2632        assert allclose(domain.min_timestep, 0.0210448446782)
2633        assert allclose(domain.max_timestep, 0.0210448446782)
2634
2635        #print domain.quantities['stage'].vertex_values[:4,0]
2636        #print domain.quantities['xmomentum'].vertex_values[:4,0]
2637        #print domain.quantities['ymomentum'].vertex_values[:4,0]
2638
2639        #FIXME: These numbers were from version before 25/10
2640        #assert allclose(domain.quantities['stage'].vertex_values[:4,0],
2641        #                [0.00101913,0.05352143,0.00104852,0.05354394])
2642
2643        #FIXME: These numbers were from version before 21/3/6 -
2644        #could be recreated by setting maximum_allowed_speed to 0 maybe
2645        #assert allclose(domain.quantities['xmomentum'].vertex_values[:4,0],
2646        #                [ 0.00064835, 0.03685719, 0.00085073, 0.03687313])
2647       
2648        assert allclose(domain.quantities['xmomentum'].vertex_values[:4,0],
2649                        [ 0.00090581, 0.03685719, 0.00088303, 0.03687313])
2650                       
2651                       
2652
2653        #assert allclose(domain.quantities['xmomentum'].vertex_values[:4,0],
2654        #                [0.00090581,0.03685719,0.00088303,0.03687313])
2655
2656        assert allclose(domain.quantities['ymomentum'].vertex_values[:4,0],
2657                        [-0.00139463,0.0006156,-0.00060364,0.00061827])
2658
2659
2660        os.remove(domain.get_name() + '.sww')
2661
2662       
2663    def test_flatbed_second_order_vmax_0(self):
2664        from mesh_factory import rectangular
2665        from Numeric import array
2666
2667        #Create basic mesh
2668        N = 8
2669        points, vertices, boundary = rectangular(N, N)
2670
2671        #Create shallow water domain
2672        domain = Domain(points, vertices, boundary)
2673        domain.smooth = False
2674        domain.default_order=2
2675        domain.beta_w      = 0.9
2676        domain.beta_w_dry  = 0.9
2677        domain.beta_uh     = 0.9
2678        domain.beta_uh_dry = 0.9
2679        domain.beta_vh     = 0.9
2680        domain.beta_vh_dry = 0.9       
2681        domain.maximum_allowed_speed = 0.0 #Makes it like the 'oldstyle'
2682        domain.H0 = 0 # Backwards compatibility (6/2/7)       
2683
2684        # Boundary conditions
2685        Br = Reflective_boundary(domain)
2686        Bd = Dirichlet_boundary([0.2,0.,0.])
2687
2688        domain.set_boundary({'left': Bd, 'right': Br, 'top': Br, 'bottom': Br})
2689        domain.check_integrity()
2690
2691        #Evolution
2692        for t in domain.evolve(yieldstep = 0.01, finaltime = 0.03):
2693            pass
2694
2695
2696        assert allclose(domain.min_timestep, 0.0210448446782)
2697        assert allclose(domain.max_timestep, 0.0210448446782)
2698
2699        #FIXME: These numbers were from version before 21/3/6 -
2700        #could be recreated by setting maximum_allowed_speed to 0 maybe
2701        assert allclose(domain.quantities['xmomentum'].vertex_values[:4,0],
2702                        [ 0.00064835, 0.03685719, 0.00085073, 0.03687313]) 
2703       
2704
2705        assert allclose(domain.quantities['ymomentum'].vertex_values[:4,0],
2706                        [-0.00139463,0.0006156,-0.00060364,0.00061827])
2707
2708
2709        os.remove(domain.get_name() + '.sww')
2710
2711       
2712
2713    def test_flatbed_second_order_distribute(self):
2714        #Use real data from anuga.abstract_2d_finite_volumes 2
2715        #painfully setup and extracted.
2716        from mesh_factory import rectangular
2717        from Numeric import array
2718
2719        #Create basic mesh
2720        N = 8
2721        points, vertices, boundary = rectangular(N, N)
2722
2723        #Create shallow water domain
2724        domain = Domain(points, vertices, boundary)
2725        domain.smooth = False
2726        domain.default_order=domain._order_=2
2727        domain.beta_w      = 0.9
2728        domain.beta_w_dry  = 0.9
2729        domain.beta_uh     = 0.9
2730        domain.beta_uh_dry = 0.9
2731        domain.beta_vh     = 0.9
2732        domain.beta_vh_dry = 0.9
2733        domain.H0 = 0 # Backwards compatibility (6/2/7)       
2734
2735        # Boundary conditions
2736        Br = Reflective_boundary(domain)
2737        Bd = Dirichlet_boundary([0.2,0.,0.])
2738
2739        domain.set_boundary({'left': Bd, 'right': Br, 'top': Br, 'bottom': Br})
2740        domain.check_integrity()
2741
2742
2743
2744        for V in [False, True]:
2745            if V:
2746                #Set centroids as if system had been evolved
2747                L = zeros(2*N*N, Float)
2748                L[:32] = [7.21205592e-003, 5.35214298e-002, 1.00910824e-002,
2749                          5.35439433e-002, 1.00910824e-002, 5.35439433e-002,
2750                          1.00910824e-002, 5.35439433e-002, 1.00910824e-002,
2751                          5.35439433e-002, 1.00910824e-002, 5.35439433e-002,
2752                          1.00910824e-002, 5.35393928e-002, 1.02344264e-002,
2753                          5.59605058e-002, 0.00000000e+000, 3.31027800e-004,
2754                          0.00000000e+000, 4.37962142e-005, 0.00000000e+000,
2755                          4.37962142e-005, 0.00000000e+000, 4.37962142e-005,
2756                          0.00000000e+000, 4.37962142e-005, 0.00000000e+000,
2757                          4.37962142e-005, 0.00000000e+000, 4.37962142e-005,
2758                          0.00000000e+000, 5.57305948e-005]
2759
2760                X = zeros(2*N*N, Float)
2761                X[:32] = [6.48351607e-003, 3.68571894e-002, 8.50733285e-003,
2762                          3.68731327e-002, 8.50733285e-003, 3.68731327e-002,
2763                          8.50733285e-003, 3.68731327e-002, 8.50733285e-003,
2764                          3.68731327e-002, 8.50733285e-003, 3.68731327e-002,
2765                          8.50733285e-003, 3.68693861e-002, 8.65220973e-003,
2766                          3.85055387e-002, 0.00000000e+000, 2.86060840e-004,
2767                          0.00000000e+000, 3.58905503e-005, 0.00000000e+000,
2768                          3.58905503e-005, 0.00000000e+000, 3.58905503e-005,
2769                          0.00000000e+000, 3.58905503e-005, 0.00000000e+000,
2770                          3.58905503e-005, 0.00000000e+000, 3.58905503e-005,
2771                          0.00000000e+000, 4.57662812e-005]
2772
2773                Y = zeros(2*N*N, Float)
2774                Y[:32]=[-1.39463104e-003, 6.15600298e-004, -6.03637382e-004,
2775                        6.18272251e-004, -6.03637382e-004, 6.18272251e-004,
2776                        -6.03637382e-004, 6.18272251e-004, -6.03637382e-004,
2777                        6.18272251e-004, -6.03637382e-004, 6.18272251e-004,
2778                        -6.03637382e-004, 6.18599320e-004, -6.74622797e-004,
2779                        -1.48934756e-004, 0.00000000e+000, -5.35079969e-005,
2780                        0.00000000e+000, -2.57264987e-005, 0.00000000e+000,
2781                        -2.57264987e-005, 0.00000000e+000, -2.57264987e-005,
2782                        0.00000000e+000, -2.57264987e-005, 0.00000000e+000,
2783                        -2.57264987e-005, 0.00000000e+000, -2.57264987e-005,
2784                        0.00000000e+000, -2.57635178e-005]
2785
2786
2787                domain.set_quantity('stage', L, location='centroids')
2788                domain.set_quantity('xmomentum', X, location='centroids')
2789                domain.set_quantity('ymomentum', Y, location='centroids')
2790
2791                domain.check_integrity()
2792            else:
2793                #Evolution
2794                for t in domain.evolve(yieldstep = 0.01, finaltime = 0.03):
2795                    pass
2796                assert allclose(domain.min_timestep, 0.0210448446782)
2797                assert allclose(domain.max_timestep, 0.0210448446782)
2798
2799
2800            #Centroids were correct but not vertices.
2801            #Hence the check of distribute below.
2802            assert allclose(domain.quantities['stage'].centroid_values[:4],
2803                            [0.00721206,0.05352143,0.01009108,0.05354394])
2804
2805            assert allclose(domain.quantities['xmomentum'].centroid_values[:4],
2806                            [0.00648352,0.03685719,0.00850733,0.03687313])
2807
2808            assert allclose(domain.quantities['ymomentum'].centroid_values[:4],
2809                            [-0.00139463,0.0006156,-0.00060364,0.00061827])
2810
2811            #print 'C17=', domain.quantities['xmomentum'].centroid_values[17]
2812            #print 'C19=', domain.quantities['xmomentum'].centroid_values[19]
2813
2814            #assert allclose(domain.quantities['xmomentum'].centroid_values[17],0.00028606084)
2815            ##print domain.quantities['xmomentum'].centroid_values[17], V
2816            ##print
2817            if not V:
2818                #FIXME: These numbers were from version before 21/3/6 -
2819                #could be recreated by setting maximum_allowed_speed to 0 maybe           
2820               
2821                #assert allclose(domain.quantities['xmomentum'].centroid_values[17], 0.0)               
2822                assert allclose(domain.quantities['xmomentum'].centroid_values[17], 0.000286060839592)                         
2823
2824            else:
2825                assert allclose(domain.quantities['xmomentum'].centroid_values[17], 0.00028606084)
2826
2827            import copy
2828            XX = copy.copy(domain.quantities['xmomentum'].centroid_values)
2829            assert allclose(domain.quantities['xmomentum'].centroid_values, XX)
2830
2831            domain.distribute_to_vertices_and_edges()
2832
2833            #assert allclose(domain.quantities['xmomentum'].centroid_values, XX)
2834
2835            #assert allclose(domain.quantities['xmomentum'].centroid_values[17],
2836            #                0.0)
2837            assert allclose(domain.quantities['xmomentum'].centroid_values[17], 0.000286060839592)                                                 
2838
2839
2840            #FIXME: These numbers were from version before 25/10
2841            #assert allclose(domain.quantities['stage'].vertex_values[:4,0],
2842            #                [0.00101913,0.05352143,0.00104852,0.05354394])
2843
2844            assert allclose(domain.quantities['ymomentum'].vertex_values[:4,0],
2845                            [-0.00139463,0.0006156,-0.00060364,0.00061827])
2846
2847
2848            assert allclose(domain.quantities['xmomentum'].vertex_values[:4,0],
2849                            [0.00090581,0.03685719,0.00088303,0.03687313])
2850
2851
2852            #NB NO longer relvant:
2853
2854            #This was the culprit. First triangles vertex 0 had an
2855            #x-momentum of 0.0064835 instead of 0.00090581 and
2856            #third triangle had 0.00850733 instead of 0.00088303
2857            #print domain.quantities['xmomentum'].vertex_values[:4,0]
2858
2859            #print domain.quantities['xmomentum'].vertex_values[:4,0]
2860            #assert allclose(domain.quantities['xmomentum'].vertex_values[:4,0],
2861            #                [0.00090581,0.03685719,0.00088303,0.03687313])
2862
2863        os.remove(domain.get_name() + '.sww')
2864
2865
2866
2867    def test_bedslope_problem_first_order(self):
2868
2869        from mesh_factory import rectangular
2870        from Numeric import array
2871
2872        #Create basic mesh
2873        points, vertices, boundary = rectangular(6, 6)
2874
2875        #Create shallow water domain
2876        domain = Domain(points, vertices, boundary)
2877        domain.smooth = False
2878        domain.default_order=1
2879
2880        #Bed-slope and friction
2881        def x_slope(x, y):
2882            return -x/3
2883
2884        domain.set_quantity('elevation', x_slope)
2885
2886        # Boundary conditions
2887        Br = Reflective_boundary(domain)
2888        domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
2889
2890        #Initial condition
2891        domain.set_quantity('stage', Constant_height(x_slope, 0.05))
2892        domain.check_integrity()
2893
2894        #Evolution
2895        for t in domain.evolve(yieldstep = 0.05, finaltime = 0.05):
2896            pass# domain.write_time()
2897
2898        assert allclose(domain.min_timestep, 0.050010003001)
2899        assert allclose(domain.max_timestep, 0.050010003001)
2900
2901
2902        os.remove(domain.get_name() + '.sww')
2903       
2904    def test_bedslope_problem_first_order_moresteps(self):
2905
2906        from mesh_factory import rectangular
2907        from Numeric import array
2908
2909        #Create basic mesh
2910        points, vertices, boundary = rectangular(6, 6)
2911
2912        #Create shallow water domain
2913        domain = Domain(points, vertices, boundary)
2914        domain.smooth = False
2915        domain.default_order=1
2916        domain.beta_h = 0.0 #Use first order in h-limiter
2917        domain.H0 = 0 # Backwards compatibility (6/2/7)       
2918
2919        #Bed-slope and friction
2920        def x_slope(x, y):
2921            return -x/3
2922
2923        domain.set_quantity('elevation', x_slope)
2924
2925        # Boundary conditions
2926        Br = Reflective_boundary(domain)
2927        domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
2928
2929        #Initial condition
2930        domain.set_quantity('stage', Constant_height(x_slope, 0.05))
2931        domain.check_integrity()
2932
2933        #Evolution
2934        for t in domain.evolve(yieldstep = 0.05, finaltime = 0.5):
2935            pass# domain.write_time()
2936
2937        #Data from earlier version of abstract_2d_finite_volumes
2938        #print domain.quantities['stage'].centroid_values
2939
2940        assert allclose(domain.quantities['stage'].centroid_values,
2941                        [-0.02998628, -0.01520652, -0.03043492,
2942                        -0.0149132, -0.03004706, -0.01476251,
2943                        -0.0298215, -0.01467976, -0.02988158,
2944                        -0.01474662, -0.03036161, -0.01442995,
2945                        -0.07624583, -0.06297061, -0.07733792,
2946                        -0.06342237, -0.07695439, -0.06289595,
2947                        -0.07635559, -0.0626065, -0.07633628,
2948                        -0.06280072, -0.07739632, -0.06386738,
2949                        -0.12161738, -0.11028239, -0.1223796,
2950                        -0.11095953, -0.12189744, -0.11048616,
2951                        -0.12074535, -0.10987605, -0.12014311,
2952                        -0.10976691, -0.12096859, -0.11087692,
2953                        -0.16868259, -0.15868061, -0.16801135,
2954                        -0.1588003, -0.16674343, -0.15813323,
2955                        -0.16457595, -0.15693826, -0.16281096,
2956                        -0.15585154, -0.16283873, -0.15540068,
2957                        -0.17450362, -0.19919913, -0.18062882,
2958                        -0.19764131, -0.17783111, -0.19407213,
2959                        -0.1736915, -0.19053624, -0.17228678,
2960                        -0.19105634, -0.17920133, -0.1968828,
2961                        -0.14244395, -0.14604641, -0.14473537,
2962                        -0.1506107, -0.14510055, -0.14919522,
2963                        -0.14175896, -0.14560798, -0.13911658,
2964                        -0.14439383, -0.13924047, -0.14829043])
2965
2966        os.remove(domain.get_name() + '.sww')
2967       
2968    def test_bedslope_problem_second_order_one_step(self):
2969
2970        from mesh_factory import rectangular
2971        from Numeric import array
2972
2973        #Create basic mesh
2974        points, vertices, boundary = rectangular(6, 6)
2975
2976        #Create shallow water domain
2977        domain = Domain(points, vertices, boundary)
2978        domain.smooth = False
2979        domain.default_order=2
2980        domain.beta_w      = 0.9
2981        domain.beta_w_dry  = 0.9
2982        domain.beta_uh     = 0.9
2983        domain.beta_uh_dry = 0.9
2984        domain.beta_vh     = 0.9
2985        domain.beta_vh_dry = 0.9
2986
2987        #Bed-slope and friction at vertices (and interpolated elsewhere)
2988        def x_slope(x, y):
2989            return -x/3
2990
2991        domain.set_quantity('elevation', x_slope)
2992
2993        # Boundary conditions
2994        Br = Reflective_boundary(domain)
2995        domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
2996
2997        #Initial condition
2998        domain.set_quantity('stage', Constant_height(x_slope, 0.05))
2999        domain.check_integrity()
3000
3001        assert allclose(domain.quantities['stage'].centroid_values,
3002                        [0.01296296, 0.03148148, 0.01296296,
3003                        0.03148148, 0.01296296, 0.03148148,
3004                        0.01296296, 0.03148148, 0.01296296,
3005                        0.03148148, 0.01296296, 0.03148148,
3006                        -0.04259259, -0.02407407, -0.04259259,
3007                        -0.02407407, -0.04259259, -0.02407407,
3008                        -0.04259259, -0.02407407, -0.04259259,
3009                        -0.02407407, -0.04259259, -0.02407407,
3010                        -0.09814815, -0.07962963, -0.09814815,
3011                        -0.07962963, -0.09814815, -0.07962963,
3012                        -0.09814815, -0.07962963, -0.09814815,
3013                        -0.07962963, -0.09814815, -0.07962963,
3014                        -0.1537037 , -0.13518519, -0.1537037,
3015                        -0.13518519, -0.1537037, -0.13518519,
3016                        -0.1537037 , -0.13518519, -0.1537037,
3017                        -0.13518519, -0.1537037, -0.13518519,
3018                        -0.20925926, -0.19074074, -0.20925926,
3019                        -0.19074074, -0.20925926, -0.19074074,
3020                        -0.20925926, -0.19074074, -0.20925926,
3021                        -0.19074074, -0.20925926, -0.19074074,
3022                        -0.26481481, -0.2462963, -0.26481481,
3023                        -0.2462963, -0.26481481, -0.2462963,
3024                        -0.26481481, -0.2462963, -0.26481481,
3025                        -0.2462963, -0.26481481, -0.2462963])
3026
3027
3028        #print domain.quantities['stage'].extrapolate_second_order()
3029        #domain.distribute_to_vertices_and_edges()
3030        #print domain.quantities['stage'].vertex_values[:,0]
3031
3032        #Evolution
3033        for t in domain.evolve(yieldstep = 0.05, finaltime = 0.05):
3034            #domain.write_time()
3035            pass
3036
3037
3038        #print domain.quantities['stage'].centroid_values
3039        assert allclose(domain.quantities['stage'].centroid_values,
3040                 [0.01290985, 0.02356019, 0.01619096, 0.02356019, 0.01619096,
3041                  0.02356019, 0.01619096, 0.02356019, 0.01619096, 0.02356019,
3042                  0.01619096, 0.0268413, -0.04411074, -0.0248011, -0.04186556,
3043                  -0.0248011, -0.04186556, -0.0248011, -0.04186556, -0.0248011,
3044                  -0.04186556, -0.0248011, -0.04186556, -0.02255593,
3045                  -0.09966629, -0.08035666, -0.09742112, -0.08035666,
3046                  -0.09742112, -0.08035666, -0.09742112, -0.08035666,
3047                  -0.09742112, -0.08035666, -0.09742112, -0.07811149,
3048                  -0.15522185, -0.13591222, -0.15297667, -0.13591222,
3049                  -0.15297667, -0.13591222, -0.15297667, -0.13591222,
3050                  -0.15297667, -0.13591222, -0.15297667, -0.13366704,
3051                  -0.2107774, -0.19146777, -0.20853223, -0.19146777,
3052                  -0.20853223, -0.19146777, -0.20853223, -0.19146777,
3053                  -0.20853223, -0.19146777, -0.20853223, -0.1892226,
3054                  -0.26120669, -0.24776246, -0.25865535, -0.24776246,
3055                  -0.25865535, -0.24776246, -0.25865535, -0.24776246,
3056                  -0.25865535, -0.24776246, -0.25865535, -0.24521113])
3057
3058        os.remove(domain.get_name() + '.sww')
3059
3060    def test_bedslope_problem_second_order_two_steps(self):
3061
3062        from mesh_factory import rectangular
3063        from Numeric import array
3064
3065        #Create basic mesh
3066        points, vertices, boundary = rectangular(6, 6)
3067
3068        #Create shallow water domain
3069        domain = Domain(points, vertices, boundary)
3070        domain.smooth = False
3071        domain.default_order=2
3072        domain.beta_w      = 0.9
3073        domain.beta_w_dry  = 0.9
3074        domain.beta_uh     = 0.9
3075        domain.beta_uh_dry = 0.9
3076        domain.beta_vh     = 0.9
3077        domain.beta_vh_dry = 0.9
3078        domain.beta_h = 0.0 #Use first order in h-limiter
3079        domain.H0 = 0 # Backwards compatibility (6/2/7)       
3080
3081        #Bed-slope and friction at vertices (and interpolated elsewhere)
3082        def x_slope(x, y):
3083            return -x/3
3084
3085        domain.set_quantity('elevation', x_slope)
3086
3087        # Boundary conditions
3088        Br = Reflective_boundary(domain)
3089        domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
3090
3091        #Initial condition
3092        domain.set_quantity('stage', Constant_height(x_slope, 0.05))
3093        domain.check_integrity()
3094
3095        assert allclose(domain.quantities['stage'].centroid_values,
3096                        [0.01296296, 0.03148148, 0.01296296,
3097                        0.03148148, 0.01296296, 0.03148148,
3098                        0.01296296, 0.03148148, 0.01296296,
3099                        0.03148148, 0.01296296, 0.03148148,
3100                        -0.04259259, -0.02407407, -0.04259259,
3101                        -0.02407407, -0.04259259, -0.02407407,
3102                        -0.04259259, -0.02407407, -0.04259259,
3103                        -0.02407407, -0.04259259, -0.02407407,
3104                        -0.09814815, -0.07962963, -0.09814815,
3105                        -0.07962963, -0.09814815, -0.07962963,
3106                        -0.09814815, -0.07962963, -0.09814815,
3107                        -0.07962963, -0.09814815, -0.07962963,
3108                        -0.1537037 , -0.13518519, -0.1537037,
3109                        -0.13518519, -0.1537037, -0.13518519,
3110                        -0.1537037 , -0.13518519, -0.1537037,
3111                        -0.13518519, -0.1537037, -0.13518519,
3112                        -0.20925926, -0.19074074, -0.20925926,
3113                        -0.19074074, -0.20925926, -0.19074074,
3114                        -0.20925926, -0.19074074, -0.20925926,
3115                        -0.19074074, -0.20925926, -0.19074074,
3116                        -0.26481481, -0.2462963, -0.26481481,
3117                        -0.2462963, -0.26481481, -0.2462963,
3118                        -0.26481481, -0.2462963, -0.26481481,
3119                        -0.2462963, -0.26481481, -0.2462963])
3120
3121
3122        #print domain.quantities['stage'].extrapolate_second_order()
3123        #domain.distribute_to_vertices_and_edges()
3124        #print domain.quantities['stage'].vertex_values[:,0]
3125
3126        #Evolution
3127        for t in domain.evolve(yieldstep = 0.05, finaltime = 0.1):
3128            pass
3129
3130
3131        #Data from earlier version of abstract_2d_finite_volumes ft=0.1
3132        assert allclose(domain.min_timestep, 0.0376895634803)
3133        assert allclose(domain.max_timestep, 0.0415635655309)
3134
3135
3136        assert allclose(domain.quantities['stage'].centroid_values,
3137                        [0.00855788, 0.01575204, 0.00994606, 0.01717072,
3138                         0.01005985, 0.01716362, 0.01005985, 0.01716299,
3139                         0.01007098, 0.01736248, 0.01216452, 0.02026776,
3140                         -0.04462374, -0.02479045, -0.04199789, -0.0229465,
3141                         -0.04184033, -0.02295693, -0.04184013, -0.02295675,
3142                         -0.04184486, -0.0228168, -0.04028876, -0.02036486,
3143                         -0.10029444, -0.08170809, -0.09772846, -0.08021704,
3144                         -0.09760006, -0.08022143, -0.09759984, -0.08022124,
3145                         -0.09760261, -0.08008893, -0.09603914, -0.07758209,
3146                         -0.15584152, -0.13723138, -0.15327266, -0.13572906,
3147                         -0.15314427, -0.13573349, -0.15314405, -0.13573331,
3148                         -0.15314679, -0.13560104, -0.15158523, -0.13310701,
3149                         -0.21208605, -0.19283913, -0.20955631, -0.19134189,
3150                         -0.20942821, -0.19134598, -0.20942799, -0.1913458,
3151                         -0.20943005, -0.19120952, -0.20781177, -0.18869401,
3152                         -0.25384082, -0.2463294, -0.25047649, -0.24464654,
3153                         -0.25031159, -0.24464253, -0.25031112, -0.24464253,
3154                         -0.25031463, -0.24454764, -0.24885323, -0.24286438])
3155
3156
3157        os.remove(domain.get_name() + '.sww')
3158
3159    def test_bedslope_problem_second_order_two_yieldsteps(self):
3160
3161        from mesh_factory import rectangular
3162        from Numeric import array
3163
3164        #Create basic mesh
3165        points, vertices, boundary = rectangular(6, 6)
3166
3167        #Create shallow water domain
3168        domain = Domain(points, vertices, boundary)
3169        domain.smooth = False
3170        domain.default_order=2
3171        domain.beta_w      = 0.9
3172        domain.beta_w_dry  = 0.9
3173        domain.beta_uh     = 0.9
3174        domain.beta_uh_dry = 0.9
3175        domain.beta_vh     = 0.9
3176        domain.beta_vh_dry = 0.9
3177        domain.beta_h = 0.0 #Use first order in h-limiter
3178        domain.H0 = 0 # Backwards compatibility (6/2/7)       
3179
3180        #Bed-slope and friction at vertices (and interpolated elsewhere)
3181        def x_slope(x, y):
3182            return -x/3
3183
3184        domain.set_quantity('elevation', x_slope)
3185
3186        # Boundary conditions
3187        Br = Reflective_boundary(domain)
3188        domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
3189
3190        #Initial condition
3191        domain.set_quantity('stage', Constant_height(x_slope, 0.05))
3192        domain.check_integrity()
3193
3194        assert allclose(domain.quantities['stage'].centroid_values,
3195                        [0.01296296, 0.03148148, 0.01296296,
3196                        0.03148148, 0.01296296, 0.03148148,
3197                        0.01296296, 0.03148148, 0.01296296,
3198                        0.03148148, 0.01296296, 0.03148148,
3199                        -0.04259259, -0.02407407, -0.04259259,
3200                        -0.02407407, -0.04259259, -0.02407407,
3201                        -0.04259259, -0.02407407, -0.04259259,
3202                        -0.02407407, -0.04259259, -0.02407407,
3203                        -0.09814815, -0.07962963, -0.09814815,
3204                        -0.07962963, -0.09814815, -0.07962963,
3205                        -0.09814815, -0.07962963, -0.09814815,
3206                        -0.07962963, -0.09814815, -0.07962963,
3207                        -0.1537037 , -0.13518519, -0.1537037,
3208                        -0.13518519, -0.1537037, -0.13518519,
3209                        -0.1537037 , -0.13518519, -0.1537037,
3210                        -0.13518519, -0.1537037, -0.13518519,
3211                        -0.20925926, -0.19074074, -0.20925926,
3212                        -0.19074074, -0.20925926, -0.19074074,
3213                        -0.20925926, -0.19074074, -0.20925926,
3214                        -0.19074074, -0.20925926, -0.19074074,
3215                        -0.26481481, -0.2462963, -0.26481481,
3216                        -0.2462963, -0.26481481, -0.2462963,
3217                        -0.26481481, -0.2462963, -0.26481481,
3218                        -0.2462963, -0.26481481, -0.2462963])
3219
3220
3221        #print domain.quantities['stage'].extrapolate_second_order()
3222        #domain.distribute_to_vertices_and_edges()
3223        #print domain.quantities['stage'].vertex_values[:,0]
3224
3225        #Evolution
3226        for t in domain.evolve(yieldstep = 0.05, finaltime = 0.1):   #0.05??
3227            #domain.write_time()
3228            pass
3229
3230
3231
3232        assert allclose(domain.quantities['stage'].centroid_values,
3233                 [0.00855788, 0.01575204, 0.00994606, 0.01717072, 0.01005985,
3234                  0.01716362, 0.01005985, 0.01716299, 0.01007098, 0.01736248,
3235                  0.01216452, 0.02026776, -0.04462374, -0.02479045, -0.04199789,
3236                  -0.0229465, -0.04184033, -0.02295693, -0.04184013,
3237                  -0.02295675, -0.04184486, -0.0228168, -0.04028876,
3238                  -0.02036486, -0.10029444, -0.08170809, -0.09772846,
3239                  -0.08021704, -0.09760006, -0.08022143, -0.09759984,
3240                  -0.08022124, -0.09760261, -0.08008893, -0.09603914,
3241                  -0.07758209, -0.15584152, -0.13723138, -0.15327266,
3242                  -0.13572906, -0.15314427, -0.13573349, -0.15314405,
3243                  -0.13573331, -0.15314679, -0.13560104, -0.15158523,
3244                  -0.13310701, -0.21208605, -0.19283913, -0.20955631,
3245                  -0.19134189, -0.20942821, -0.19134598, -0.20942799,
3246                  -0.1913458, -0.20943005, -0.19120952, -0.20781177,
3247                  -0.18869401, -0.25384082, -0.2463294, -0.25047649,
3248                  -0.24464654, -0.25031159, -0.24464253, -0.25031112,
3249                  -0.24464253, -0.25031463, -0.24454764, -0.24885323,
3250                  -0.24286438])
3251
3252        os.remove(domain.get_name() + '.sww')
3253
3254    def test_bedslope_problem_second_order_more_steps(self):
3255
3256        from mesh_factory import rectangular
3257        from Numeric import array
3258
3259        #Create basic mesh
3260        points, vertices, boundary = rectangular(6, 6)
3261
3262        #Create shallow water domain
3263        domain = Domain(points, vertices, boundary)
3264        domain.smooth = False
3265        domain.default_order=2
3266        domain.beta_w      = 0.9
3267        domain.beta_w_dry  = 0.9
3268        domain.beta_uh     = 0.9
3269        domain.beta_uh_dry = 0.9
3270        domain.beta_vh     = 0.9
3271        domain.beta_vh_dry = 0.9
3272        domain.beta_h = 0.0 #Use first order in h-limiter
3273        domain.H0 = 0 # Backwards compatibility (6/2/7)       
3274
3275        #Bed-slope and friction at vertices (and interpolated elsewhere)
3276        def x_slope(x, y):
3277            return -x/3
3278
3279        domain.set_quantity('elevation', x_slope)
3280
3281        # Boundary conditions
3282        Br = Reflective_boundary(domain)
3283        domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
3284
3285        #Initial condition
3286        domain.set_quantity('stage', expression = 'elevation + 0.05')
3287        domain.check_integrity()
3288
3289        assert allclose(domain.quantities['stage'].centroid_values,
3290                        [0.01296296, 0.03148148, 0.01296296,
3291                        0.03148148, 0.01296296, 0.03148148,
3292                        0.01296296, 0.03148148, 0.01296296,
3293                        0.03148148, 0.01296296, 0.03148148,
3294                        -0.04259259, -0.02407407, -0.04259259,
3295                        -0.02407407, -0.04259259, -0.02407407,
3296                        -0.04259259, -0.02407407, -0.04259259,
3297                        -0.02407407, -0.04259259, -0.02407407,
3298                        -0.09814815, -0.07962963, -0.09814815,
3299                        -0.07962963, -0.09814815, -0.07962963,
3300                        -0.09814815, -0.07962963, -0.09814815,
3301                        -0.07962963, -0.09814815, -0.07962963,
3302                        -0.1537037 , -0.13518519, -0.1537037,
3303                        -0.13518519, -0.1537037, -0.13518519,
3304                        -0.1537037 , -0.13518519, -0.1537037,
3305                        -0.13518519, -0.1537037, -0.13518519,
3306                        -0.20925926, -0.19074074, -0.20925926,
3307                        -0.19074074, -0.20925926, -0.19074074,
3308                        -0.20925926, -0.19074074, -0.20925926,
3309                        -0.19074074, -0.20925926, -0.19074074,
3310                        -0.26481481, -0.2462963, -0.26481481,
3311                        -0.2462963, -0.26481481, -0.2462963,
3312                        -0.26481481, -0.2462963, -0.26481481,
3313                        -0.2462963, -0.26481481, -0.2462963])
3314
3315
3316        #print domain.quantities['stage'].extrapolate_second_order()
3317        #domain.distribute_to_vertices_and_edges()
3318        #print domain.quantities['stage'].vertex_values[:,0]
3319
3320        #Evolution
3321        for t in domain.evolve(yieldstep = 0.05, finaltime = 0.5):
3322            pass
3323
3324
3325        assert allclose(domain.quantities['stage'].centroid_values,
3326     [-0.02907028, -0.01475478, -0.02973417, -0.01447186, -0.02932665, -0.01428285,
3327      -0.02901975, -0.0141361,  -0.02898816, -0.01418135, -0.02961409, -0.01403487,
3328      -0.07597998, -0.06252591, -0.07664854, -0.06312532, -0.07638287, -0.06265139,
3329      -0.07571145, -0.06235231, -0.0756817,  -0.06245309, -0.07652292, -0.06289946,
3330      -0.12367464, -0.11088981, -0.12237277, -0.11115338, -0.1218934,  -0.1107174,
3331      -0.12081485, -0.11000491, -0.12038451, -0.11010335, -0.12102113, -0.11012105,
3332      -0.16909116, -0.15831543, -0.16730214, -0.15786249, -0.1665493,  -0.15697919,
3333      -0.16496618, -0.15559852, -0.16338679, -0.15509088, -0.16364092, -0.15424423,
3334      -0.18771107, -0.19903904, -0.18903759, -0.19858437, -0.18701552, -0.19697797,
3335      -0.1833593,  -0.19505871, -0.1818806,  -0.19418042, -0.18586159, -0.19576946,
3336      -0.13986873, -0.14170053, -0.14132188, -0.14560674, -0.14095617, -0.14373292,
3337      -0.13785933, -0.14033364, -0.13592955, -0.13936356, -0.13596008, -0.14216296])
3338
3339        assert allclose(domain.quantities['xmomentum'].centroid_values,
3340     [ 0.00831121,  0.00317948,  0.00731797,  0.00334939,  0.00764717,  0.00348053,
3341       0.00788729,  0.00356522,  0.00780649,  0.00341919,  0.00693525,  0.00310375,
3342       0.02166196,  0.01421475,  0.02017737,  0.01316839,  0.02037015,  0.01368659,
3343       0.02106,     0.01399161,  0.02074514,  0.01354935,  0.01887407,  0.0123113,
3344       0.03775083,  0.02855197,  0.03689337,  0.02759782,  0.03732848,  0.02812072,
3345       0.03872545,  0.02913348,  0.03880939,  0.02803804,  0.03546499,  0.0260039,
3346       0.0632131,   0.04730634,  0.0576324,   0.04592336,  0.05790921,  0.04690514,
3347       0.05986467,  0.04871165,  0.06170068,  0.04811572,  0.05657041,  0.04416292,
3348       0.08489642,  0.07188097,  0.07835261,  0.06843406,  0.07986412,  0.0698247,
3349       0.08201071,  0.07216756,  0.08378418,  0.07273624,  0.080399,    0.06645841,
3350       0.01631548,  0.04691608,  0.0206632,   0.044409,    0.02115518,  0.04560305,
3351       0.02160608,  0.04663725,  0.02174734,  0.04795559,  0.02281427,  0.05667111])
3352
3353
3354        assert allclose(domain.quantities['ymomentum'].centroid_values,
3355     [ 1.45876601e-004, -3.24627393e-004, -1.57572719e-004, -2.92790187e-004,
3356      -9.90988382e-005, -3.06677335e-004, -1.62493106e-004, -3.71310004e-004,
3357      -1.99445058e-004, -3.28493467e-004,  6.68217349e-005, -8.42042805e-006,
3358       5.05093371e-004, -1.42842214e-004, -6.81454718e-005, -5.02084057e-004,
3359      -8.50583861e-005, -4.65443981e-004, -1.96406564e-004, -5.88889562e-004,
3360      -2.70160173e-004, -5.35485454e-004,  2.60780997e-004,  3.12145471e-005,
3361       5.16189608e-004,  1.07069062e-004,  9.29989252e-005, -3.71211119e-004,
3362       1.16350246e-004, -3.82407830e-004, -1.62077969e-004, -6.30906636e-004,
3363      -4.74025708e-004, -6.94463009e-004,  6.15092843e-005,  2.22106820e-004,
3364      -6.29589294e-004,  2.43611937e-004, -5.88125094e-004, -6.94293192e-005,
3365      -4.17914641e-004,  6.64609019e-005, -7.68334577e-004, -3.40232101e-004,
3366      -1.67424308e-003, -7.39485066e-004, -1.59966988e-003,  5.68262838e-005,
3367      -1.48470633e-003, -1.84554882e-003, -2.27200099e-003, -1.67506848e-003,
3368      -1.95610258e-003, -1.47638801e-003, -1.73779477e-003, -1.85498791e-003,
3369      -2.01357843e-003, -2.17675471e-003, -1.65783870e-003, -1.15818681e-003,
3370      -1.18663036e-003, -2.94229849e-003, -3.59309018e-003, -5.13496584e-003,
3371      -6.17359400e-003, -5.98761937e-003, -6.00540116e-003, -5.01121966e-003,
3372      -4.50964850e-003, -3.06319963e-003,  6.08950810e-004, -4.79537921e-004])
3373
3374        os.remove(domain.get_name() + '.sww')
3375
3376
3377
3378    def test_bedslope_problem_second_order_more_steps_feb_2007(self):
3379        """test_bedslope_problem_second_order_more_steps_feb_2007
3380
3381        Test shallow water finite volumes, using parameters from feb 2007 rather
3382        than backward compatibility ad infinitum
3383       
3384        """
3385        from mesh_factory import rectangular
3386        from Numeric import array
3387
3388        #Create basic mesh
3389        points, vertices, boundary = rectangular(6, 6)
3390
3391        #Create shallow water domain
3392        domain = Domain(points, vertices, boundary)
3393        domain.smooth = False
3394        domain.default_order = 2
3395        domain.beta_w      = 0.9
3396        domain.beta_w_dry  = 0.9
3397        domain.beta_uh     = 0.9
3398        domain.beta_uh_dry = 0.9
3399        domain.beta_vh     = 0.9
3400        domain.beta_vh_dry = 0.9
3401        domain.beta_h = 0.0 #Use first order in h-limiter
3402        domain.H0 = 0.001 
3403
3404        #Bed-slope and friction at vertices (and interpolated elsewhere)
3405        def x_slope(x, y):
3406            return -x/3
3407
3408        domain.set_quantity('elevation', x_slope)
3409
3410        # Boundary conditions
3411        Br = Reflective_boundary(domain)
3412        domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
3413
3414        #Initial condition
3415        domain.set_quantity('stage', expression = 'elevation + 0.05')
3416        domain.check_integrity()
3417
3418        assert allclose(domain.quantities['stage'].centroid_values,
3419                        [0.01296296, 0.03148148, 0.01296296,
3420                        0.03148148, 0.01296296, 0.03148148,
3421                        0.01296296, 0.03148148, 0.01296296,
3422                        0.03148148, 0.01296296, 0.03148148,
3423                        -0.04259259, -0.02407407, -0.04259259,
3424                        -0.02407407, -0.04259259, -0.02407407,
3425                        -0.04259259, -0.02407407, -0.04259259,
3426                        -0.02407407, -0.04259259, -0.02407407,
3427                        -0.09814815, -0.07962963, -0.09814815,
3428                        -0.07962963, -0.09814815, -0.07962963,
3429                        -0.09814815, -0.07962963, -0.09814815,
3430                        -0.07962963, -0.09814815, -0.07962963,
3431                        -0.1537037 , -0.13518519, -0.1537037,
3432                        -0.13518519, -0.1537037, -0.13518519,
3433                        -0.1537037 , -0.13518519, -0.1537037,
3434                        -0.13518519, -0.1537037, -0.13518519,
3435                        -0.20925926, -0.19074074, -0.20925926,
3436                        -0.19074074, -0.20925926, -0.19074074,
3437                        -0.20925926, -0.19074074, -0.20925926,
3438                        -0.19074074, -0.20925926, -0.19074074,
3439                        -0.26481481, -0.2462963, -0.26481481,
3440                        -0.2462963, -0.26481481, -0.2462963,
3441                        -0.26481481, -0.2462963, -0.26481481,
3442                        -0.2462963, -0.26481481, -0.2462963])
3443
3444
3445        #print domain.quantities['stage'].extrapolate_second_order()
3446        #domain.distribute_to_vertices_and_edges()
3447        #print domain.quantities['stage'].vertex_values[:,0]
3448
3449        #Evolution
3450        for t in domain.evolve(yieldstep = 0.05, finaltime = 0.5):
3451            pass
3452
3453
3454        assert allclose(domain.quantities['stage'].centroid_values,
3455     [-0.02907614, -0.01462948, -0.02971293, -0.01435568, -0.02930822, -0.0141715,
3456      -0.02900256, -0.01402774, -0.02897007, -0.01407298, -0.02958336, -0.01393198,
3457      -0.07599907, -0.06253965, -0.07666738, -0.06312924, -0.07639753, -0.06265574,
3458      -0.07572658, -0.06235626, -0.07569544, -0.0624551,  -0.07653298, -0.06289883,
3459      -0.1236842,  -0.11090096, -0.12238737, -0.11116539, -0.12190705, -0.11072785,
3460      -0.1208281,  -0.11001521, -0.12039674, -0.11011275, -0.1210331,  -0.11013222,
3461      -0.16910086, -0.1583269,  -0.16731364, -0.15787261, -0.16655826, -0.15698876,
3462      -0.16497539, -0.15560667, -0.16339389, -0.15509626, -0.16364638, -0.15425332,
3463      -0.18773952, -0.19904678, -0.18906095, -0.1985912,  -0.18703781, -0.19698407,
3464      -0.18337965, -0.19506463, -0.18190047, -0.19418413, -0.18587491, -0.19576587,
3465      -0.1398922,  -0.14172812, -0.14134412, -0.14563187, -0.14097668, -0.14375562,
3466      -0.13787839, -0.14035428, -0.13594691, -0.13938278, -0.13597595, -0.14218274])
3467                     
3468
3469        assert allclose(domain.quantities['xmomentum'].centroid_values,
3470     [ 0.00831991,  0.00327584,  0.0073476,   0.0034367,   0.00767331,  0.00356348,
3471       0.00791201,  0.00364529,  0.00783242,  0.00349728,  0.0069742,   0.0031689,
3472       0.02165378,  0.01421696,  0.02016895,  0.01318091,  0.02036575,  0.01369801,
3473       0.02105492,  0.01400313,  0.02074176,  0.01356261,  0.01887467,  0.01232703,
3474       0.03774493,  0.02854758,  0.03688063,  0.02759182,  0.03731697,  0.02811662,
3475       0.03871488,  0.02912926,  0.03880077,  0.02803529,  0.0354565,   0.02599893,
3476       0.06319926,  0.04729427,  0.05761625,  0.04591713,  0.05789906,  0.0468986,
3477       0.05985286,  0.04870638,  0.06169141,  0.04811799,  0.05656696,  0.04415568,
3478       0.08488718,  0.07187458,  0.07834844,  0.06842993,  0.07985979,  0.06981954,
3479       0.08200942,  0.07216429,  0.08378261,  0.07273359,  0.08040488,  0.06646477,
3480       0.01631518,  0.04691654,  0.02066439,  0.04441294,  0.02115705,  0.04560776,
3481       0.02160783,  0.04664204,  0.02174952,  0.0479616,   0.02281756,  0.05667927])
3482
3483
3484        assert allclose(domain.quantities['ymomentum'].centroid_values,
3485     [  1.49908296e-04,  -3.32118806e-04,  -1.55139120e-04,  -2.97772609e-04,
3486       -9.57477241e-05,  -3.11011790e-04,  -1.58896911e-04,  -3.76997605e-04,
3487       -1.97659519e-04,  -3.34831296e-04,   6.54935308e-05,  -8.43493883e-06,
3488        5.05063334e-04,  -1.43305846e-04,  -6.76061638e-05,  -5.01728590e-04,
3489       -8.39003270e-05,  -4.64804117e-04,  -1.95135035e-04,  -5.88384357e-04,
3490       -2.69800068e-04,  -5.35718006e-04,   2.59588334e-04,   3.00642498e-05,
3491        5.15520349e-04,   1.05711270e-04,   9.26087960e-05,  -3.71855339e-04,
3492        1.16386690e-04,  -3.82345749e-04,  -1.61741416e-04,  -6.31090292e-04,
3493       -4.74530925e-04,  -6.95229436e-04,   6.08967544e-05,   2.20974020e-04,
3494       -6.30000309e-04,   2.42320158e-04,  -5.89290588e-04,  -7.03283441e-05,
3495       -4.18421888e-04,   6.62703090e-05,  -7.68647846e-04,  -3.40294284e-04,
3496       -1.67585557e-03,  -7.40500723e-04,  -1.60020305e-03,   5.62070746e-05,
3497       -1.48807206e-03,  -1.84455791e-03,  -2.27365819e-03,  -1.67381169e-03,
3498       -1.95607481e-03,  -1.47497442e-03,  -1.73851968e-03,  -1.85314716e-03,
3499       -2.01489344e-03,  -2.17608206e-03,  -1.66072261e-03,  -1.15856505e-03,
3500       -1.18717624e-03,  -2.94595857e-03,  -3.59685615e-03,  -5.13811671e-03,
3501       -6.17481232e-03,  -5.98871894e-03,  -6.00593324e-03,  -5.01282532e-03,
3502       -4.51124363e-03,  -3.06579417e-03,   6.07680580e-04,  -4.80786234e-04])
3503
3504        os.remove(domain.get_name() + '.sww')
3505
3506
3507    def test_temp_play(self):
3508
3509        from mesh_factory import rectangular
3510        from Numeric import array
3511
3512        #Create basic mesh
3513        points, vertices, boundary = rectangular(5, 5)
3514
3515        #Create shallow water domain
3516        domain = Domain(points, vertices, boundary)
3517        domain.smooth = False
3518        domain.default_order=2
3519        domain.beta_w      = 0.9
3520        domain.beta_w_dry  = 0.9
3521        domain.beta_uh     = 0.9
3522        domain.beta_uh_dry = 0.9
3523        domain.beta_vh     = 0.9
3524        domain.beta_vh_dry = 0.9
3525        domain.beta_h = 0.0 #Use first order in h-limiter
3526        domain.H0 = 0 # Backwards compatibility (6/2/7)
3527       
3528
3529        #Bed-slope and friction at vertices (and interpolated elsewhere)
3530        def x_slope(x, y):
3531            return -x/3
3532
3533        domain.set_quantity('elevation', x_slope)
3534
3535        # Boundary conditions
3536        Br = Reflective_boundary(domain)
3537        domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
3538
3539        #Initial condition
3540        domain.set_quantity('stage', Constant_height(x_slope, 0.05))
3541        domain.check_integrity()
3542
3543        #Evolution
3544        for t in domain.evolve(yieldstep = 0.05, finaltime = 0.1):
3545            pass
3546
3547        assert allclose(domain.quantities['stage'].centroid_values[:4],
3548                        [0.00206836, 0.01296714, 0.00363415, 0.01438924])
3549        #print domain.quantities['xmomentum'].centroid_values[:4]
3550        assert allclose(domain.quantities['xmomentum'].centroid_values[:4],
3551                        [0.01360154, 0.00671133, 0.01264578, 0.00648503])
3552        assert allclose(domain.quantities['ymomentum'].centroid_values[:4],
3553                        [-1.19201077e-003, -7.23647546e-004,
3554                        -6.39083123e-005, 6.29815168e-005])
3555
3556        os.remove(domain.get_name() + '.sww')
3557
3558    def test_complex_bed(self):
3559        #No friction is tested here
3560
3561        from mesh_factory import rectangular
3562        from Numeric import array
3563
3564        N = 12
3565        points, vertices, boundary = rectangular(N, N/2, len1=1.2,len2=0.6,
3566                                                 origin=(-0.07, 0))
3567
3568
3569        domain = Domain(points, vertices, boundary)
3570        domain.smooth = False
3571        domain.default_order=2
3572
3573
3574        inflow_stage = 0.1
3575        Z = Weir(inflow_stage)
3576        domain.set_quantity('elevation', Z)
3577
3578        Br = Reflective_boundary(domain)
3579        Bd = Dirichlet_boundary([inflow_stage, 0.0, 0.0])
3580        domain.set_boundary({'left': Bd, 'right': Br, 'bottom': Br, 'top': Br})
3581
3582        domain.set_quantity('stage', Constant_height(Z, 0.))
3583
3584        for t in domain.evolve(yieldstep = 0.02, finaltime = 0.2):
3585            pass
3586
3587
3588        #print domain.quantities['stage'].centroid_values
3589
3590        #FIXME: These numbers were from version before 25/10
3591        #assert allclose(domain.quantities['stage'].centroid_values,
3592# [3.95822638e-002,  5.61022588e-002,  4.66437868e-002,  5.73081011e-002,
3593#  4.72394613e-002,  5.74684939e-002,  4.74309483e-002,  5.77458084e-002,
3594#  4.80628177e-002,  5.85656225e-002,  4.90498542e-002,  6.02609831e-002,
3595#  1.18470315e-002,  1.75136443e-002,  1.18035266e-002,  2.15565695e-002,
3596#  1.31620268e-002,  2.14351640e-002,  1.32351076e-002,  2.15450687e-002,
3597#  1.36414028e-002,  2.24274619e-002,  1.51689511e-002,  2.21789655e-002,
3598# -7.54337535e-003, -6.86362021e-004, -7.74146760e-003, -1.83756530e-003,
3599# -8.16773628e-003, -4.49916813e-004, -8.08202599e-003, -3.91118720e-004,
3600# -8.10292716e-003, -3.88584984e-004, -7.35226124e-003,  2.73985295e-004,
3601#  1.86166683e-001,  8.74070369e-002,  1.86166712e-001,  8.74035875e-002,
3602#  6.11666935e-002, -3.76173225e-002, -6.38333276e-002, -3.76147365e-002,
3603#  6.11666725e-002,  8.73846774e-002,  1.86166697e-001,  8.74171550e-002,
3604# -4.83333333e-002,  1.18333333e-001, -4.83333333e-002,  1.18333333e-001,
3605# -4.83333333e-002, -6.66666667e-003, -1.73333333e-001, -1.31666667e-001,
3606# -1.73333333e-001, -6.66666667e-003, -4.83333333e-002,  1.18333333e-001,
3607# -2.48333333e-001, -2.31666667e-001, -2.48333333e-001, -2.31666667e-001,
3608# -2.48333333e-001, -2.31666667e-001, -2.48333333e-001, -2.31666667e-001,
3609# -2.48333333e-001, -2.31666667e-001, -2.48333333e-001, -2.31666667e-001,
3610# -4.65000000e-001, -3.65000000e-001, -4.65000000e-001, -3.65000000e-001,
3611# -4.65000000e-001, -3.65000000e-001, -4.65000000e-001, -3.65000000e-001,
3612# -4.65000000e-001, -3.65000000e-001, -4.65000000e-001, -3.65000000e-001,
3613# -5.98333333e-001, -5.81666667e-001, -5.98333333e-001, -5.81666667e-001,
3614# -5.98333333e-001, -5.81666667e-001, -5.98333333e-001, -5.81666667e-001,
3615# -5.98333333e-001, -5.81666667e-001, -5.98333333e-001, -5.81666667e-001,
3616# -6.48333333e-001, -6.31666667e-001, -6.48333333e-001, -6.31666667e-001,
3617# -6.48333333e-001, -6.31666667e-001, -6.48333333e-001, -6.31666667e-001,
3618# -6.48333333e-001, -6.31666667e-001, -6.48333333e-001, -6.31666667e-001,
3619# -5.31666667e-001, -5.98333333e-001, -5.31666667e-001, -5.98333333e-001,
3620# -5.31666667e-001, -5.98333333e-001, -5.31666667e-001, -5.98333333e-001,
3621# -5.31666667e-001, -5.98333333e-001, -5.31666667e-001, -5.98333333e-001,
3622# -4.98333333e-001, -4.81666667e-001, -4.98333333e-001, -4.81666667e-001,
3623# -4.98333333e-001, -4.81666667e-001, -4.98333333e-001, -4.81666667e-001,
3624# -4.98333333e-001, -4.81666667e-001, -4.98333333e-001, -4.81666667e-001,
3625# -5.48333333e-001, -5.31666667e-001, -5.48333333e-001, -5.31666667e-001,
3626# -5.48333333e-001, -5.31666667e-001, -5.48333333e-001, -5.31666667e-001,
3627# -5.48333333e-001, -5.31666667e-001, -5.48333333e-001, -5.31666667e-001])
3628
3629        os.remove(domain.get_name() + '.sww')
3630
3631    def test_spatio_temporal_boundary_1(self):
3632        """Test that boundary values can be read from file and interpolated
3633        in both time and space.
3634
3635        Verify that the same steady state solution is arrived at and that
3636        time interpolation works.
3637
3638        The full solution history is not exactly the same as
3639        file boundary must read and interpolate from *smoothed* version
3640        as stored in sww.
3641        """
3642        import time
3643
3644        #Create sww file of simple propagation from left to right
3645        #through rectangular domain
3646
3647        from mesh_factory import rectangular
3648
3649        #Create basic mesh
3650        points, vertices, boundary = rectangular(3, 3)
3651
3652        #Create shallow water domain
3653        domain1 = Domain(points, vertices, boundary)
3654
3655        domain1.reduction = mean
3656        domain1.smooth = False  #Exact result
3657
3658        domain1.default_order = 2
3659        domain1.store = True
3660        domain1.set_datadir('.')
3661        domain1.set_name('spatio_temporal_boundary_source' + str(time.time()))
3662
3663        #FIXME: This is extremely important!
3664        #How can we test if they weren't stored?
3665        domain1.quantities_to_be_stored = ['stage', 'xmomentum', 'ymomentum']
3666
3667
3668        #Bed-slope and friction at vertices (and interpolated elsewhere)
3669        domain1.set_quantity('elevation', 0)
3670        domain1.set_quantity('friction', 0)
3671
3672        # Boundary conditions
3673        Br = Reflective_boundary(domain1)
3674        Bd = Dirichlet_boundary([0.3,0,0])
3675        domain1.set_boundary({'left': Bd, 'top': Bd, 'right': Br, 'bottom': Br})
3676        #Initial condition
3677        domain1.set_quantity('stage', 0)
3678        domain1.check_integrity()
3679
3680        finaltime = 5
3681        #Evolution  (full domain - large steps)
3682        for t in domain1.evolve(yieldstep = 0.671, finaltime = finaltime):
3683            pass
3684            #domain1.write_time()
3685
3686        cv1 = domain1.quantities['stage'].centroid_values
3687
3688
3689        #Create a triangle shaped domain (reusing coordinates from domain 1),
3690        #formed from the lower and right hand  boundaries and
3691        #the sw-ne diagonal
3692        #from domain 1. Call it domain2
3693
3694        points = [ [0,0], [1.0/3,0], [1.0/3,1.0/3],
3695                   [2.0/3,0], [2.0/3,1.0/3], [2.0/3,2.0/3],
3696                   [1,0], [1,1.0/3], [1,2.0/3], [1,1]]
3697
3698        vertices = [ [1,2,0], [3,4,1], [2,1,4], [4,5,2],
3699                     [6,7,3], [4,3,7], [7,8,4], [5,4,8], [8,9,5]]
3700
3701        boundary = { (0,1):'bottom', (1,1):'bottom', (4,1): 'bottom',
3702                     (4,2):'right', (6,2):'right', (8,2):'right',
3703                     (0,0):'diagonal', (3,0):'diagonal', (8,0):'diagonal'}
3704
3705        domain2 = Domain(points, vertices, boundary)
3706
3707        domain2.reduction = domain1.reduction
3708        domain2.smooth = False
3709        domain2.default_order = 2
3710
3711        #Bed-slope and friction at vertices (and interpolated elsewhere)
3712        domain2.set_quantity('elevation', 0)
3713        domain2.set_quantity('friction', 0)
3714        domain2.set_quantity('stage', 0)
3715
3716        # Boundary conditions
3717        Br = Reflective_boundary(domain2)
3718        #Bf = Spatio_temporal_boundary(domain1.get_name() + '.' +\
3719        #                              domain1.format, domain2)
3720        Bf = Field_boundary(domain1.get_name() + '.' +\
3721                            domain1.format, domain2)       
3722        domain2.set_boundary({'right':Br, 'bottom':Br, 'diagonal':Bf})
3723        domain2.check_integrity()
3724
3725
3726
3727        #Evolution (small steps)
3728        for t in domain2.evolve(yieldstep = 0.0711, finaltime = finaltime):
3729            pass
3730
3731
3732        #Use output from domain1 as spatio-temporal boundary for domain2
3733        #and verify that results at right hand side are close.
3734
3735        cv2 = domain2.quantities['stage'].centroid_values
3736
3737        #print take(cv1, (12,14,16))  #Right
3738        #print take(cv2, (4,6,8))
3739        #print take(cv1, (0,6,12))  #Bottom
3740        #print take(cv2, (0,1,4))
3741        #print take(cv1, (0,8,16))  #Diag
3742        #print take(cv2, (0,3,8))
3743
3744        assert allclose( take(cv1, (0,8,16)), take(cv2, (0,3,8))) #Diag
3745        assert allclose( take(cv1, (0,6,12)), take(cv2, (0,1,4))) #Bottom
3746        assert allclose( take(cv1, (12,14,16)), take(cv2, (4,6,8))) #RHS
3747
3748        #Cleanup
3749        os.remove(domain1.get_name() + '.' + domain1.format)
3750        os.remove(domain2.get_name() + '.' + domain2.format)       
3751
3752
3753
3754    def test_spatio_temporal_boundary_2(self):
3755        """Test that boundary values can be read from file and interpolated
3756        in both time and space.
3757        This is a more basic test, verifying that boundary object
3758        produces the expected results
3759
3760
3761        """
3762        import time
3763
3764        #Create sww file of simple propagation from left to right
3765        #through rectangular domain
3766
3767        from mesh_factory import rectangular
3768
3769        #Create basic mesh
3770        points, vertices, boundary = rectangular(3, 3)
3771
3772        #Create shallow water domain
3773        domain1 = Domain(points, vertices, boundary)
3774
3775        domain1.reduction = mean
3776        domain1.smooth = True #To mimic MOST output
3777
3778        domain1.default_order = 2
3779        domain1.store = True
3780        domain1.set_datadir('.')
3781        domain1.set_name('spatio_temporal_boundary_source' + str(time.time()))
3782
3783        #FIXME: This is extremely important!
3784        #How can we test if they weren't stored?
3785        domain1.quantities_to_be_stored = ['stage', 'xmomentum', 'ymomentum']
3786
3787
3788        #Bed-slope and friction at vertices (and interpolated elsewhere)
3789        domain1.set_quantity('elevation', 0)
3790        domain1.set_quantity('friction', 0)
3791
3792        # Boundary conditions
3793        Br = Reflective_boundary(domain1)
3794        Bd = Dirichlet_boundary([0.3,0,0])
3795        domain1.set_boundary({'left': Bd, 'top': Bd, 'right': Br, 'bottom': Br})
3796        #Initial condition
3797        domain1.set_quantity('stage', 0)
3798        domain1.check_integrity()
3799
3800        finaltime = 5
3801        #Evolution  (full domain - large steps)
3802        for t in domain1.evolve(yieldstep = 1, finaltime = finaltime):
3803            pass
3804            #domain1.write_time()
3805
3806
3807        #Create an triangle shaped domain (coinciding with some
3808        #coordinates from domain 1),
3809        #formed from the lower and right hand  boundaries and
3810        #the sw-ne diagonal
3811        #from domain 1. Call it domain2
3812
3813        points = [ [0,0], [1.0/3,0], [1.0/3,1.0/3],
3814                   [2.0/3,0], [2.0/3,1.0/3], [2.0/3,2.0/3],
3815                   [1,0], [1,1.0/3], [1,2.0/3], [1,1]]
3816
3817        vertices = [ [1,2,0],
3818                     [3,4,1], [2,1,4], [4,5,2],
3819                     [6,7,3], [4,3,7], [7,8,4], [5,4,8], [8,9,5]]
3820
3821        boundary = { (0,1):'bottom', (1,1):'bottom', (4,1): 'bottom',
3822                     (4,2):'right', (6,2):'right', (8,2):'right',
3823                     (0,0):'diagonal', (3,0):'diagonal', (8,0):'diagonal'}
3824
3825        domain2 = Domain(points, vertices, boundary)
3826
3827        domain2.reduction = domain1.reduction
3828        domain2.smooth = False
3829        domain2.default_order = 2
3830
3831        #Bed-slope and friction at vertices (and interpolated elsewhere)
3832        domain2.set_quantity('elevation', 0)
3833        domain2.set_quantity('friction', 0)
3834        domain2.set_quantity('stage', 0)
3835
3836
3837        #Read results for specific timesteps t=1 and t=2
3838        from Scientific.IO.NetCDF import NetCDFFile
3839        fid = NetCDFFile(domain1.get_name() + '.' + domain1.format)
3840
3841        x = fid.variables['x'][:]
3842        y = fid.variables['y'][:]
3843        s1 = fid.variables['stage'][1,:]
3844        s2 = fid.variables['stage'][2,:]
3845        fid.close()
3846
3847        from Numeric import take, reshape, concatenate
3848        shp = (len(x), 1)
3849        points = concatenate( (reshape(x, shp), reshape(y, shp)), axis=1)
3850        #The diagonal points of domain 1 are 0, 5, 10, 15
3851
3852        #print points[0], points[5], points[10], points[15]
3853        assert allclose( take(points, [0,5,10,15]),
3854                         [[0,0], [1.0/3, 1.0/3], [2.0/3, 2.0/3], [1,1]])
3855
3856
3857        # Boundary conditions
3858        Br = Reflective_boundary(domain2)
3859        #Bf = Spatio_temporal_boundary(domain1.get_name() + '.' + domain1.format,
3860        #                              domain2)
3861        Bf = Field_boundary(domain1.get_name() + '.' + domain1.format,
3862                            domain2)       
3863        domain2.set_boundary({'right':Br, 'bottom':Br, 'diagonal':Bf})
3864        domain2.check_integrity()
3865
3866        #Test that interpolation points are the mid points of the all boundary
3867        #segments
3868
3869        boundary_midpoints = [[1.0/6, 0], [1.0/2, 0], [5.0/6,0],
3870                              [1.0, 1.0/6], [1.0, 1.0/2], [1.0, 5.0/6],
3871                              [1.0/6, 1.0/6], [0.5, 0.5], [5.0/6, 5.0/6]]
3872
3873        boundary_midpoints.sort()
3874        R = Bf.F.interpolation_points.tolist()
3875        R.sort()
3876        assert allclose(boundary_midpoints, R)
3877
3878        #Check spatially interpolated output at time == 1
3879        domain2.time = 1
3880
3881        #First diagonal midpoint
3882        R0 = Bf.evaluate(0,0)
3883        assert allclose(R0[0], (s1[0] + s1[5])/2)
3884
3885        #Second diagonal midpoint
3886        R0 = Bf.evaluate(3,0)
3887        assert allclose(R0[0], (s1[5] + s1[10])/2)
3888
3889        #First diagonal midpoint
3890        R0 = Bf.evaluate(8,0)
3891        assert allclose(R0[0], (s1[10] + s1[15])/2)
3892
3893        #Check spatially interpolated output at time == 2
3894        domain2.time = 2
3895
3896        #First diagonal midpoint
3897        R0 = Bf.evaluate(0,0)
3898        assert allclose(R0[0], (s2[0] + s2[5])/2)
3899
3900        #Second diagonal midpoint
3901        R0 = Bf.evaluate(3,0)
3902        assert allclose(R0[0], (s2[5] + s2[10])/2)
3903
3904        #First diagonal midpoint
3905        R0 = Bf.evaluate(8,0)
3906        assert allclose(R0[0], (s2[10] + s2[15])/2)
3907
3908
3909        #Now check temporal interpolation
3910
3911        domain2.time = 1 + 2.0/3
3912
3913        #First diagonal midpoint
3914        R0 = Bf.evaluate(0,0)
3915        assert allclose(R0[0], ((s1[0] + s1[5])/2 + 2.0*(s2[0] + s2[5])/2)/3)
3916
3917        #Second diagonal midpoint
3918        R0 = Bf.evaluate(3,0)
3919        assert allclose(R0[0], ((s1[5] + s1[10])/2 + 2.0*(s2[5] + s2[10])/2)/3)
3920
3921        #First diagonal midpoint
3922        R0 = Bf.evaluate(8,0)
3923        assert allclose(R0[0], ((s1[10] + s1[15])/2 + 2.0*(s2[10] + s2[15])/2)/3)
3924
3925
3926
3927        #Cleanup
3928        os.remove(domain1.get_name() + '.' + domain1.format)
3929
3930
3931    def test_spatio_temporal_boundary_3(self):
3932        """Test that boundary values can be read from file and interpolated
3933        in both time and space.
3934        This is a more basic test, verifying that boundary object
3935        produces the expected results
3936
3937        This tests adjusting using mean_stage
3938
3939        """
3940
3941        import time
3942
3943        mean_stage = 5.2 # Adjust stage by this amount in boundary
3944
3945        #Create sww file of simple propagation from left to right
3946        #through rectangular domain
3947
3948        from mesh_factory import rectangular
3949
3950        #Create basic mesh
3951        points, vertices, boundary = rectangular(3, 3)
3952
3953        #Create shallow water domain
3954        domain1 = Domain(points, vertices, boundary)
3955
3956        domain1.reduction = mean
3957        domain1.smooth = True #To mimic MOST output
3958
3959        domain1.default_order = 2
3960        domain1.store = True
3961        domain1.set_datadir('.')
3962        domain1.set_name('spatio_temporal_boundary_source' + str(time.time()))
3963
3964        #FIXME: This is extremely important!
3965        #How can we test if they weren't stored?
3966        domain1.quantities_to_be_stored = ['stage', 'xmomentum', 'ymomentum']
3967
3968
3969        #Bed-slope and friction at vertices (and interpolated elsewhere)
3970        domain1.set_quantity('elevation', 0)
3971        domain1.set_quantity('friction', 0)
3972
3973        # Boundary conditions
3974        Br = Reflective_boundary(domain1)
3975        Bd = Dirichlet_boundary([0.3,0,0])
3976        domain1.set_boundary({'left': Bd, 'top': Bd, 'right': Br, 'bottom': Br})
3977        #Initial condition
3978        domain1.set_quantity('stage', 0)
3979        domain1.check_integrity()
3980
3981        finaltime = 5
3982        #Evolution  (full domain - large steps)
3983        for t in domain1.evolve(yieldstep = 1, finaltime = finaltime):
3984            pass
3985            #domain1.write_time()
3986
3987
3988        #Create an triangle shaped domain (coinciding with some
3989        #coordinates from domain 1),
3990        #formed from the lower and right hand  boundaries and
3991        #the sw-ne diagonal
3992        #from domain 1. Call it domain2
3993
3994        points = [ [0,0], [1.0/3,0], [1.0/3,1.0/3],
3995                   [2.0/3,0], [2.0/3,1.0/3], [2.0/3,2.0/3],
3996                   [1,0], [1,1.0/3], [1,2.0/3], [1,1]]
3997
3998        vertices = [ [1,2,0],
3999                     [3,4,1], [2,1,4], [4,5,2],
4000                     [6,7,3], [4,3,7], [7,8,4], [5,4,8], [8,9,5]]
4001
4002        boundary = { (0,1):'bottom', (1,1):'bottom', (4,1): 'bottom',
4003                     (4,2):'right', (6,2):'right', (8,2):'right',
4004                     (0,0):'diagonal', (3,0):'diagonal', (8,0):'diagonal'}
4005
4006        domain2 = Domain(points, vertices, boundary)
4007
4008        domain2.reduction = domain1.reduction
4009        domain2.smooth = False
4010        domain2.default_order = 2
4011
4012        #Bed-slope and friction at vertices (and interpolated elsewhere)
4013        domain2.set_quantity('elevation', 0)
4014        domain2.set_quantity('friction', 0)
4015        domain2.set_quantity('stage', 0)
4016
4017
4018        #Read results for specific timesteps t=1 and t=2
4019        from Scientific.IO.NetCDF import NetCDFFile
4020        fid = NetCDFFile(domain1.get_name() + '.' + domain1.format)
4021
4022        x = fid.variables['x'][:]
4023        y = fid.variables['y'][:]
4024        s1 = fid.variables['stage'][1,:]
4025        s2 = fid.variables['stage'][2,:]
4026        fid.close()
4027
4028        from Numeric import take, reshape, concatenate
4029        shp = (len(x), 1)
4030        points = concatenate( (reshape(x, shp), reshape(y, shp)), axis=1)
4031        #The diagonal points of domain 1 are 0, 5, 10, 15
4032
4033        #print points[0], points[5], points[10], points[15]
4034        assert allclose( take(points, [0,5,10,15]),
4035                         [[0,0], [1.0/3, 1.0/3], [2.0/3, 2.0/3], [1,1]])
4036
4037
4038        # Boundary conditions
4039        Br = Reflective_boundary(domain2)
4040        #Bf = Spatio_temporal_boundary(domain1.get_name() + '.' + domain1.format,
4041        #                              domain2)
4042        Bf = Field_boundary(domain1.get_name() + '.' + domain1.format,
4043                            domain2, mean_stage=mean_stage)
4044       
4045        domain2.set_boundary({'right':Br, 'bottom':Br, 'diagonal':Bf})
4046        domain2.check_integrity()
4047
4048        #Test that interpolation points are the mid points of the all boundary
4049        #segments
4050
4051        boundary_midpoints = [[1.0/6, 0], [1.0/2, 0], [5.0/6,0],
4052                              [1.0, 1.0/6], [1.0, 1.0/2], [1.0, 5.0/6],
4053                              [1.0/6, 1.0/6], [0.5, 0.5], [5.0/6, 5.0/6]]
4054
4055        boundary_midpoints.sort()
4056        R = Bf.F.interpolation_points.tolist()
4057        R.sort()
4058        assert allclose(boundary_midpoints, R)
4059
4060        #Check spatially interpolated output at time == 1
4061        domain2.time = 1
4062
4063        #First diagonal midpoint
4064        R0 = Bf.evaluate(0,0)
4065        assert allclose(R0[0], (s1[0] + s1[5])/2 + mean_stage)
4066
4067        #Second diagonal midpoint
4068        R0 = Bf.evaluate(3,0)
4069        assert allclose(R0[0], (s1[5] + s1[10])/2 + mean_stage)
4070
4071        #First diagonal midpoint
4072        R0 = Bf.evaluate(8,0)
4073        assert allclose(R0[0], (s1[10] + s1[15])/2 + mean_stage)
4074
4075        #Check spatially interpolated output at time == 2
4076        domain2.time = 2
4077
4078        #First diagonal midpoint
4079        R0 = Bf.evaluate(0,0)
4080        assert allclose(R0[0], (s2[0] + s2[5])/2 + mean_stage)
4081
4082        #Second diagonal midpoint
4083        R0 = Bf.evaluate(3,0)
4084        assert allclose(R0[0], (s2[5] + s2[10])/2 + mean_stage)
4085
4086        #First diagonal midpoint
4087        R0 = Bf.evaluate(8,0)
4088        assert allclose(R0[0], (s2[10] + s2[15])/2 + mean_stage)
4089
4090
4091        #Now check temporal interpolation
4092
4093        domain2.time = 1 + 2.0/3
4094
4095        #First diagonal midpoint
4096        R0 = Bf.evaluate(0,0)
4097        assert allclose(R0[0], ((s1[0] + s1[5])/2 + 2.0*(s2[0] + s2[5])/2)/3 + mean_stage)
4098
4099        #Second diagonal midpoint
4100        R0 = Bf.evaluate(3,0)
4101        assert allclose(R0[0], ((s1[5] + s1[10])/2 + 2.0*(s2[5] + s2[10])/2)/3 + mean_stage)
4102
4103        #First diagonal midpoint
4104        R0 = Bf.evaluate(8,0)
4105        assert allclose(R0[0], ((s1[10] + s1[15])/2 + 2.0*(s2[10] + s2[15])/2)/3 + mean_stage)
4106
4107
4108        #Cleanup
4109        os.remove(domain1.get_name() + '.' + domain1.format)
4110
4111
4112    def test_pmesh2Domain(self):
4113         import os
4114         import tempfile
4115
4116         fileName = tempfile.mktemp(".tsh")
4117         file = open(fileName,"w")
4118         file.write("4 3 # <vertex #> <x> <y> [attributes]\n \
41190 0.0 0.0 0.0 0.0 0.01 \n \
41201 1.0 0.0 10.0 10.0 0.02  \n \
41212 0.0 1.0 0.0 10.0 0.03  \n \
41223 0.5 0.25 8.0 12.0 0.04  \n \
4123# Vert att title  \n \
4124elevation  \n \
4125stage  \n \
4126friction  \n \
41272 # <triangle #> [<vertex #>] [<neigbouring triangle #>]  \n\
41280 0 3 2 -1  -1  1 dsg\n\
41291 0 1 3 -1  0 -1   ole nielsen\n\
41304 # <segment #> <vertex #>  <vertex #> [boundary tag] \n\
41310 1 0 2 \n\
41321 0 2 3 \n\
41332 2 3 \n\
41343 3 1 1 \n\
41353 0 # <x> <y> [attributes] ...Mesh Vertices... \n \
41360 216.0 -86.0 \n \
41371 160.0 -167.0 \n \
41382 114.0 -91.0 \n \
41393 # <vertex #>  <vertex #> [boundary tag] ...Mesh Segments... \n \
41400 0 1 0 \n \
41411 1 2 0 \n \
41422 2 0 0 \n \
41430 # <x> <y> ...Mesh Holes... \n \
41440 # <x> <y> <attribute>...Mesh Regions... \n \
41450 # <x> <y> <attribute>...Mesh Regions, area... \n\
4146#Geo reference \n \
414756 \n \
4148140 \n \
4149120 \n")
4150         file.close()
4151
4152         tags = {}
4153         b1 =  Dirichlet_boundary(conserved_quantities = array([0.0]))
4154         b2 =  Dirichlet_boundary(conserved_quantities = array([1.0]))
4155         b3 =  Dirichlet_boundary(conserved_quantities = array([2.0]))
4156         tags["1"] = b1
4157         tags["2"] = b2
4158         tags["3"] = b3
4159
4160         #from anuga.abstract_2d_finite_volumes.pmesh2domain import pmesh_to_domain_instance
4161         #domain = pmesh_to_domain_instance(fileName, Domain)
4162
4163         domain = Domain(mesh_filename=fileName)
4164                         #verbose=True, use_cache=True)
4165         
4166         #print "domain.tagged_elements", domain.tagged_elements
4167         ## check the quantities
4168         #print domain.quantities['elevation'].vertex_values
4169         answer = [[0., 8., 0.],
4170                   [0., 10., 8.]]
4171         assert allclose(domain.quantities['elevation'].vertex_values,
4172                        answer)
4173
4174         #print domain.quantities['stage'].vertex_values
4175         answer = [[0., 12., 10.],
4176                   [0., 10., 12.]]
4177         assert allclose(domain.quantities['stage'].vertex_values,
4178                        answer)
4179
4180         #print domain.quantities['friction'].vertex_values
4181         answer = [[0.01, 0.04, 0.03],
4182                   [0.01, 0.02, 0.04]]
4183         assert allclose(domain.quantities['friction'].vertex_values,
4184                        answer)
4185
4186         #print domain.quantities['friction'].vertex_values
4187         assert allclose(domain.tagged_elements['dsg'][0],0)
4188         assert allclose(domain.tagged_elements['ole nielsen'][0],1)
4189
4190         self.failUnless( domain.boundary[(1, 0)]  == '1',
4191                          "test_tags_to_boundaries  failed. Single boundary wasn't added.")
4192         self.failUnless( domain.boundary[(1, 2)]  == '2',
4193                          "test_tags_to_boundaries  failed. Single boundary wasn't added.")
4194         self.failUnless( domain.boundary[(0, 1)]  == '3',
4195                          "test_tags_to_boundaries  failed. Single boundary wasn't added.")
4196         self.failUnless( domain.boundary[(0, 0)]  == 'exterior',
4197                          "test_tags_to_boundaries  failed. Single boundary wasn't added.")
4198         #print "domain.boundary",domain.boundary
4199         self.failUnless( len(domain.boundary)  == 4,
4200                          "test_pmesh2Domain Too many boundaries")
4201         #FIXME change to use get_xllcorner
4202         #print "d.geo_reference.xllcorner",domain.geo_reference.xllcorner
4203         self.failUnless(domain.geo_reference.xllcorner  == 140.0,
4204                          "bad geo_referece")
4205
4206
4207         #************
4208
4209   
4210         domain = Domain(fileName)
4211         
4212         #print "domain.tagged_elements", domain.tagged_elements
4213         ## check the quantities
4214         #print domain.quantities['elevation'].vertex_values
4215         answer = [[0., 8., 0.],
4216                   [0., 10., 8.]]
4217         assert allclose(domain.quantities['elevation'].vertex_values,
4218                        answer)
4219
4220         #print domain.quantities['stage'].vertex_values
4221         answer = [[0., 12., 10.],
4222                   [0., 10., 12.]]
4223         assert allclose(domain.quantities['stage'].vertex_values,
4224                        answer)
4225
4226         #print domain.quantities['friction'].vertex_values
4227         answer = [[0.01, 0.04, 0.03],
4228                   [0.01, 0.02, 0.04]]
4229         assert allclose(domain.quantities['friction'].vertex_values,
4230                        answer)
4231
4232         #print domain.quantities['friction'].vertex_values
4233         assert allclose(domain.tagged_elements['dsg'][0],0)
4234         assert allclose(domain.tagged_elements['ole nielsen'][0],1)
4235
4236         self.failUnless( domain.boundary[(1, 0)]  == '1',
4237                          "test_tags_to_boundaries  failed. Single boundary wasn't added.")
4238         self.failUnless( domain.boundary[(1, 2)]  == '2',
4239                          "test_tags_to_boundaries  failed. Single boundary wasn't added.")
4240         self.failUnless( domain.boundary[(0, 1)]  == '3',
4241                          "test_tags_to_boundaries  failed. Single boundary wasn't added.")
4242         self.failUnless( domain.boundary[(0, 0)]  == 'exterior',
4243                          "test_tags_to_boundaries  failed. Single boundary wasn't added.")
4244         #print "domain.boundary",domain.boundary
4245         self.failUnless( len(domain.boundary)  == 4,
4246                          "test_pmesh2Domain Too many boundaries")
4247         #FIXME change to use get_xllcorner
4248         #print "d.geo_reference.xllcorner",domain.geo_reference.xllcorner
4249         self.failUnless(domain.geo_reference.xllcorner  == 140.0,
4250                          "bad geo_referece")
4251         #************
4252         os.remove(fileName)
4253
4254        #-------------------------------------------------------------
4255       
4256if __name__ == "__main__":
4257    suite = unittest.makeSuite(Test_Shallow_Water,'test')
4258    #suite = unittest.makeSuite(Test_Shallow_Water,'test_get_maximum_inundation_3')
4259    runner = unittest.TextTestRunner(verbosity=1)   
4260    runner.run(suite)
Note: See TracBrowser for help on using the repository browser.