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

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

New test of runup timeseries - are they consistent across computer architectures?

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