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

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

Implemented optimisation excluding dry cells from flux calculation.
All tests pass and the script run_okushiri_profile.py reports an improvement
in compute_fluxes from 11.25s to 8.58s or 24% faster.
The overall computation was about 40s, so this optimisation improved the
total running time for the problem in question by 7%.

This partially addresses ticket:135

File size: 177.4 KB
Line 
1#!/usr/bin/env python
2
3import unittest, os
4from math import sqrt, pi
5
6from anuga.config import g, epsilon
7from Numeric import allclose, alltrue, array, zeros, ones, Float, take
8from anuga.utilities.numerical_tools import mean
9from anuga.utilities.polygon import is_inside_polygon
10
11from shallow_water_domain import *
12from shallow_water_domain import flux_function_central as flux_function
13
14#Variable windfield implemented using functions
15def speed(t,x,y):
16    """Large speeds halfway between center and edges
17    Low speeds at center and edges
18    """
19
20    from math import sqrt, exp, cos, pi
21
22    x = array(x)
23    y = array(y)
24
25    N = len(x)
26    s = 0*#New array
27
28    for k in range(N):
29
30        r = sqrt(x[k]**2 + y[k]**2)
31
32        factor = exp( -(r-0.15)**2 )
33
34        s[k] = 4000 * factor * (cos(t*2*pi/150) + 2)
35
36    return s
37
38
39def scalar_func(t,x,y):
40    """Function that returns a scalar.
41    Used to test error message when Numeric array is expected
42    """
43
44    return 17.7
45
46
47def angle(t,x,y):
48    """Rotating field
49    """
50    from math import sqrt, atan, pi
51
52    x = array(x)
53    y = array(y)
54
55    N = len(x)
56    a = 0*#New array
57
58    for k in range(N):
59        r = sqrt(x[k]**2 + y[k]**2)
60
61        angle = atan(y[k]/x[k])
62
63        if x[k] < 0:
64            angle+=pi #atan in ]-pi/2; pi/2[
65
66        #Take normal direction
67        angle -= pi/2
68
69        #Ensure positive radians
70        if angle < 0:
71            angle += 2*pi
72
73        a[k] = angle/pi*180
74
75    return a
76
77
78class Test_Shallow_Water(unittest.TestCase):
79    def setUp(self):
80        pass
81
82    def tearDown(self):
83        pass
84
85    def test_rotate(self):
86        normal = array([0.0,-1.0])
87
88        q = array([1.0,2.0,3.0])
89
90        r = rotate(q, normal, direction = 1)
91        assert r[0] == 1
92        assert r[1] == -3
93        assert r[2] == 2
94
95        w = rotate(r, normal, direction = -1)
96        assert allclose(w, q)
97
98        #Check error check
99        try:
100            rotate(r, array([1,1,1]) )
101        except:
102            pass
103        else:
104            raise 'Should have raised an exception'
105
106
107    #FIXME (Ole): Individual flux tests do NOT test C implementation directly.   
108    def test_flux_zero_case(self):
109        ql = zeros( 3, Float )
110        qr = zeros( 3, Float )
111        normal = zeros( 2, Float )
112        zl = zr = 0.
113        flux, max_speed = flux_function(normal, ql, qr, zl, zr)
114
115        assert allclose(flux, [0,0,0])
116        assert max_speed == 0.
117
118    def test_flux_constants(self):
119        w = 2.0
120
121        normal = array([1.,0])
122        ql = array([w, 0, 0])
123        qr = array([w, 0, 0])
124        zl = zr = 0.
125        h = w - (zl+zr)/2
126
127        flux, max_speed = flux_function(normal, ql, qr, zl, zr)
128
129        assert allclose(flux, [0., 0.5*g*h**2, 0.])
130        assert max_speed == sqrt(g*h)
131
132    #def test_flux_slope(self):
133    #    #FIXME: TODO
134    #    w = 2.0
135    #
136    #    normal = array([1.,0])
137    #    ql = array([w, 0, 0])
138    #    qr = array([w, 0, 0])
139    #    zl = zr = 0.
140    #    h = w - (zl+zr)/2
141    #
142    #    flux, max_speed = flux_function(normal, ql, qr, zl, zr)
143    #
144    #    assert allclose(flux, [0., 0.5*g*h**2, 0.])
145    #    assert max_speed == sqrt(g*h)
146
147
148    def test_flux1(self):
149        #Use data from previous version of abstract_2d_finite_volumes
150        normal = array([1.,0])
151        ql = array([-0.2, 2, 3])
152        qr = array([-0.2, 2, 3])
153        zl = zr = -0.5
154        flux, max_speed = flux_function(normal, ql, qr, zl, zr)
155
156        assert allclose(flux, [2.,13.77433333, 20.])
157        assert allclose(max_speed, 8.38130948661)
158
159
160    def test_flux2(self):
161        #Use data from previous version of abstract_2d_finite_volumes
162        normal = array([0., -1.])
163        ql = array([-0.075, 2, 3])
164        qr = array([-0.075, 2, 3])
165        zl = zr = -0.375
166        flux, max_speed = flux_function(normal, ql, qr, zl, zr)
167
168        assert allclose(flux, [-3.,-20.0, -30.441])
169        assert allclose(max_speed, 11.7146428199)
170
171    def test_flux3(self):
172        #Use data from previous version of abstract_2d_finite_volumes
173        normal = array([-sqrt(2)/2, sqrt(2)/2])
174        ql = array([-0.075, 2, 3])
175        qr = array([-0.075, 2, 3])
176        zl = zr = -0.375
177        flux, max_speed = flux_function(normal, ql, qr, zl, zr)
178
179        assert allclose(flux, [sqrt(2)/2, 4.40221112, 7.3829019])
180        assert allclose(max_speed, 4.0716654239)
181
182    def test_flux4(self):
183        #Use data from previous version of abstract_2d_finite_volumes
184        normal = array([-sqrt(2)/2, sqrt(2)/2])
185        ql = array([-0.34319278, 0.10254161, 0.07273855])
186        qr = array([-0.30683287, 0.1071986, 0.05930515])
187        zl = zr = -0.375
188        flux, max_speed = flux_function(normal, ql, qr, zl, zr)
189
190        assert allclose(flux, [-0.04072676, -0.07096636, -0.01604364])
191        assert allclose(max_speed, 1.31414103233)
192
193    def test_flux_computation(self):   
194        """test_flux_computation - test flux calculation (actual C implementation)
195        This one tests the constant case where only the pressure term contributes to each edge and cancels out
196        once the total flux has been summed up.
197        """
198               
199        a = [0.0, 0.0]
200        b = [0.0, 2.0]
201        c = [2.0,0.0]
202        d = [0.0, 4.0]
203        e = [2.0, 2.0]
204        f = [4.0,0.0]
205
206        points = [a, b, c, d, e, f]
207        #bac, bce, ecf, dbe, daf, dae
208        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
209
210        domain = Domain(points, vertices)
211        domain.check_integrity()
212
213        # The constant case             
214        domain.set_quantity('elevation', -1)
215        domain.set_quantity('stage', 1) 
216       
217        domain.compute_fluxes()
218        assert allclose(domain.get_quantity('stage').explicit_update[1], 0) # Central triangle
219       
220
221        # The more general case                 
222        def surface(x,y):
223            return -x/2                   
224       
225        domain.set_quantity('elevation', -10)
226        domain.set_quantity('stage', surface)   
227        domain.set_quantity('xmomentum', 1)             
228       
229        domain.compute_fluxes()
230       
231        #print domain.get_quantity('stage').explicit_update
232        # FIXME (Ole): TODO the general case
233        #assert allclose(domain.get_quantity('stage').explicit_update[1], ........??)
234               
235       
236               
237    def test_sw_domain_simple(self):
238        a = [0.0, 0.0]
239        b = [0.0, 2.0]
240        c = [2.0,0.0]
241        d = [0.0, 4.0]
242        e = [2.0, 2.0]
243        f = [4.0,0.0]
244
245        points = [a, b, c, d, e, f]
246        #bac, bce, ecf, dbe, daf, dae
247        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
248
249
250        #from anuga.abstract_2d_finite_volumes.domain import Domain as Generic_domain
251        #msg = 'The class %s is not a subclass of the generic domain class %s'\
252        #      %(DomainClass, Domain)
253        #assert issubclass(DomainClass, Domain), msg
254
255        domain = Domain(points, vertices)
256        domain.check_integrity()
257
258        for name in ['stage', 'xmomentum', 'ymomentum',
259                     'elevation', 'friction']:
260            assert domain.quantities.has_key(name)
261
262
263        assert domain.get_conserved_quantities(0, edge=1) == 0.
264
265
266    def test_boundary_conditions(self):
267
268        a = [0.0, 0.0]
269        b = [0.0, 2.0]
270        c = [2.0,0.0]
271        d = [0.0, 4.0]
272        e = [2.0, 2.0]
273        f = [4.0,0.0]
274
275        points = [a, b, c, d, e, f]
276        #bac, bce, ecf, dbe
277        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4] ]
278        boundary = { (0, 0): 'Third',
279                     (0, 2): 'First',
280                     (2, 0): 'Second',
281                     (2, 1): 'Second',
282                     (3, 1): 'Second',
283                     (3, 2): 'Third'}
284
285
286        domain = Domain(points, vertices, boundary)
287        domain.check_integrity()
288
289
290        domain.set_quantity('stage', [[1,2,3], [5,5,5],
291                                      [0,0,9], [-6, 3, 3]])
292
293        domain.set_quantity('xmomentum', [[1,1,1], [2,2,2],
294                                          [3,3,3], [4, 4, 4]])
295
296        domain.set_quantity('ymomentum', [[10,10,10], [20,20,20],
297                                          [30,30,30], [40, 40, 40]])
298
299
300        D = Dirichlet_boundary([5,2,1])
301        T = Transmissive_boundary(domain)
302        R = Reflective_boundary(domain)
303        domain.set_boundary( {'First': D, 'Second': T, 'Third': R})
304
305        domain.update_boundary()
306
307        #Stage
308        assert domain.quantities['stage'].boundary_values[0] == 2.5
309        assert domain.quantities['stage'].boundary_values[0] ==\
310               domain.get_conserved_quantities(0, edge=0)[0] #Reflective (2.5)
311        assert domain.quantities['stage'].boundary_values[1] == 5. #Dirichlet
312        assert domain.quantities['stage'].boundary_values[2] ==\
313               domain.get_conserved_quantities(2, edge=0)[0] #Transmissive (4.5)
314        assert domain.quantities['stage'].boundary_values[3] ==\
315               domain.get_conserved_quantities(2, edge=1)[0] #Transmissive (4.5)
316        assert domain.quantities['stage'].boundary_values[4] ==\
317               domain.get_conserved_quantities(3, edge=1)[0] #Transmissive (-1.5)
318        assert domain.quantities['stage'].boundary_values[5] ==\
319               domain.get_conserved_quantities(3, edge=2)[0] #Reflective (-1.5)
320
321        #Xmomentum
322        assert domain.quantities['xmomentum'].boundary_values[0] == 1.0 #Reflective
323        assert domain.quantities['xmomentum'].boundary_values[1] == 2. #Dirichlet
324        assert domain.quantities['xmomentum'].boundary_values[2] ==\
325               domain.get_conserved_quantities(2, edge=0)[1] #Transmissive
326        assert domain.quantities['xmomentum'].boundary_values[3] ==\
327               domain.get_conserved_quantities(2, edge=1)[1] #Transmissive
328        assert domain.quantities['xmomentum'].boundary_values[4] ==\
329               domain.get_conserved_quantities(3, edge=1)[1] #Transmissive
330        assert domain.quantities['xmomentum'].boundary_values[5] == -4.0  #Reflective
331
332        #Ymomentum
333        assert domain.quantities['ymomentum'].boundary_values[0] == -10.0 #Reflective
334        assert domain.quantities['ymomentum'].boundary_values[1] == 1.  #Dirichlet
335        assert domain.quantities['ymomentum'].boundary_values[2] == 30. #Transmissive
336        assert domain.quantities['ymomentum'].boundary_values[3] == 30. #Transmissive
337        assert domain.quantities['ymomentum'].boundary_values[4] == 40. #Transmissive
338        assert domain.quantities['ymomentum'].boundary_values[5] == 40. #Reflective
339
340
341    def test_boundary_conditionsII(self):
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        boundary = { (0, 0): 'Third',
354                     (0, 2): 'First',
355                     (2, 0): 'Second',
356                     (2, 1): 'Second',
357                     (3, 1): 'Second',
358                     (3, 2): 'Third',
359                     (0, 1): 'Internal'}
360
361
362        domain = Domain(points, vertices, boundary)
363        domain.check_integrity()
364
365
366        domain.set_quantity('stage', [[1,2,3], [5,5,5],
367                                      [0,0,9], [-6, 3, 3]])
368
369        domain.set_quantity('xmomentum', [[1,1,1], [2,2,2],
370                                          [3,3,3], [4, 4, 4]])
371
372        domain.set_quantity('ymomentum', [[10,10,10], [20,20,20],
373                                          [30,30,30], [40, 40, 40]])
374
375
376        D = Dirichlet_boundary([5,2,1])
377        T = Transmissive_boundary(domain)
378        R = Reflective_boundary(domain)
379        domain.set_boundary( {'First': D, 'Second': T,
380                              'Third': R, 'Internal': None})
381
382        domain.update_boundary()
383        domain.check_integrity()
384
385
386    def test_compute_fluxes0(self):
387        #Do a full triangle and check that fluxes cancel out for
388        #the constant stage case
389
390        a = [0.0, 0.0]
391        b = [0.0, 2.0]
392        c = [2.0,0.0]
393        d = [0.0, 4.0]
394        e = [2.0, 2.0]
395        f = [4.0,0.0]
396
397        points = [a, b, c, d, e, f]
398        #bac, bce, ecf, dbe
399        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
400
401        domain = Domain(points, vertices)
402        domain.set_quantity('stage', [[2,2,2], [2,2,2],
403                                      [2,2,2], [2,2,2]])
404        domain.check_integrity()
405
406        assert allclose(domain.neighbours, [[-1,1,-1], [2,3,0], [-1,-1,1],[1,-1,-1]])
407        assert allclose(domain.neighbour_edges, [[-1,2,-1], [2,0,1], [-1,-1,0],[1,-1,-1]])
408
409        zl=zr=0. #Assume flat bed
410
411        #Flux across right edge of volume 1
412        normal = domain.get_normal(1,0)
413        ql = domain.get_conserved_quantities(vol_id=1, edge=0)
414        qr = domain.get_conserved_quantities(vol_id=2, edge=2)
415        flux0, max_speed = flux_function(normal, ql, qr, zl, zr)
416
417        #Check that flux seen from other triangles is inverse
418        tmp = qr; qr=ql; ql=tmp
419        normal = domain.get_normal(2,2)
420        flux, max_speed = flux_function(normal, ql, qr, zl, zr)
421        assert allclose(flux + flux0, 0.)
422
423        #Flux across upper edge of volume 1
424        normal = domain.get_normal(1,1)
425        ql = domain.get_conserved_quantities(vol_id=1, edge=1)
426        qr = domain.get_conserved_quantities(vol_id=3, edge=0)
427        flux1, max_speed = flux_function(normal, ql, qr, zl, zr)
428
429        #Check that flux seen from other triangles is inverse
430        tmp = qr; qr=ql; ql=tmp
431        normal = domain.get_normal(3,0)
432        flux, max_speed = flux_function(normal, ql, qr, zl, zr)
433        assert allclose(flux + flux1, 0.)
434
435        #Flux across lower left hypotenuse of volume 1
436        normal = domain.get_normal(1,2)
437        ql = domain.get_conserved_quantities(vol_id=1, edge=2)
438        qr = domain.get_conserved_quantities(vol_id=0, edge=1)
439        flux2, max_speed = flux_function(normal, ql, qr, zl, zr)
440
441        #Check that flux seen from other triangles is inverse
442        tmp = qr; qr=ql; ql=tmp
443        normal = domain.get_normal(0,1)
444        flux, max_speed = flux_function(normal, ql, qr, zl, zr)
445        assert allclose(flux + flux2, 0.)
446
447
448        #Scale by edgelengths, add up anc check that total flux is zero
449        e0 = domain.edgelengths[1, 0]
450        e1 = domain.edgelengths[1, 1]
451        e2 = domain.edgelengths[1, 2]
452
453        assert allclose(e0*flux0+e1*flux1+e2*flux2, 0.)
454
455        #Now check that compute_flux yields zeros as well
456        domain.compute_fluxes()
457
458        for name in ['stage', 'xmomentum', 'ymomentum']:
459            #print name, domain.quantities[name].explicit_update
460            assert allclose(domain.quantities[name].explicit_update[1], 0)
461
462
463
464    def test_compute_fluxes1(self):
465        #Use values from previous version
466
467        a = [0.0, 0.0]
468        b = [0.0, 2.0]
469        c = [2.0,0.0]
470        d = [0.0, 4.0]
471        e = [2.0, 2.0]
472        f = [4.0,0.0]
473
474        points = [a, b, c, d, e, f]
475        #bac, bce, ecf, dbe
476        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
477
478        domain = Domain(points, vertices)
479        val0 = 2.+2.0/3
480        val1 = 4.+4.0/3
481        val2 = 8.+2.0/3
482        val3 = 2.+8.0/3
483
484        domain.set_quantity('stage', [[val0, val0, val0], [val1, val1, val1],
485                                      [val2, val2, val2], [val3, val3, val3]])
486        domain.check_integrity()
487
488        zl=zr=0. #Assume flat bed
489
490        #Flux across right edge of volume 1
491        normal = domain.get_normal(1,0) #Get normal 0 of triangle 1
492        assert allclose(normal, [1, 0])
493       
494        ql = domain.get_conserved_quantities(vol_id=1, edge=0)
495        assert allclose(ql, [val1, 0, 0])
496       
497        qr = domain.get_conserved_quantities(vol_id=2, edge=2)
498        assert allclose(qr, [val2, 0, 0])
499       
500        flux0, max_speed = flux_function(normal, ql, qr, zl, zr)
501
502        #Flux across edge in the east direction (as per normal vector)
503        assert allclose(flux0, [-15.3598804, 253.71111111, 0.])
504        assert allclose(max_speed, 9.21592824046)
505
506
507        #Flux across edge in the west direction (opposite sign for xmomentum)
508        normal_opposite = domain.get_normal(2,2) #Get normal 2 of triangle 2
509        assert allclose(normal_opposite, [-1, 0])
510        flux_opposite, max_speed = flux_function([-1, 0], ql, qr, zl, zr)
511        assert allclose(flux_opposite, [-15.3598804, -253.71111111, 0.])
512       
513
514        #Flux across upper edge of volume 1
515        normal = domain.get_normal(1,1)
516        ql = domain.get_conserved_quantities(vol_id=1, edge=1)
517        qr = domain.get_conserved_quantities(vol_id=3, edge=0)
518        flux1, max_speed = flux_function(normal, ql, qr, zl, zr)
519        assert allclose(flux1, [2.4098563, 0., 123.04444444])
520        assert allclose(max_speed, 7.22956891292)
521
522        #Flux across lower left hypotenuse of volume 1
523        normal = domain.get_normal(1,2)
524        ql = domain.get_conserved_quantities(vol_id=1, edge=2)
525        qr = domain.get_conserved_quantities(vol_id=0, edge=1)
526        flux2, max_speed = flux_function(normal, ql, qr, zl, zr)
527
528        assert allclose(flux2, [9.63942522, -61.59685738, -61.59685738])
529        assert allclose(max_speed, 7.22956891292)
530
531        #Scale, add up and check that compute_fluxes is correct for vol 1
532        e0 = domain.edgelengths[1, 0]
533        e1 = domain.edgelengths[1, 1]
534        e2 = domain.edgelengths[1, 2]
535
536        total_flux = -(e0*flux0+e1*flux1+e2*flux2)/domain.areas[1]
537        assert allclose(total_flux, [-0.68218178, -166.6, -35.93333333])
538
539
540        domain.compute_fluxes()
541
542        #assert allclose(total_flux, domain.explicit_update[1,:])
543        for i, name in enumerate(['stage', 'xmomentum', 'ymomentum']):
544            assert allclose(total_flux[i],
545                            domain.quantities[name].explicit_update[1])
546
547        #assert allclose(domain.explicit_update, [
548        #    [0., -69.68888889, -69.68888889],
549        #    [-0.68218178, -166.6, -35.93333333],
550        #    [-111.77316251, 69.68888889, 0.],
551        #    [-35.68522449, 0., 69.68888889]])
552
553        assert allclose(domain.quantities['stage'].explicit_update,
554                        [0., -0.68218178, -111.77316251, -35.68522449])
555        assert allclose(domain.quantities['xmomentum'].explicit_update,
556                        [-69.68888889, -166.6, 69.68888889, 0])
557        assert allclose(domain.quantities['ymomentum'].explicit_update,
558                        [-69.68888889, -35.93333333, 0., 69.68888889])
559
560
561        #assert allclose(domain.quantities[name].explicit_update
562
563
564
565
566
567    def test_compute_fluxes2(self):
568        #Random values, incl momentum
569
570        a = [0.0, 0.0]
571        b = [0.0, 2.0]
572        c = [2.0,0.0]
573        d = [0.0, 4.0]
574        e = [2.0, 2.0]
575        f = [4.0,0.0]
576
577        points = [a, b, c, d, e, f]
578        #bac, bce, ecf, dbe
579        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
580
581        domain = Domain(points, vertices)
582        val0 = 2.+2.0/3
583        val1 = 4.+4.0/3
584        val2 = 8.+2.0/3
585        val3 = 2.+8.0/3
586
587        zl=zr=0 #Assume flat zero bed
588
589        domain.set_quantity('elevation', zl*ones( (4,3) ))
590
591
592        domain.set_quantity('stage', [[val0, val0-1, val0-2],
593                                      [val1, val1+1, val1],
594                                      [val2, val2-2, val2],
595                                      [val3-0.5, val3, val3]])
596
597        domain.set_quantity('xmomentum',
598                            [[1, 2, 3], [3, 4, 5],
599                             [1, -1, 0], [0, -2, 2]])
600
601        domain.set_quantity('ymomentum',
602                            [[1, -1, 0], [0, -3, 2],
603                             [0, 1, 0], [-1, 2, 2]])
604
605
606        domain.check_integrity()
607
608
609
610        #Flux across right edge of volume 1
611        normal = domain.get_normal(1,0)
612        ql = domain.get_conserved_quantities(vol_id=1, edge=0)
613        qr = domain.get_conserved_quantities(vol_id=2, edge=2)
614        flux0, max_speed = flux_function(normal, ql, qr, zl, zr)
615
616        #Flux across upper edge of volume 1
617        normal = domain.get_normal(1,1)
618        ql = domain.get_conserved_quantities(vol_id=1, edge=1)
619        qr = domain.get_conserved_quantities(vol_id=3, edge=0)
620        flux1, max_speed = flux_function(normal, ql, qr, zl, zr)
621
622        #Flux across lower left hypotenuse of volume 1
623        normal = domain.get_normal(1,2)
624        ql = domain.get_conserved_quantities(vol_id=1, edge=2)
625        qr = domain.get_conserved_quantities(vol_id=0, edge=1)
626        flux2, max_speed = flux_function(normal, ql, qr, zl, zr)
627
628        #Scale, add up and check that compute_fluxes is correct for vol 1
629        e0 = domain.edgelengths[1, 0]
630        e1 = domain.edgelengths[1, 1]
631        e2 = domain.edgelengths[1, 2]
632
633        total_flux = -(e0*flux0+e1*flux1+e2*flux2)/domain.areas[1]
634
635
636        domain.compute_fluxes()
637        for i, name in enumerate(['stage', 'xmomentum', 'ymomentum']):
638            assert allclose(total_flux[i],
639                            domain.quantities[name].explicit_update[1])
640        #assert allclose(total_flux, domain.explicit_update[1,:])
641
642
643    # FIXME (Ole): Need test like this for fluxes in very shallow water.   
644    def test_compute_fluxes3(self):
645        #Random values, incl momentum
646
647        a = [0.0, 0.0]
648        b = [0.0, 2.0]
649        c = [2.0,0.0]
650        d = [0.0, 4.0]
651        e = [2.0, 2.0]
652        f = [4.0,0.0]
653
654        points = [a, b, c, d, e, f]
655        #bac, bce, ecf, dbe
656        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
657
658        domain = Domain(points, vertices)
659        val0 = 2.+2.0/3
660        val1 = 4.+4.0/3
661        val2 = 8.+2.0/3
662        val3 = 2.+8.0/3
663
664        zl=zr=-3.75 #Assume constant bed (must be less than stage)
665        domain.set_quantity('elevation', zl*ones( (4,3) ))
666
667
668        domain.set_quantity('stage', [[val0, val0-1, val0-2],
669                                      [val1, val1+1, val1],
670                                      [val2, val2-2, val2],
671                                      [val3-0.5, val3, val3]])
672
673        domain.set_quantity('xmomentum',
674                            [[1, 2, 3], [3, 4, 5],
675                             [1, -1, 0], [0, -2, 2]])
676
677        domain.set_quantity('ymomentum',
678                            [[1, -1, 0], [0, -3, 2],
679                             [0, 1, 0], [-1, 2, 2]])
680
681
682        domain.check_integrity()
683
684
685
686        #Flux across right edge of volume 1
687        normal = domain.get_normal(1,0)
688        ql = domain.get_conserved_quantities(vol_id=1, edge=0)
689        qr = domain.get_conserved_quantities(vol_id=2, edge=2)
690        flux0, max_speed = flux_function(normal, ql, qr, zl, zr)
691
692        #Flux across upper edge of volume 1
693        normal = domain.get_normal(1,1)
694        ql = domain.get_conserved_quantities(vol_id=1, edge=1)
695        qr = domain.get_conserved_quantities(vol_id=3, edge=0)
696        flux1, max_speed = flux_function(normal, ql, qr, zl, zr)
697
698        #Flux across lower left hypotenuse of volume 1
699        normal = domain.get_normal(1,2)
700        ql = domain.get_conserved_quantities(vol_id=1, edge=2)
701        qr = domain.get_conserved_quantities(vol_id=0, edge=1)
702        flux2, max_speed = flux_function(normal, ql, qr, zl, zr)
703
704        #Scale, add up and check that compute_fluxes is correct for vol 1
705        e0 = domain.edgelengths[1, 0]
706        e1 = domain.edgelengths[1, 1]
707        e2 = domain.edgelengths[1, 2]
708
709        total_flux = -(e0*flux0+e1*flux1+e2*flux2)/domain.areas[1]
710
711        domain.compute_fluxes()
712        for i, name in enumerate(['stage', 'xmomentum', 'ymomentum']):
713            assert allclose(total_flux[i],
714                            domain.quantities[name].explicit_update[1])
715
716
717
718    def xtest_catching_negative_heights(self):
719
720        #OBSOLETE
721       
722        a = [0.0, 0.0]
723        b = [0.0, 2.0]
724        c = [2.0,0.0]
725        d = [0.0, 4.0]
726        e = [2.0, 2.0]
727        f = [4.0,0.0]
728
729        points = [a, b, c, d, e, f]
730        #bac, bce, ecf, dbe
731        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
732
733        domain = Domain(points, vertices)
734        val0 = 2.+2.0/3
735        val1 = 4.+4.0/3
736        val2 = 8.+2.0/3
737        val3 = 2.+8.0/3
738
739        zl=zr=4  #Too large
740        domain.set_quantity('elevation', zl*ones( (4,3) ))
741        domain.set_quantity('stage', [[val0, val0-1, val0-2],
742                                      [val1, val1+1, val1],
743                                      [val2, val2-2, val2],
744                                      [val3-0.5, val3, val3]])
745
746        #Should fail
747        try:
748            domain.check_integrity()
749        except:
750            pass
751
752
753
754    def test_get_wet_elements(self):
755
756        a = [0.0, 0.0]
757        b = [0.0, 2.0]
758        c = [2.0,0.0]
759        d = [0.0, 4.0]
760        e = [2.0, 2.0]
761        f = [4.0,0.0]
762
763        points = [a, b, c, d, e, f]
764        #bac, bce, ecf, dbe
765        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
766
767        domain = Domain(points, vertices)
768        val0 = 2.+2.0/3
769        val1 = 4.+4.0/3
770        val2 = 8.+2.0/3
771        val3 = 2.+8.0/3
772
773        zl=zr=5
774        domain.set_quantity('elevation', zl*ones( (4,3) ))
775        domain.set_quantity('stage', [[val0, val0-1, val0-2],
776                                      [val1, val1+1, val1],
777                                      [val2, val2-2, val2],
778                                      [val3-0.5, val3, val3]])
779
780
781
782        #print domain.get_quantity('elevation').get_values(location='centroids')
783        #print domain.get_quantity('stage').get_values(location='centroids')       
784        domain.check_integrity()
785
786        indices = domain.get_wet_elements()
787        assert allclose(indices, [1,2])
788
789        indices = domain.get_wet_elements(indices=[0,1,3])
790        assert allclose(indices, [1])
791       
792
793
794    def test_get_maximum_inundation_1(self):
795
796        a = [0.0, 0.0]
797        b = [0.0, 2.0]
798        c = [2.0,0.0]
799        d = [0.0, 4.0]
800        e = [2.0, 2.0]
801        f = [4.0,0.0]
802
803        points = [a, b, c, d, e, f]
804        #bac, bce, ecf, dbe
805        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
806
807        domain = Domain(points, vertices)
808
809        domain.set_quantity('elevation', lambda x, y: x+2*y) #2 4 4 6
810        domain.set_quantity('stage', 3)
811
812        domain.check_integrity()
813
814        indices = domain.get_wet_elements()
815        assert allclose(indices, [0])
816
817        q = domain.get_maximum_inundation_elevation()
818        assert allclose(q, domain.get_quantity('elevation').get_values(location='centroids')[0])
819
820        x, y = domain.get_maximum_inundation_location()
821        assert allclose([x, y], domain.get_centroid_coordinates()[0])
822
823
824    def test_get_maximum_inundation_2(self):
825        """test_get_maximum_inundation_2(self)
826
827        Test multiple wet cells with same elevation
828        """
829       
830        a = [0.0, 0.0]
831        b = [0.0, 2.0]
832        c = [2.0,0.0]
833        d = [0.0, 4.0]
834        e = [2.0, 2.0]
835        f = [4.0,0.0]
836
837        points = [a, b, c, d, e, f]
838        #bac, bce, ecf, dbe
839        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
840
841        domain = Domain(points, vertices)
842
843        domain.set_quantity('elevation', lambda x, y: x+2*y) #2 4 4 6
844        domain.set_quantity('stage', 4.1)
845
846        domain.check_integrity()
847
848        indices = domain.get_wet_elements()
849        assert allclose(indices, [0,1,2])
850
851        q = domain.get_maximum_inundation_elevation()
852        assert allclose(q, 4)       
853
854        x, y = domain.get_maximum_inundation_location()
855        assert allclose([x, y], domain.get_centroid_coordinates()[1])       
856       
857
858    def test_get_maximum_inundation_3(self):
859        """test_get_maximum_inundation_3(self)
860
861        Test of real runup example:
862        """
863
864        from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular_cross
865
866        initial_runup_height = -0.4
867        final_runup_height = -0.3
868
869
870        #--------------------------------------------------------------
871        # Setup computational domain
872        #--------------------------------------------------------------
873        N = 5
874        points, vertices, boundary = rectangular_cross(N, N) 
875        domain = Domain(points, vertices, boundary)           
876
877        #--------------------------------------------------------------
878        # Setup initial conditions
879        #--------------------------------------------------------------
880        def topography(x,y):
881            return -x/2                             # linear bed slope
882           
883
884        domain.set_quantity('elevation', topography)       # Use function for elevation
885        domain.set_quantity('friction', 0.)                # Zero friction
886        domain.set_quantity('stage', initial_runup_height) # Constant negative initial stage
887
888
889        #--------------------------------------------------------------
890        # Setup boundary conditions
891        #--------------------------------------------------------------
892        Br = Reflective_boundary(domain)              # Reflective wall
893        Bd = Dirichlet_boundary([final_runup_height,  # Constant inflow
894                                 0,
895                                 0])
896
897        # All reflective to begin with (still water)
898        domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
899
900
901        #--------------------------------------------------------------
902        # Test initial inundation height
903        #--------------------------------------------------------------
904
905        indices = domain.get_wet_elements()
906        z = domain.get_quantity('elevation').\
907            get_values(location='centroids', indices=indices)
908        assert alltrue(z<initial_runup_height)
909
910        q = domain.get_maximum_inundation_elevation()
911        assert allclose(q, initial_runup_height, rtol = 1.0/N) # First order accuracy
912
913        x, y = domain.get_maximum_inundation_location()
914
915        qref = domain.get_quantity('elevation').get_values(interpolation_points = [[x, y]])
916        assert allclose(q, qref)
917
918
919        wet_elements = domain.get_wet_elements()
920        wet_elevations = domain.get_quantity('elevation').get_values(location='centroids',
921                                                                     indices=wet_elements)
922        assert alltrue(wet_elevations<initial_runup_height)
923        assert allclose(wet_elevations, z)       
924
925
926        #print domain.get_quantity('elevation').get_maximum_value(indices=wet_elements)
927        #print domain.get_quantity('elevation').get_maximum_location(indices=wet_elements)
928        #print domain.get_quantity('elevation').get_maximum_index(indices=wet_elements)
929
930       
931        #--------------------------------------------------------------
932        # Let triangles adjust
933        #--------------------------------------------------------------
934        for t in domain.evolve(yieldstep = 0.1, finaltime = 1.0):
935            pass
936
937
938        #--------------------------------------------------------------
939        # Test inundation height again
940        #--------------------------------------------------------------
941
942        indices = domain.get_wet_elements()
943        z = domain.get_quantity('elevation').\
944            get_values(location='centroids', indices=indices)
945
946        assert alltrue(z<initial_runup_height)
947
948        q = domain.get_maximum_inundation_elevation()
949        assert allclose(q, initial_runup_height, rtol = 1.0/N) # First order accuracy
950       
951        x, y = domain.get_maximum_inundation_location()
952        qref = domain.get_quantity('elevation').get_values(interpolation_points = [[x, y]])
953        assert allclose(q, qref)       
954
955
956        #--------------------------------------------------------------
957        # Update boundary to allow inflow
958        #--------------------------------------------------------------
959        domain.set_boundary({'right': Bd})
960
961       
962        #--------------------------------------------------------------
963        # Evolve system through time
964        #--------------------------------------------------------------
965        for t in domain.evolve(yieldstep = 0.1, finaltime = 3.0):
966            #domain.write_time()
967            pass
968   
969        #--------------------------------------------------------------
970        # Test inundation height again
971        #--------------------------------------------------------------
972
973        indices = domain.get_wet_elements()
974        z = domain.get_quantity('elevation').\
975            get_values(location='centroids', indices=indices)
976
977        assert alltrue(z<final_runup_height)
978
979        q = domain.get_maximum_inundation_elevation()
980        assert allclose(q, final_runup_height, rtol = 1.0/N) # First order accuracy
981
982        x, y = domain.get_maximum_inundation_location()
983        qref = domain.get_quantity('elevation').get_values(interpolation_points = [[x, y]])
984        assert allclose(q, qref)       
985
986
987        wet_elements = domain.get_wet_elements()
988        wet_elevations = domain.get_quantity('elevation').get_values(location='centroids',
989                                                                     indices=wet_elements)
990        assert alltrue(wet_elevations<final_runup_height)
991        assert allclose(wet_elevations, z)       
992       
993
994
995    def test_get_maximum_inundation_from_sww(self):
996        """test_get_maximum_inundation_from_sww(self)
997
998        Test of get_maximum_inundation_elevation()
999        and get_maximum_inundation_location() from data_manager.py
1000       
1001        This is based on test_get_maximum_inundation_3(self) but works with the
1002        stored results instead of with the internal data structure.
1003
1004        This test uses the underlying get_maximum_inundation_data for tests
1005        """
1006
1007        from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular_cross
1008        from data_manager import get_maximum_inundation_elevation, get_maximum_inundation_location, get_maximum_inundation_data
1009       
1010
1011        initial_runup_height = -0.4
1012        final_runup_height = -0.3
1013
1014
1015        #--------------------------------------------------------------
1016        # Setup computational domain
1017        #--------------------------------------------------------------
1018        N = 10
1019        points, vertices, boundary = rectangular_cross(N, N) 
1020        domain = Domain(points, vertices, boundary)
1021        domain.set_name('runup_test')
1022        #domain.tight_slope_limiters = 1 #FIXME: This works better with old limiters
1023
1024        #--------------------------------------------------------------
1025        # Setup initial conditions
1026        #--------------------------------------------------------------
1027        def topography(x,y):
1028            return -x/2                             # linear bed slope
1029           
1030
1031        domain.set_quantity('elevation', topography)       # Use function for elevation
1032        domain.set_quantity('friction', 0.)                # Zero friction
1033        domain.set_quantity('stage', initial_runup_height) # Constant negative initial stage
1034
1035
1036        #--------------------------------------------------------------
1037        # Setup boundary conditions
1038        #--------------------------------------------------------------
1039        Br = Reflective_boundary(domain)              # Reflective wall
1040        Bd = Dirichlet_boundary([final_runup_height,  # Constant inflow
1041                                 0,
1042                                 0])
1043
1044        # All reflective to begin with (still water)
1045        domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
1046
1047
1048        #--------------------------------------------------------------
1049        # Test initial inundation height
1050        #--------------------------------------------------------------
1051
1052        indices = domain.get_wet_elements()
1053        z = domain.get_quantity('elevation').\
1054            get_values(location='centroids', indices=indices)
1055        assert alltrue(z<initial_runup_height)
1056
1057        q_ref = domain.get_maximum_inundation_elevation()
1058        assert allclose(q_ref, initial_runup_height, rtol = 1.0/N) # First order accuracy
1059
1060       
1061        #--------------------------------------------------------------
1062        # Let triangles adjust
1063        #--------------------------------------------------------------
1064        for t in domain.evolve(yieldstep = 0.1, finaltime = 1.0):
1065            pass
1066
1067
1068        #--------------------------------------------------------------
1069        # Test inundation height again
1070        #--------------------------------------------------------------
1071       
1072        q_ref = domain.get_maximum_inundation_elevation()
1073        q = get_maximum_inundation_elevation('runup_test.sww')
1074        msg = 'We got %f, should have been %f' %(q, q_ref)
1075        assert allclose(q, q_ref, rtol=1.0/N), msg
1076
1077
1078        q = get_maximum_inundation_elevation('runup_test.sww')
1079        msg = 'We got %f, should have been %f' %(q, initial_runup_height)
1080        assert allclose(q, initial_runup_height, rtol = 1.0/N), msg
1081
1082
1083        # Test error condition if time interval is out
1084        try:
1085            q = get_maximum_inundation_elevation('runup_test.sww',
1086                                                 time_interval=[2.0, 3.0])
1087        except ValueError:
1088            pass
1089        else:
1090            msg = 'should have caught wrong time interval'
1091            raise Exception, msg
1092
1093        # Check correct time interval
1094        q, loc = get_maximum_inundation_data('runup_test.sww',
1095                                             time_interval=[0.0, 3.0])       
1096        msg = 'We got %f, should have been %f' %(q, initial_runup_height)
1097        assert allclose(q, initial_runup_height, rtol = 1.0/N), msg
1098        assert allclose(-loc[0]/2, q) # From topography formula         
1099       
1100
1101        #--------------------------------------------------------------
1102        # Update boundary to allow inflow
1103        #--------------------------------------------------------------
1104        domain.set_boundary({'right': Bd})
1105
1106       
1107        #--------------------------------------------------------------
1108        # Evolve system through time
1109        #--------------------------------------------------------------
1110        q_max = None
1111        for t in domain.evolve(yieldstep = 0.1, finaltime = 3.0,
1112                               skip_initial_step = True): 
1113            q = domain.get_maximum_inundation_elevation()
1114            if q > q_max: q_max = q
1115
1116   
1117        #--------------------------------------------------------------
1118        # Test inundation height again
1119        #--------------------------------------------------------------
1120
1121        indices = domain.get_wet_elements()
1122        z = domain.get_quantity('elevation').\
1123            get_values(location='centroids', indices=indices)
1124
1125        assert alltrue(z<final_runup_height)
1126
1127        q = domain.get_maximum_inundation_elevation()
1128        assert allclose(q, final_runup_height, rtol = 1.0/N) # First order accuracy
1129
1130        q, loc = get_maximum_inundation_data('runup_test.sww', time_interval=[3.0, 3.0])
1131        msg = 'We got %f, should have been %f' %(q, final_runup_height)
1132        assert allclose(q, final_runup_height, rtol=1.0/N), msg
1133        #print 'loc', loc, q       
1134        assert allclose(-loc[0]/2, q) # From topography formula         
1135
1136        q = get_maximum_inundation_elevation('runup_test.sww')
1137        loc = get_maximum_inundation_location('runup_test.sww')       
1138        msg = 'We got %f, should have been %f' %(q, q_max)
1139        assert allclose(q, q_max, rtol=1.0/N), msg
1140        #print 'loc', loc, q
1141        assert allclose(-loc[0]/2, q) # From topography formula
1142
1143       
1144
1145        q = get_maximum_inundation_elevation('runup_test.sww', time_interval=[0, 3])
1146        msg = 'We got %f, should have been %f' %(q, q_max)
1147        assert allclose(q, q_max, rtol=1.0/N), msg
1148
1149
1150        # Check polygon mode
1151        polygon = [[0.3, 0.0], [0.9, 0.0], [0.9, 1.0], [0.3, 1.0]] # Runup region
1152        q = get_maximum_inundation_elevation('runup_test.sww',
1153                                             polygon = polygon,
1154                                             time_interval=[0, 3])
1155        msg = 'We got %f, should have been %f' %(q, q_max)
1156        assert allclose(q, q_max, rtol=1.0/N), msg
1157
1158       
1159        polygon = [[0.9, 0.0], [1.0, 0.0], [1.0, 1.0], [0.9, 1.0]] # Offshore region
1160        q, loc = get_maximum_inundation_data('runup_test.sww',
1161                                             polygon = polygon,
1162                                             time_interval=[0, 3])
1163        msg = 'We got %f, should have been %f' %(q, -0.475)
1164        assert allclose(q, -0.475, rtol=1.0/N), msg
1165        assert is_inside_polygon(loc, polygon)
1166        assert allclose(-loc[0]/2, q) # From topography formula         
1167
1168
1169        polygon = [[0.0, 0.0], [0.4, 0.0], [0.4, 1.0], [0.0, 1.0]] # Dry region
1170        q, loc = get_maximum_inundation_data('runup_test.sww',
1171                                             polygon = polygon,
1172                                             time_interval=[0, 3])
1173        msg = 'We got %s, should have been None' %(q)
1174        assert q is None, msg
1175        msg = 'We got %s, should have been None' %(loc)
1176        assert loc is None, msg       
1177
1178        # Check what happens if no time point is within interval
1179        try:
1180            q = get_maximum_inundation_elevation('runup_test.sww', time_interval=[2.75, 2.75])
1181        except AssertionError:
1182            pass
1183        else:
1184            msg = 'Time interval should have raised an exception'
1185            raise msg
1186           
1187
1188    def test_another_runup_example(self):
1189        """test_another_runup_example
1190
1191        Test runup example where actual timeseries at interpolated
1192        points are tested.
1193        """
1194
1195        #-----------------------------------------------------------------
1196        # Import necessary modules
1197        #-----------------------------------------------------------------
1198
1199        from anuga.pmesh.mesh_interface import create_mesh_from_regions
1200        from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular_cross
1201        from anuga.shallow_water import Domain
1202        from anuga.shallow_water import Reflective_boundary
1203        from anuga.shallow_water import Dirichlet_boundary
1204
1205
1206        #-----------------------------------------------------------------
1207        # Setup computational domain
1208        #-----------------------------------------------------------------
1209        points, vertices, boundary = rectangular_cross(10, 10) # Basic mesh
1210        domain = Domain(points, vertices, boundary) # Create domain
1211        domain.set_quantities_to_be_stored(None)
1212        domain.set_maximum_allowed_speed(100) #
1213       
1214        # FIXME (Ole): Need tests where this is commented out
1215        domain.tight_slope_limiters = 0 # Backwards compatibility (14/4/7)                 
1216        domain.H0 = 0 # Backwards compatibility (6/2/7)
1217        domain.beta_h = 0.2 # Backwards compatibility (14/2/7)
1218
1219        #-----------------------------------------------------------------
1220        # Setup initial conditions
1221        #-----------------------------------------------------------------
1222
1223        def topography(x,y): 
1224            return -x/2                              # linear bed slope
1225
1226        domain.set_quantity('elevation', topography) 
1227        domain.set_quantity('friction', 0.0)         
1228        domain.set_quantity('stage', expression='elevation')           
1229
1230
1231        #----------------------------------------------------------------
1232        # Setup boundary conditions
1233        #----------------------------------------------------------------
1234
1235        Br = Reflective_boundary(domain)      # Solid reflective wall
1236        Bd = Dirichlet_boundary([-0.2,0.,0.]) # Constant boundary values
1237        domain.set_boundary({'left': Br, 'right': Bd, 'top': Br, 'bottom': Br})
1238
1239
1240        #----------------------------------------------------------------
1241        # Evolve system through time
1242        #----------------------------------------------------------------
1243
1244        interpolation_points = [[0.4,0.5], [0.6,0.5], [0.8,0.5], [0.9,0.5]]
1245        gauge_values = []
1246        for _ in interpolation_points:
1247            gauge_values.append([])
1248
1249        time = []
1250        for t in domain.evolve(yieldstep = 0.1, finaltime = 5.0):
1251            # Record time series at known points
1252            time.append(domain.get_time())
1253           
1254            stage = domain.get_quantity('stage')
1255            w = stage.get_values(interpolation_points=interpolation_points)
1256           
1257            for i, _ in enumerate(interpolation_points):
1258                gauge_values[i].append(w[i])
1259
1260
1261        #print
1262        #print time
1263        #print
1264        #for i, (x,y) in enumerate(interpolation_points):
1265        #    print i, gauge_values[i]
1266        #    print
1267
1268        #Reference (nautilus 13/10/2006)
1269
1270        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]
1271
1272
1273        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]
1274
1275        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]
1276
1277        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]
1278       
1279        #FIXME (DSG):This is a hack so the anuga install, not precompiled
1280        # works on DSG's win2000, python 2.3
1281        #The problem is the gauge_values[X] are 52 long, not 51.
1282        #
1283        # This was probably fixed by Stephen in changeset:3804
1284        if len(gauge_values[0]) == 52: gauge_values[0].pop()
1285        if len(gauge_values[1]) == 52: gauge_values[1].pop()
1286        if len(gauge_values[2]) == 52: gauge_values[2].pop()
1287        if len(gauge_values[3]) == 52: gauge_values[3].pop()
1288
1289##         print len(G0), len(gauge_values[0])
1290##         print len(G1), len(gauge_values[1])
1291       
1292        #print array(gauge_values[0])-array(G0)
1293
1294       
1295       
1296        assert allclose(gauge_values[0], G0)
1297        assert allclose(gauge_values[1], G1)
1298        assert allclose(gauge_values[2], G2)
1299        assert allclose(gauge_values[3], G3)       
1300
1301
1302
1303
1304
1305
1306
1307    #####################################################
1308
1309    def test_flux_optimisation(self):
1310        """test_flux_optimisation
1311        Test that fluxes are correctly computed using
1312        dry and still cell exclusions
1313        """
1314
1315        from anuga.config import g
1316        import copy
1317
1318        a = [0.0, 0.0]
1319        b = [0.0, 2.0]
1320        c = [2.0, 0.0]
1321        d = [0.0, 4.0]
1322        e = [2.0, 2.0]
1323        f = [4.0, 0.0]
1324
1325        points = [a, b, c, d, e, f]
1326        #bac, bce, ecf, dbe
1327        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
1328
1329        domain = Domain(points, vertices)
1330
1331        #Set up for a gradient of (3,0) at mid triangle (bce)
1332        def slope(x, y):
1333            return 3*x
1334
1335        h = 0.1
1336        def stage(x,y):
1337            return slope(x,y)+h
1338
1339        domain.set_quantity('elevation', slope)
1340        domain.set_quantity('stage', stage)
1341
1342        # Allow slope limiters to work (FIXME (Ole): Shouldn't this be automatic in ANUGA?)     
1343        domain.distribute_to_vertices_and_edges()       
1344
1345        initial_stage = copy.copy(domain.quantities['stage'].vertex_values)
1346
1347        domain.set_boundary({'exterior': Reflective_boundary(domain)})
1348
1349
1350        #  Check that update arrays are initialised to zero
1351        assert allclose(domain.get_quantity('stage').explicit_update, 0)
1352        assert allclose(domain.get_quantity('xmomentum').explicit_update, 0)
1353        assert allclose(domain.get_quantity('ymomentum').explicit_update, 0)               
1354
1355
1356        # Get true values
1357        domain.optimise_dry_cells = False
1358        domain.compute_fluxes()
1359        stage_ref = copy.copy(domain.get_quantity('stage').explicit_update)
1360        xmom_ref = copy.copy(domain.get_quantity('xmomentum').explicit_update)
1361        ymom_ref = copy.copy(domain.get_quantity('ymomentum').explicit_update)       
1362
1363        # Try with flux optimisation
1364        domain.optimise_dry_cells = True
1365        domain.compute_fluxes()
1366
1367        assert allclose(stage_ref, domain.get_quantity('stage').explicit_update)
1368        assert allclose(xmom_ref, domain.get_quantity('xmomentum').explicit_update)
1369        assert allclose(ymom_ref, domain.get_quantity('ymomentum').explicit_update)
1370       
1371   
1372       
1373    def test_initial_condition(self):
1374        """test_initial_condition
1375        Test that initial condition is output at time == 0 and that
1376        computed values change as system evolves
1377        """
1378
1379        from anuga.config import g
1380        import copy
1381
1382        a = [0.0, 0.0]
1383        b = [0.0, 2.0]
1384        c = [2.0, 0.0]
1385        d = [0.0, 4.0]
1386        e = [2.0, 2.0]
1387        f = [4.0, 0.0]
1388
1389        points = [a, b, c, d, e, f]
1390        #bac, bce, ecf, dbe
1391        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
1392
1393        domain = Domain(points, vertices)
1394
1395        #Set up for a gradient of (3,0) at mid triangle (bce)
1396        def slope(x, y):
1397            return 3*x
1398
1399        h = 0.1
1400        def stage(x,y):
1401            return slope(x,y)+h
1402
1403        domain.set_quantity('elevation', slope)
1404        domain.set_quantity('stage', stage)
1405
1406        # Allow slope limiters to work (FIXME (Ole): Shouldn't this be automatic in ANUGA?)     
1407        domain.distribute_to_vertices_and_edges()       
1408
1409        initial_stage = copy.copy(domain.quantities['stage'].vertex_values)
1410
1411        domain.set_boundary({'exterior': Reflective_boundary(domain)})
1412
1413        domain.optimise_dry_cells = True
1414        #Evolution
1415        for t in domain.evolve(yieldstep = 0.5, finaltime = 2.0):
1416            stage = domain.quantities['stage'].vertex_values
1417
1418            if t == 0.0:
1419                assert allclose(stage, initial_stage)
1420            else:
1421                assert not allclose(stage, initial_stage)
1422
1423
1424        os.remove(domain.get_name() + '.sww')
1425
1426
1427
1428    #####################################################
1429    def test_gravity(self):
1430        #Assuming no friction
1431
1432        from anuga.config import g
1433
1434        a = [0.0, 0.0]
1435        b = [0.0, 2.0]
1436        c = [2.0, 0.0]
1437        d = [0.0, 4.0]
1438        e = [2.0, 2.0]
1439        f = [4.0, 0.0]
1440
1441        points = [a, b, c, d, e, f]
1442        #bac, bce, ecf, dbe
1443        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
1444
1445        domain = Domain(points, vertices)
1446
1447        #Set up for a gradient of (3,0) at mid triangle (bce)
1448        def slope(x, y):
1449            return 3*x
1450
1451        h = 0.1
1452        def stage(x,y):
1453            return slope(x,y)+h
1454
1455        domain.set_quantity('elevation', slope)
1456        domain.set_quantity('stage', stage)
1457
1458        for name in domain.conserved_quantities:
1459            assert allclose(domain.quantities[name].explicit_update, 0)
1460            assert allclose(domain.quantities[name].semi_implicit_update, 0)
1461
1462        domain.compute_forcing_terms()
1463
1464        assert allclose(domain.quantities['stage'].explicit_update, 0)
1465        assert allclose(domain.quantities['xmomentum'].explicit_update, -g*h*3)
1466        assert allclose(domain.quantities['ymomentum'].explicit_update, 0)
1467
1468
1469    def test_manning_friction(self):
1470        from anuga.config import g
1471
1472        a = [0.0, 0.0]
1473        b = [0.0, 2.0]
1474        c = [2.0, 0.0]
1475        d = [0.0, 4.0]
1476        e = [2.0, 2.0]
1477        f = [4.0, 0.0]
1478
1479        points = [a, b, c, d, e, f]
1480        #bac, bce, ecf, dbe
1481        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
1482
1483        domain = Domain(points, vertices)
1484
1485        #Set up for a gradient of (3,0) at mid triangle (bce)
1486        def slope(x, y):
1487            return 3*x
1488
1489        h = 0.1
1490        def stage(x,y):
1491            return slope(x,y)+h
1492
1493        eta = 0.07
1494        domain.set_quantity('elevation', slope)
1495        domain.set_quantity('stage', stage)
1496        domain.set_quantity('friction', eta)
1497
1498        for name in domain.conserved_quantities:
1499            assert allclose(domain.quantities[name].explicit_update, 0)
1500            assert allclose(domain.quantities[name].semi_implicit_update, 0)
1501
1502        domain.compute_forcing_terms()
1503
1504        assert allclose(domain.quantities['stage'].explicit_update, 0)
1505        assert allclose(domain.quantities['xmomentum'].explicit_update, -g*h*3)
1506        assert allclose(domain.quantities['ymomentum'].explicit_update, 0)
1507
1508        assert allclose(domain.quantities['stage'].semi_implicit_update, 0)
1509        assert allclose(domain.quantities['xmomentum'].semi_implicit_update, 0)
1510        assert allclose(domain.quantities['ymomentum'].semi_implicit_update, 0)
1511
1512        #Create some momentum for friction to work with
1513        domain.set_quantity('xmomentum', 1)
1514        S = -g * eta**2 / h**(7.0/3)
1515
1516        domain.compute_forcing_terms()
1517        assert allclose(domain.quantities['stage'].semi_implicit_update, 0)
1518        assert allclose(domain.quantities['xmomentum'].semi_implicit_update, S)
1519        assert allclose(domain.quantities['ymomentum'].semi_implicit_update, 0)
1520
1521        #A more complex example
1522        domain.quantities['stage'].semi_implicit_update[:] = 0.0
1523        domain.quantities['xmomentum'].semi_implicit_update[:] = 0.0
1524        domain.quantities['ymomentum'].semi_implicit_update[:] = 0.0
1525
1526        domain.set_quantity('xmomentum', 3)
1527        domain.set_quantity('ymomentum', 4)
1528
1529        S = -g * eta**2 * 5 / h**(7.0/3)
1530
1531
1532        domain.compute_forcing_terms()
1533
1534        assert allclose(domain.quantities['stage'].semi_implicit_update, 0)
1535        assert allclose(domain.quantities['xmomentum'].semi_implicit_update, 3*S)
1536        assert allclose(domain.quantities['ymomentum'].semi_implicit_update, 4*S)
1537
1538    def test_constant_wind_stress(self):
1539        from anuga.config import rho_a, rho_w, eta_w
1540        from math import pi, cos, sin, sqrt
1541
1542        a = [0.0, 0.0]
1543        b = [0.0, 2.0]
1544        c = [2.0, 0.0]
1545        d = [0.0, 4.0]
1546        e = [2.0, 2.0]
1547        f = [4.0, 0.0]
1548
1549        points = [a, b, c, d, e, f]
1550        #bac, bce, ecf, dbe
1551        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
1552
1553
1554        domain = Domain(points, vertices)
1555
1556        #Flat surface with 1m of water
1557        domain.set_quantity('elevation', 0)
1558        domain.set_quantity('stage', 1.0)
1559        domain.set_quantity('friction', 0)
1560
1561        Br = Reflective_boundary(domain)
1562        domain.set_boundary({'exterior': Br})
1563
1564        #Setup only one forcing term, constant wind stress
1565        s = 100
1566        phi = 135
1567        domain.forcing_terms = []
1568        domain.forcing_terms.append( Wind_stress(s, phi) )
1569
1570        domain.compute_forcing_terms()
1571
1572
1573        const = eta_w*rho_a/rho_w
1574
1575        #Convert to radians
1576        phi = phi*pi/180
1577
1578        #Compute velocity vector (u, v)
1579        u = s*cos(phi)
1580        v = s*sin(phi)
1581
1582        #Compute wind stress
1583        S = const * sqrt(u**2 + v**2)
1584
1585        assert allclose(domain.quantities['stage'].explicit_update, 0)
1586        assert allclose(domain.quantities['xmomentum'].explicit_update, S*u)
1587        assert allclose(domain.quantities['ymomentum'].explicit_update, S*v)
1588
1589
1590    def test_variable_wind_stress(self):
1591        from anuga.config import rho_a, rho_w, eta_w
1592        from math import pi, cos, sin, sqrt
1593
1594        a = [0.0, 0.0]
1595        b = [0.0, 2.0]
1596        c = [2.0, 0.0]
1597        d = [0.0, 4.0]
1598        e = [2.0, 2.0]
1599        f = [4.0, 0.0]
1600
1601        points = [a, b, c, d, e, f]
1602        #bac, bce, ecf, dbe
1603        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
1604
1605        domain = Domain(points, vertices)
1606
1607        #Flat surface with 1m of water
1608        domain.set_quantity('elevation', 0)
1609        domain.set_quantity('stage', 1.0)
1610        domain.set_quantity('friction', 0)
1611
1612        Br = Reflective_boundary(domain)
1613        domain.set_boundary({'exterior': Br})
1614
1615
1616        domain.time = 5.54 #Take a random time (not zero)
1617
1618        #Setup only one forcing term, constant wind stress
1619        s = 100
1620        phi = 135
1621        domain.forcing_terms = []
1622        domain.forcing_terms.append( Wind_stress(s = speed, phi = angle) )
1623
1624        domain.compute_forcing_terms()
1625
1626        #Compute reference solution
1627        const = eta_w*rho_a/rho_w
1628
1629        N = len(domain) # number_of_triangles
1630
1631        xc = domain.get_centroid_coordinates()
1632        t = domain.time
1633
1634        x = xc[:,0]
1635        y = xc[:,1]
1636        s_vec = speed(t,x,y)
1637        phi_vec = angle(t,x,y)
1638
1639
1640        for k in range(N):
1641            #Convert to radians
1642            phi = phi_vec[k]*pi/180
1643            s = s_vec[k]
1644
1645            #Compute velocity vector (u, v)
1646            u = s*cos(phi)
1647            v = s*sin(phi)
1648
1649            #Compute wind stress
1650            S = const * sqrt(u**2 + v**2)
1651
1652            assert allclose(domain.quantities['stage'].explicit_update[k], 0)
1653            assert allclose(domain.quantities['xmomentum'].explicit_update[k], S*u)
1654            assert allclose(domain.quantities['ymomentum'].explicit_update[k], S*v)
1655
1656
1657
1658
1659    def test_windfield_from_file(self):
1660        from anuga.config import rho_a, rho_w, eta_w
1661        from math import pi, cos, sin, sqrt
1662        from anuga.config import time_format
1663        from anuga.abstract_2d_finite_volumes.util import file_function
1664        import time
1665
1666
1667        a = [0.0, 0.0]
1668        b = [0.0, 2.0]
1669        c = [2.0, 0.0]
1670        d = [0.0, 4.0]
1671        e = [2.0, 2.0]
1672        f = [4.0, 0.0]
1673
1674        points = [a, b, c, d, e, f]
1675        #bac, bce, ecf, dbe
1676        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
1677
1678        domain = Domain(points, vertices)
1679
1680        #Flat surface with 1m of water
1681        domain.set_quantity('elevation', 0)
1682        domain.set_quantity('stage', 1.0)
1683        domain.set_quantity('friction', 0)
1684
1685        Br = Reflective_boundary(domain)
1686        domain.set_boundary({'exterior': Br})
1687
1688
1689        domain.time = 7 #Take a time that is represented in file (not zero)
1690
1691        #Write wind stress file (ensure that domain.time is covered)
1692        #Take x=1 and y=0
1693        filename = 'test_windstress_from_file'
1694        start = time.mktime(time.strptime('2000', '%Y'))
1695        fid = open(filename + '.txt', 'w')
1696        dt = 1  #One second interval
1697        t = 0.0
1698        while t <= 10.0:
1699            t_string = time.strftime(time_format, time.gmtime(t+start))
1700
1701            fid.write('%s, %f %f\n' %(t_string,
1702                                      speed(t,[1],[0])[0],
1703                                      angle(t,[1],[0])[0]))
1704            t += dt
1705
1706        fid.close()
1707
1708
1709        #Convert ASCII file to NetCDF (Which is what we really like!)
1710        from data_manager import timefile2netcdf       
1711        timefile2netcdf(filename)
1712        os.remove(filename + '.txt')
1713
1714       
1715        #Setup wind stress
1716        F = file_function(filename + '.tms', quantities = ['Attribute0',
1717                                                           'Attribute1'])
1718        os.remove(filename + '.tms')
1719       
1720
1721        #print 'F(5)', F(5)
1722
1723        #print 'F(5,x,y)', F(5,x=zeros(3),y=zeros(3))
1724       
1725        #print dir(F)
1726        #print F.T
1727        #print F.precomputed_values
1728        #
1729        #F = file_function(filename + '.txt')       
1730        #
1731        #print dir(F)
1732        #print F.T       
1733        #print F.Q
1734       
1735        W = Wind_stress(F)
1736
1737        domain.forcing_terms = []
1738        domain.forcing_terms.append(W)
1739
1740        domain.compute_forcing_terms()
1741
1742        #Compute reference solution
1743        const = eta_w*rho_a/rho_w
1744
1745        N = len(domain) # number_of_triangles
1746
1747        t = domain.time
1748
1749        s = speed(t,[1],[0])[0]
1750        phi = angle(t,[1],[0])[0]
1751
1752        #Convert to radians
1753        phi = phi*pi/180
1754
1755
1756        #Compute velocity vector (u, v)
1757        u = s*cos(phi)
1758        v = s*sin(phi)
1759
1760        #Compute wind stress
1761        S = const * sqrt(u**2 + v**2)
1762
1763        for k in range(N):
1764            assert allclose(domain.quantities['stage'].explicit_update[k], 0)
1765            assert allclose(domain.quantities['xmomentum'].explicit_update[k], S*u)
1766            assert allclose(domain.quantities['ymomentum'].explicit_update[k], S*v)
1767
1768
1769    def test_windfield_from_file_seconds(self):
1770        from anuga.config import rho_a, rho_w, eta_w
1771        from math import pi, cos, sin, sqrt
1772        from anuga.config import time_format
1773        from anuga.abstract_2d_finite_volumes.util import file_function
1774        import time
1775
1776
1777        a = [0.0, 0.0]
1778        b = [0.0, 2.0]
1779        c = [2.0, 0.0]
1780        d = [0.0, 4.0]
1781        e = [2.0, 2.0]
1782        f = [4.0, 0.0]
1783
1784        points = [a, b, c, d, e, f]
1785        #bac, bce, ecf, dbe
1786        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
1787
1788        domain = Domain(points, vertices)
1789
1790        #Flat surface with 1m of water
1791        domain.set_quantity('elevation', 0)
1792        domain.set_quantity('stage', 1.0)
1793        domain.set_quantity('friction', 0)
1794
1795        Br = Reflective_boundary(domain)
1796        domain.set_boundary({'exterior': Br})
1797
1798
1799        domain.time = 7 #Take a time that is represented in file (not zero)
1800
1801        #Write wind stress file (ensure that domain.time is covered)
1802        #Take x=1 and y=0
1803        filename = 'test_windstress_from_file'
1804        start = time.mktime(time.strptime('2000', '%Y'))
1805        fid = open(filename + '.txt', 'w')
1806        dt = 0.5 #1  #One second interval
1807        t = 0.0
1808        while t <= 10.0:
1809            fid.write('%s, %f %f\n' %(str(t),
1810                                      speed(t,[1],[0])[0],
1811                                      angle(t,[1],[0])[0]))
1812            t += dt
1813
1814        fid.close()
1815
1816
1817        #Convert ASCII file to NetCDF (Which is what we really like!)
1818        from data_manager import timefile2netcdf       
1819        timefile2netcdf(filename, time_as_seconds=True)
1820        os.remove(filename + '.txt')
1821
1822       
1823        #Setup wind stress
1824        F = file_function(filename + '.tms', quantities = ['Attribute0',
1825                                                           'Attribute1'])
1826        os.remove(filename + '.tms')
1827       
1828
1829        #print 'F(5)', F(5)
1830
1831        #print 'F(5,x,y)', F(5,x=zeros(3),y=zeros(3))
1832       
1833        #print dir(F)
1834        #print F.T
1835        #print F.precomputed_values
1836        #
1837        #F = file_function(filename + '.txt')       
1838        #
1839        #print dir(F)
1840        #print F.T       
1841        #print F.Q
1842       
1843        W = Wind_stress(F)
1844
1845        domain.forcing_terms = []
1846        domain.forcing_terms.append(W)
1847
1848        domain.compute_forcing_terms()
1849
1850        #Compute reference solution
1851        const = eta_w*rho_a/rho_w
1852
1853        N = len(domain) # number_of_triangles
1854
1855        t = domain.time
1856
1857        s = speed(t,[1],[0])[0]
1858        phi = angle(t,[1],[0])[0]
1859
1860        #Convert to radians
1861        phi = phi*pi/180
1862
1863
1864        #Compute velocity vector (u, v)
1865        u = s*cos(phi)
1866        v = s*sin(phi)
1867
1868        #Compute wind stress
1869        S = const * sqrt(u**2 + v**2)
1870
1871        for k in range(N):
1872            assert allclose(domain.quantities['stage'].explicit_update[k], 0)
1873            assert allclose(domain.quantities['xmomentum'].explicit_update[k], S*u)
1874            assert allclose(domain.quantities['ymomentum'].explicit_update[k], S*v)
1875
1876
1877       
1878
1879    def test_wind_stress_error_condition(self):
1880        """Test that windstress reacts properly when forcing functions
1881        are wrong - e.g. returns a scalar
1882        """
1883
1884        from anuga.config import rho_a, rho_w, eta_w
1885        from math import pi, cos, sin, sqrt
1886
1887        a = [0.0, 0.0]
1888        b = [0.0, 2.0]
1889        c = [2.0, 0.0]
1890        d = [0.0, 4.0]
1891        e = [2.0, 2.0]
1892        f = [4.0, 0.0]
1893
1894        points = [a, b, c, d, e, f]
1895        #bac, bce, ecf, dbe
1896        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
1897
1898        domain = Domain(points, vertices)
1899
1900        #Flat surface with 1m of water
1901        domain.set_quantity('elevation', 0)
1902        domain.set_quantity('stage', 1.0)
1903        domain.set_quantity('friction', 0)
1904
1905        Br = Reflective_boundary(domain)
1906        domain.set_boundary({'exterior': Br})
1907
1908
1909        domain.time = 5.54 #Take a random time (not zero)
1910
1911        #Setup only one forcing term, bad func
1912        domain.forcing_terms = []
1913
1914        try:
1915            domain.forcing_terms.append(Wind_stress(s = scalar_func,
1916                                                    phi = angle))
1917        except AssertionError:
1918            pass
1919        else:
1920            msg = 'Should have raised exception'
1921            raise msg
1922
1923
1924        try:
1925            domain.forcing_terms.append(Wind_stress(s = speed,
1926                                                    phi = scalar_func))
1927        except AssertionError:
1928            pass
1929        else:
1930            msg = 'Should have raised exception'
1931            raise msg
1932
1933        try:
1934            domain.forcing_terms.append(Wind_stress(s = speed,
1935                                                    phi = 'xx'))
1936        except:
1937            pass
1938        else:
1939            msg = 'Should have raised exception'
1940            raise msg
1941
1942
1943    #####################################################
1944    def test_first_order_extrapolator_const_z(self):
1945
1946        a = [0.0, 0.0]
1947        b = [0.0, 2.0]
1948        c = [2.0, 0.0]
1949        d = [0.0, 4.0]
1950        e = [2.0, 2.0]
1951        f = [4.0, 0.0]
1952
1953        points = [a, b, c, d, e, f]
1954        #bac, bce, ecf, dbe
1955        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
1956
1957        domain = Domain(points, vertices)
1958        val0 = 2.+2.0/3
1959        val1 = 4.+4.0/3
1960        val2 = 8.+2.0/3
1961        val3 = 2.+8.0/3
1962
1963        zl=zr=-3.75 #Assume constant bed (must be less than stage)
1964        domain.set_quantity('elevation', zl*ones( (4,3) ))
1965        domain.set_quantity('stage', [[val0, val0-1, val0-2],
1966                                      [val1, val1+1, val1],
1967                                      [val2, val2-2, val2],
1968                                      [val3-0.5, val3, val3]])
1969
1970
1971
1972        domain._order_ = 1
1973        domain.distribute_to_vertices_and_edges()
1974
1975        #Check that centroid values were distributed to vertices
1976        C = domain.quantities['stage'].centroid_values
1977        for i in range(3):
1978            assert allclose( domain.quantities['stage'].vertex_values[:,i], C)
1979
1980
1981    def test_first_order_limiter_variable_z(self):
1982        #Check that first order limiter follows bed_slope
1983        from Numeric import alltrue, greater_equal
1984        from anuga.config import epsilon
1985
1986        a = [0.0, 0.0]
1987        b = [0.0, 2.0]
1988        c = [2.0,0.0]
1989        d = [0.0, 4.0]
1990        e = [2.0, 2.0]
1991        f = [4.0,0.0]
1992
1993        points = [a, b, c, d, e, f]
1994        #bac, bce, ecf, dbe
1995        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
1996
1997        domain = Domain(points, vertices)
1998        val0 = 2.+2.0/3
1999        val1 = 4.+4.0/3
2000        val2 = 8.+2.0/3
2001        val3 = 2.+8.0/3
2002
2003        domain.set_quantity('elevation', [[0,0,0], [6,0,0],
2004                                          [6,6,6], [6,6,6]])
2005        domain.set_quantity('stage', [[val0, val0, val0],
2006                                      [val1, val1, val1],
2007                                      [val2, val2, val2],
2008                                      [val3, val3, val3]])
2009
2010        E = domain.quantities['elevation'].vertex_values
2011        L = domain.quantities['stage'].vertex_values
2012
2013
2014        #Check that some stages are not above elevation (within eps)
2015        #- so that the limiter has something to work with
2016        assert not alltrue(alltrue(greater_equal(L,E-epsilon)))
2017
2018        domain._order_ = 1
2019        domain.distribute_to_vertices_and_edges()
2020
2021        #Check that all stages are above elevation (within eps)
2022        assert alltrue(alltrue(greater_equal(L,E-epsilon)))
2023
2024
2025    #####################################################
2026    def test_distribute_basic(self):
2027        #Using test data generated by abstract_2d_finite_volumes-2
2028        #Assuming no friction and flat bed (0.0)
2029
2030        a = [0.0, 0.0]
2031        b = [0.0, 2.0]
2032        c = [2.0, 0.0]
2033        d = [0.0, 4.0]
2034        e = [2.0, 2.0]
2035        f = [4.0, 0.0]
2036
2037        points = [a, b, c, d, e, f]
2038        #bac, bce, ecf, dbe
2039        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
2040
2041        domain = Domain(points, vertices)
2042
2043        val0 = 2.
2044        val1 = 4.
2045        val2 = 8.
2046        val3 = 2.
2047
2048        domain.set_quantity('stage', [val0, val1, val2, val3],
2049                            location='centroids')
2050        L = domain.quantities['stage'].vertex_values
2051
2052        #First order
2053        domain._order_ = 1
2054        domain.distribute_to_vertices_and_edges()
2055        assert allclose(L[1], val1)
2056
2057        #Second order
2058        domain._order_ = 2
2059        domain.beta_w      = 0.9
2060        domain.beta_w_dry  = 0.9
2061        domain.beta_uh     = 0.9
2062        domain.beta_uh_dry = 0.9
2063        domain.beta_vh     = 0.9
2064        domain.beta_vh_dry = 0.9
2065        domain.distribute_to_vertices_and_edges()
2066        assert allclose(L[1], [2.2, 4.9, 4.9])
2067
2068
2069
2070    def test_distribute_away_from_bed(self):
2071        #Using test data generated by abstract_2d_finite_volumes-2
2072        #Assuming no friction and flat bed (0.0)
2073
2074        a = [0.0, 0.0]
2075        b = [0.0, 2.0]
2076        c = [2.0, 0.0]
2077        d = [0.0, 4.0]
2078        e = [2.0, 2.0]
2079        f = [4.0, 0.0]
2080
2081        points = [a, b, c, d, e, f]
2082        #bac, bce, ecf, dbe
2083        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
2084
2085        domain = Domain(points, vertices)
2086        L = domain.quantities['stage'].vertex_values
2087
2088        def stage(x,y):
2089            return x**2
2090
2091        domain.set_quantity('stage', stage, location='centroids')
2092
2093        a, b = domain.quantities['stage'].compute_gradients()
2094        assert allclose(a[1], 3.33333334)
2095        assert allclose(b[1], 0.0)
2096
2097        domain._order_ = 1
2098        domain.distribute_to_vertices_and_edges()
2099        assert allclose(L[1], 1.77777778)
2100
2101        domain._order_ = 2
2102        domain.beta_w      = 0.9
2103        domain.beta_w_dry  = 0.9
2104        domain.beta_uh     = 0.9
2105        domain.beta_uh_dry = 0.9
2106        domain.beta_vh     = 0.9
2107        domain.beta_vh_dry = 0.9
2108        domain.distribute_to_vertices_and_edges()
2109        assert allclose(L[1], [0.57777777, 2.37777778, 2.37777778])
2110
2111
2112
2113    def test_distribute_away_from_bed1(self):
2114        #Using test data generated by abstract_2d_finite_volumes-2
2115        #Assuming no friction and flat bed (0.0)
2116
2117        a = [0.0, 0.0]
2118        b = [0.0, 2.0]
2119        c = [2.0, 0.0]
2120        d = [0.0, 4.0]
2121        e = [2.0, 2.0]
2122        f = [4.0, 0.0]
2123
2124        points = [a, b, c, d, e, f]
2125        #bac, bce, ecf, dbe
2126        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
2127
2128        domain = Domain(points, vertices)
2129        L = domain.quantities['stage'].vertex_values
2130
2131        def stage(x,y):
2132            return x**4+y**2
2133
2134        domain.set_quantity('stage', stage, location='centroids')
2135        #print domain.quantities['stage'].centroid_values
2136
2137        a, b = domain.quantities['stage'].compute_gradients()
2138        assert allclose(a[1], 25.18518519)
2139        assert allclose(b[1], 3.33333333)
2140
2141        domain._order_ = 1
2142        domain.distribute_to_vertices_and_edges()
2143        assert allclose(L[1], 4.9382716)
2144
2145        domain._order_ = 2
2146        domain.beta_w      = 0.9
2147        domain.beta_w_dry  = 0.9
2148        domain.beta_uh     = 0.9
2149        domain.beta_uh_dry = 0.9
2150        domain.beta_vh     = 0.9
2151        domain.beta_vh_dry = 0.9
2152        domain.distribute_to_vertices_and_edges()
2153        assert allclose(L[1], [1.07160494, 6.46058131, 7.28262855])
2154
2155
2156
2157    def test_distribute_near_bed(self):
2158
2159        a = [0.0, 0.0]
2160        b = [0.0, 2.0]
2161        c = [2.0, 0.0]
2162        d = [0.0, 4.0]
2163        e = [2.0, 2.0]
2164        f = [4.0, 0.0]
2165
2166        points = [a, b, c, d, e, f]
2167        #bac, bce, ecf, dbe
2168        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
2169
2170        domain = Domain(points, vertices)
2171
2172
2173        #Set up for a gradient of (10,0) at mid triangle (bce)
2174        def slope(x, y):
2175            return 10*x
2176
2177        h = 0.1
2178        def stage(x, y):
2179            return slope(x, y) + h
2180
2181        domain.set_quantity('elevation', slope)
2182        domain.set_quantity('stage', stage, location='centroids')
2183
2184        #print domain.quantities['elevation'].centroid_values
2185        #print domain.quantities['stage'].centroid_values
2186
2187        E = domain.quantities['elevation'].vertex_values
2188        L = domain.quantities['stage'].vertex_values
2189
2190        # Get reference values
2191        volumes = []
2192        for i in range(len(L)):
2193            volumes.append(sum(L[i])/3)
2194            assert allclose(volumes[i], domain.quantities['stage'].centroid_values[i]) 
2195       
2196       
2197        domain._order_ = 1
2198       
2199        domain.tight_slope_limiters = 0
2200        domain.distribute_to_vertices_and_edges()
2201        assert allclose(L[1], [0.1, 20.1, 20.1])
2202        for i in range(len(L)):
2203            assert allclose(volumes[i], sum(L[i])/3)                   
2204       
2205        domain.tight_slope_limiters = 1 # Allow triangle to be flatter (closer to bed)
2206        domain.distribute_to_vertices_and_edges()
2207        assert allclose(L[1], [0.298, 20.001, 20.001])
2208        for i in range(len(L)):
2209            assert allclose(volumes[i], sum(L[i])/3)   
2210
2211        domain._order_ = 2
2212       
2213        domain.tight_slope_limiters = 0
2214        domain.distribute_to_vertices_and_edges()
2215        assert allclose(L[1], [0.1, 20.1, 20.1])       
2216        for i in range(len(L)):
2217            assert allclose(volumes[i], sum(L[i])/3)           
2218       
2219        domain.tight_slope_limiters = 1 # Allow triangle to be flatter (closer to bed)
2220        domain.distribute_to_vertices_and_edges()
2221        assert allclose(L[1], [0.298, 20.001, 20.001])
2222        for i in range(len(L)):
2223            assert allclose(volumes[i], sum(L[i])/3)   
2224       
2225
2226
2227    def test_distribute_near_bed1(self):
2228
2229        a = [0.0, 0.0]
2230        b = [0.0, 2.0]
2231        c = [2.0, 0.0]
2232        d = [0.0, 4.0]
2233        e = [2.0, 2.0]
2234        f = [4.0, 0.0]
2235
2236        points = [a, b, c, d, e, f]
2237        #bac, bce, ecf, dbe
2238        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
2239
2240        domain = Domain(points, vertices)
2241
2242
2243        #Set up for a gradient of (8,2) at mid triangle (bce)
2244        def slope(x, y):
2245            return x**4+y**2
2246
2247        h = 0.1
2248        def stage(x,y):
2249            return slope(x,y)+h
2250
2251        domain.set_quantity('elevation', slope)
2252        domain.set_quantity('stage', stage)
2253
2254        #print domain.quantities['elevation'].centroid_values
2255        #print domain.quantities['stage'].centroid_values
2256
2257        E = domain.quantities['elevation'].vertex_values
2258        L = domain.quantities['stage'].vertex_values
2259
2260        # Get reference values
2261        volumes = []
2262        for i in range(len(L)):
2263            volumes.append(sum(L[i])/3)
2264            assert allclose(volumes[i], domain.quantities['stage'].centroid_values[i]) 
2265       
2266        #print E
2267        domain._order_ = 1
2268       
2269        domain.tight_slope_limiters = 0
2270        domain.distribute_to_vertices_and_edges()
2271        assert allclose(L[1], [4.1, 16.1, 20.1])       
2272        for i in range(len(L)):
2273            assert allclose(volumes[i], sum(L[i])/3)
2274       
2275               
2276        domain.tight_slope_limiters = 1 # Allow triangle to be flatter (closer to bed)
2277        domain.distribute_to_vertices_and_edges()
2278        assert allclose(L[1], [4.2386, 16.0604, 20.001])
2279        for i in range(len(L)):
2280            assert allclose(volumes[i], sum(L[i])/3)   
2281       
2282
2283        domain._order_ = 2
2284       
2285        domain.tight_slope_limiters = 0   
2286        domain.distribute_to_vertices_and_edges()
2287        assert allclose(L[1], [4.1, 16.1, 20.1])
2288        for i in range(len(L)):
2289            assert allclose(volumes[i], sum(L[i])/3)   
2290       
2291        domain.tight_slope_limiters = 1 # Allow triangle to be flatter (closer to bed)
2292        domain.distribute_to_vertices_and_edges()
2293        assert allclose(L[1], [4.23370103, 16.06529897, 20.001])
2294        for i in range(len(L)):
2295            assert allclose(volumes[i], sum(L[i])/3)
2296
2297
2298    def test_second_order_distribute_real_data(self):
2299        #Using test data generated by abstract_2d_finite_volumes-2
2300        #Assuming no friction and flat bed (0.0)
2301
2302        a = [0.0, 0.0]
2303        b = [0.0, 1.0/5]
2304        c = [0.0, 2.0/5]
2305        d = [1.0/5, 0.0]
2306        e = [1.0/5, 1.0/5]
2307        f = [1.0/5, 2.0/5]
2308        g = [2.0/5, 2.0/5]
2309
2310        points = [a, b, c, d, e, f, g]
2311        #bae, efb, cbf, feg
2312        vertices = [ [1,0,4], [4,5,1], [2,1,5], [5,4,6]]
2313
2314        domain = Domain(points, vertices)
2315
2316        def slope(x, y):
2317            return -x/3
2318
2319        domain.set_quantity('elevation', slope)
2320        domain.set_quantity('stage',
2321                            [0.01298164, 0.00365611,
2322                             0.01440365, -0.0381856437096],
2323                            location='centroids')
2324        domain.set_quantity('xmomentum',
2325                            [0.00670439, 0.01263789,
2326                             0.00647805, 0.0178180740668],
2327                            location='centroids')
2328        domain.set_quantity('ymomentum',
2329                            [-7.23510980e-004, -6.30413883e-005,
2330                             6.30413883e-005, 0.000200907255866],
2331                            location='centroids')
2332
2333        E = domain.quantities['elevation'].vertex_values
2334        L = domain.quantities['stage'].vertex_values
2335        X = domain.quantities['xmomentum'].vertex_values
2336        Y = domain.quantities['ymomentum'].vertex_values
2337
2338        #print E
2339        domain._order_ = 2
2340        domain.beta_w      = 0.9
2341        domain.beta_w_dry  = 0.9
2342        domain.beta_uh     = 0.9
2343        domain.beta_uh_dry = 0.9
2344        domain.beta_vh     = 0.9
2345        domain.beta_vh_dry = 0.9
2346        domain.beta_h = 0.0 #Use first order in h-limiter
2347       
2348        # FIXME (Ole): Need tests where this is commented out
2349        domain.tight_slope_limiters = 0 # Backwards compatibility (14/4/7)                 
2350       
2351               
2352        domain.distribute_to_vertices_and_edges()
2353
2354        #print L[1,:]
2355        #print X[1,:]
2356        #print Y[1,:]
2357
2358        assert allclose(L[1,:], [-0.00825735775384,
2359                                 -0.00801881482869,
2360                                 0.0272445025825])
2361        assert allclose(X[1,:], [0.0143507718962,
2362                                 0.0142502147066,
2363                                 0.00931268339717])
2364        assert allclose(Y[1,:], [-0.000117062180693,
2365                                 7.94434448109e-005,
2366                                 -0.000151505429018])
2367
2368
2369
2370    def test_balance_deep_and_shallow(self):
2371        """Test that balanced limiters preserve conserved quantites.
2372        """
2373        import copy
2374
2375        a = [0.0, 0.0]
2376        b = [0.0, 2.0]
2377        c = [2.0, 0.0]
2378        d = [0.0, 4.0]
2379        e = [2.0, 2.0]
2380        f = [4.0, 0.0]
2381
2382        points = [a, b, c, d, e, f]
2383
2384        #bac, bce, ecf, dbe
2385        elements = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4] ]
2386
2387        mesh = Domain(points, elements)
2388        mesh.check_integrity()
2389
2390        #Create a deliberate overshoot
2391        mesh.set_quantity('stage', [[3,0,3], [2,2,6], [5,3,8], [8,3,5]])
2392        mesh.set_quantity('elevation', 0) #Flat bed
2393        stage = mesh.quantities['stage']
2394
2395        ref_centroid_values = copy.copy(stage.centroid_values[:]) #Copy
2396
2397        #Limit
2398        mesh.distribute_to_vertices_and_edges()
2399
2400        #Assert that quantities are conserved
2401        from Numeric import sum
2402        for k in range(len(mesh)):
2403            assert allclose (ref_centroid_values[k],
2404                             sum(stage.vertex_values[k,:])/3)
2405
2406
2407        #Now try with a non-flat bed - closely hugging initial stage in places
2408        #This will create alphas in the range [0, 0.478260, 1]
2409        mesh.set_quantity('stage', [[3,0,3], [2,2,6], [5,3,8], [8,3,5]])
2410        mesh.set_quantity('elevation', [[0,0,0],
2411                                        [1.8,1.9,5.9],
2412                                        [4.6,0,0],
2413                                        [0,2,4]])
2414        stage = mesh.quantities['stage']
2415
2416        ref_centroid_values = copy.copy(stage.centroid_values[:]) #Copy
2417        ref_vertex_values = copy.copy(stage.vertex_values[:]) #Copy
2418
2419        #Limit
2420        mesh.distribute_to_vertices_and_edges()
2421
2422
2423        #Assert that all vertex quantities have changed
2424        for k in range(len(mesh)):
2425            #print ref_vertex_values[k,:], stage.vertex_values[k,:]
2426            assert not allclose (ref_vertex_values[k,:], stage.vertex_values[k,:])
2427        #and assert that quantities are still conserved
2428        from Numeric import sum
2429        for k in range(len(mesh)):
2430            assert allclose (ref_centroid_values[k],
2431                             sum(stage.vertex_values[k,:])/3)
2432
2433
2434        #Also check that Python and C version produce the same
2435        # No longer applicable if tight_slope_limiters == 1
2436        #print stage.vertex_values
2437        #assert allclose (stage.vertex_values,
2438        #                 [[2,2,2],
2439        #                  [1.93333333, 2.03333333, 6.03333333],
2440        #                  [6.93333333, 4.53333333, 4.53333333],
2441        #                  [5.33333333, 5.33333333, 5.33333333]])
2442
2443
2444
2445
2446    def test_conservation_1(self):
2447        """Test that stage is conserved globally
2448
2449        This one uses a flat bed, reflective bdries and a suitable
2450        initial condition
2451        """
2452        from mesh_factory import rectangular
2453        from Numeric import array
2454
2455        #Create basic mesh
2456        points, vertices, boundary = rectangular(6, 6)
2457
2458        #Create shallow water domain
2459        domain = Domain(points, vertices, boundary)
2460        domain.smooth = False
2461        domain.default_order=2
2462
2463        #IC
2464        def x_slope(x, y):
2465            return x/3
2466
2467        domain.set_quantity('elevation', 0)
2468        domain.set_quantity('friction', 0)
2469        domain.set_quantity('stage', x_slope)
2470
2471        # Boundary conditions (reflective everywhere)
2472        Br = Reflective_boundary(domain)
2473        domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
2474
2475        domain.check_integrity()
2476
2477        initial_volume = domain.quantities['stage'].get_integral()
2478        initial_xmom = domain.quantities['xmomentum'].get_integral()
2479
2480        #print initial_xmom
2481
2482        #Evolution
2483        for t in domain.evolve(yieldstep = 0.05, finaltime = 5.0):
2484            volume =  domain.quantities['stage'].get_integral()
2485            assert allclose (volume, initial_volume)
2486
2487            #I don't believe that the total momentum should be the same
2488            #It starts with zero and ends with zero though
2489            #xmom = domain.quantities['xmomentum'].get_integral()
2490            #print xmom
2491            #assert allclose (xmom, initial_xmom)
2492
2493        os.remove(domain.get_name() + '.sww')
2494
2495
2496    def test_conservation_2(self):
2497        """Test that stage is conserved globally
2498
2499        This one uses a slopy bed, reflective bdries and a suitable
2500        initial condition
2501        """
2502        from mesh_factory import rectangular
2503        from Numeric import array
2504
2505        #Create basic mesh
2506        points, vertices, boundary = rectangular(6, 6)
2507
2508        #Create shallow water domain
2509        domain = Domain(points, vertices, boundary)
2510        domain.smooth = False
2511        domain.default_order=2
2512
2513        #IC
2514        def x_slope(x, y):
2515            return x/3
2516
2517        domain.set_quantity('elevation', x_slope)
2518        domain.set_quantity('friction', 0)
2519        domain.set_quantity('stage', 0.4) #Steady
2520
2521        # Boundary conditions (reflective everywhere)
2522        Br = Reflective_boundary(domain)
2523        domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
2524
2525        domain.check_integrity()
2526
2527        initial_volume = domain.quantities['stage'].get_integral()
2528        initial_xmom = domain.quantities['xmomentum'].get_integral()
2529
2530        #print initial_xmom
2531
2532        #Evolution
2533        for t in domain.evolve(yieldstep = 0.05, finaltime = 5.0):
2534            volume =  domain.quantities['stage'].get_integral()
2535            assert allclose (volume, initial_volume)
2536
2537            #FIXME: What would we expect from momentum
2538            #xmom = domain.quantities['xmomentum'].get_integral()
2539            #print xmom
2540            #assert allclose (xmom, initial_xmom)
2541
2542        os.remove(domain.get_name() + '.sww')
2543
2544    def test_conservation_3(self):
2545        """Test that stage is conserved globally
2546
2547        This one uses a larger grid, convoluted bed, reflective bdries and a suitable
2548        initial condition
2549        """
2550        from mesh_factory import rectangular
2551        from Numeric import array
2552
2553        #Create basic mesh
2554        points, vertices, boundary = rectangular(2, 1)
2555
2556        #Create shallow water domain
2557        domain = Domain(points, vertices, boundary)
2558        domain.smooth = False
2559        domain.default_order = 2
2560        domain.beta_h = 0.2
2561        domain.set_quantities_to_be_stored(['stage'])
2562
2563        #IC
2564        def x_slope(x, y):
2565            z = 0*x
2566            for i in range(len(x)):
2567                if x[i] < 0.3:
2568                    z[i] = x[i]/3
2569                if 0.3 <= x[i] < 0.5:
2570                    z[i] = -0.5
2571                if 0.5 <= x[i] < 0.7:
2572                    z[i] = 0.39
2573                if 0.7 <= x[i]:
2574                    z[i] = x[i]/3
2575            return z
2576
2577
2578
2579        domain.set_quantity('elevation', x_slope)
2580        domain.set_quantity('friction', 0)
2581        domain.set_quantity('stage', 0.4) #Steady
2582
2583        # Boundary conditions (reflective everywhere)
2584        Br = Reflective_boundary(domain)
2585        domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
2586
2587        domain.check_integrity()
2588
2589        initial_volume = domain.quantities['stage'].get_integral()
2590        initial_xmom = domain.quantities['xmomentum'].get_integral()
2591
2592        import copy
2593        ref_centroid_values =\
2594                 copy.copy(domain.quantities['stage'].centroid_values)
2595
2596        #print 'ORG', domain.quantities['stage'].centroid_values
2597        domain.distribute_to_vertices_and_edges()
2598
2599
2600        #print domain.quantities['stage'].centroid_values
2601        assert allclose(domain.quantities['stage'].centroid_values,
2602                        ref_centroid_values)
2603
2604
2605        #Check that initial limiter doesn't violate cons quan
2606        assert allclose (domain.quantities['stage'].get_integral(),
2607                         initial_volume)
2608
2609        #Evolution
2610        for t in domain.evolve(yieldstep = 0.05, finaltime = 10):
2611            volume =  domain.quantities['stage'].get_integral()
2612            #print t, volume, initial_volume
2613            assert allclose (volume, initial_volume)
2614
2615        os.remove(domain.get_name() + '.sww')
2616
2617    def test_conservation_4(self):
2618        """Test that stage is conserved globally
2619
2620        This one uses a larger grid, convoluted bed, reflective bdries and a suitable
2621        initial condition
2622        """
2623        from mesh_factory import rectangular
2624        from Numeric import array
2625
2626        #Create basic mesh
2627        points, vertices, boundary = rectangular(6, 6)
2628
2629        #Create shallow water domain
2630        domain = Domain(points, vertices, boundary)
2631        domain.smooth = False
2632        domain.default_order=2
2633        domain.beta_h = 0.0
2634
2635        #IC
2636        def x_slope(x, y):
2637            z = 0*x
2638            for i in range(len(x)):
2639                if x[i] < 0.3:
2640                    z[i] = x[i]/3
2641                if 0.3 <= x[i] < 0.5:
2642                    z[i] = -0.5
2643                if 0.5 <= x[i] < 0.7:
2644                    #z[i] = 0.3 #OK with beta == 0.2
2645                    z[i] = 0.34 #OK with beta == 0.0
2646                    #z[i] = 0.35#Fails after 80 timesteps with an error
2647                                #of the order 1.0e-5
2648                if 0.7 <= x[i]:
2649                    z[i] = x[i]/3
2650            return z
2651
2652
2653
2654        domain.set_quantity('elevation', x_slope)
2655        domain.set_quantity('friction', 0)
2656        domain.set_quantity('stage', 0.4) #Steady
2657
2658        # Boundary conditions (reflective everywhere)
2659        Br = Reflective_boundary(domain)
2660        domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
2661
2662        domain.check_integrity()
2663
2664        initial_volume = domain.quantities['stage'].get_integral()
2665        initial_xmom = domain.quantities['xmomentum'].get_integral()
2666
2667        import copy
2668        ref_centroid_values =\
2669                 copy.copy(domain.quantities['stage'].centroid_values)
2670
2671        #Test limiter by itself
2672        domain.distribute_to_vertices_and_edges()
2673
2674        #Check that initial limiter doesn't violate cons quan
2675        assert allclose (domain.quantities['stage'].get_integral(),
2676                         initial_volume)
2677        #NOTE: This would fail if any initial stage was less than the
2678        #corresponding bed elevation - but that is reasonable.
2679
2680
2681        #Evolution
2682        for t in domain.evolve(yieldstep = 0.05, finaltime = 10.0):
2683            volume =  domain.quantities['stage'].get_integral()
2684
2685            #print t, volume, initial_volume
2686
2687            assert allclose (volume, initial_volume)
2688
2689
2690        os.remove(domain.get_name() + '.sww')
2691
2692
2693    def test_conservation_5(self):
2694        """Test that momentum is conserved globally in
2695        steady state scenario
2696
2697        This one uses a slopy bed, dirichlet and reflective bdries
2698        """
2699        from mesh_factory import rectangular
2700        from Numeric import array
2701
2702        #Create basic mesh
2703        points, vertices, boundary = rectangular(6, 6)
2704
2705        #Create shallow water domain
2706        domain = Domain(points, vertices, boundary)
2707        domain.smooth = False
2708        domain.default_order=2
2709
2710        #IC
2711        def x_slope(x, y):
2712            return x/3
2713
2714        domain.set_quantity('elevation', x_slope)
2715        domain.set_quantity('friction', 0)
2716        domain.set_quantity('stage', 0.4) #Steady
2717
2718        # Boundary conditions (reflective everywhere)
2719        Br = Reflective_boundary(domain)
2720        Bleft = Dirichlet_boundary([0.5,0,0])
2721        Bright = Dirichlet_boundary([0.1,0,0])
2722        domain.set_boundary({'left': Bleft, 'right': Bright,
2723                             'top': Br, 'bottom': Br})
2724
2725        domain.check_integrity()
2726
2727        initial_volume = domain.quantities['stage'].get_integral()
2728        initial_xmom = domain.quantities['xmomentum'].get_integral()
2729
2730
2731        #Evolution
2732        for t in domain.evolve(yieldstep = 0.05, finaltime = 15.0):
2733            stage =  domain.quantities['stage'].get_integral()
2734            xmom = domain.quantities['xmomentum'].get_integral()
2735            ymom = domain.quantities['ymomentum'].get_integral()
2736
2737            if allclose(t, 6):  #Steady state reached
2738                steady_xmom = domain.quantities['xmomentum'].get_integral()
2739                steady_ymom = domain.quantities['ymomentum'].get_integral()
2740                steady_stage = domain.quantities['stage'].get_integral()
2741
2742            if t > 6:
2743                #print '%.2f %14.8f %14.8f' %(t, ymom, steady_ymom)
2744                assert allclose(xmom, steady_xmom)
2745                assert allclose(ymom, steady_ymom)
2746                assert allclose(stage, steady_stage)
2747
2748
2749        os.remove(domain.get_name() + '.sww')
2750
2751
2752
2753
2754
2755    def test_conservation_real(self):
2756        """Test that momentum is conserved globally
2757       
2758        Stephen finally made a test that revealed the problem.
2759        This test failed with code prior to 25 July 2005
2760        """
2761
2762        yieldstep = 0.01
2763        finaltime = 0.05
2764        min_depth = 1.0e-2
2765
2766
2767        import sys
2768        from os import sep; sys.path.append('..'+sep+'abstract_2d_finite_volumes')
2769        from mesh_factory import rectangular
2770
2771
2772        #Create shallow water domain
2773        points, vertices, boundary = rectangular(10, 10, len1=500, len2=500)
2774        domain = Domain(points, vertices, boundary)
2775        domain.smooth = False
2776        domain.default_order = 1
2777        domain.minimum_allowed_height = min_depth
2778
2779        # Set initial condition
2780        class Set_IC:
2781            """Set an initial condition with a constant value, for x0<x<x1
2782            """
2783
2784            def __init__(self, x0=0.25, x1=0.5, h=1.0):
2785                self.x0 = x0
2786                self.x1 = x1
2787                self.= h
2788
2789            def __call__(self, x, y):
2790                return self.h*((x>self.x0)&(x<self.x1))
2791
2792
2793        domain.set_quantity('stage', Set_IC(200.0,300.0,5.0))
2794
2795
2796        #Boundaries
2797        R = Reflective_boundary(domain)
2798        domain.set_boundary( {'left': R, 'right': R, 'top':R, 'bottom': R})
2799
2800        ref = domain.quantities['stage'].get_integral()
2801
2802        # Evolution
2803        for t in domain.evolve(yieldstep = yieldstep, finaltime = finaltime):
2804            pass
2805            #print 'Integral stage = ',\
2806            #      domain.quantities['stage'].get_integral(),\
2807            #      ' Time = ',domain.time
2808
2809
2810        now = domain.quantities['stage'].get_integral()
2811
2812        msg = 'Stage not conserved: was %f, now %f' %(ref, now)
2813        assert allclose(ref, now), msg
2814
2815        os.remove(domain.get_name() + '.sww')
2816
2817    def test_second_order_flat_bed_onestep(self):
2818
2819        from mesh_factory import rectangular
2820        from Numeric import array
2821
2822        #Create basic mesh
2823        points, vertices, boundary = rectangular(6, 6)
2824
2825        #Create shallow water domain
2826        domain = Domain(points, vertices, boundary)
2827        domain.smooth = False
2828        domain.default_order=2
2829        domain.beta_w      = 0.9
2830        domain.beta_w_dry  = 0.9
2831        domain.beta_uh     = 0.9
2832        domain.beta_uh_dry = 0.9
2833        domain.beta_vh     = 0.9
2834        domain.beta_vh_dry = 0.9
2835        domain.H0 = 0 # Backwards compatibility (6/2/7)       
2836       
2837        # Boundary conditions
2838        Br = Reflective_boundary(domain)
2839        Bd = Dirichlet_boundary([0.1, 0., 0.])
2840        domain.set_boundary({'left': Bd, 'right': Br, 'top': Br, 'bottom': Br})
2841
2842        domain.check_integrity()
2843
2844        #Evolution
2845        for t in domain.evolve(yieldstep = 0.05, finaltime = 0.05):
2846            pass# domain.write_time()
2847
2848        #Data from earlier version of abstract_2d_finite_volumes
2849        assert allclose(domain.min_timestep, 0.0396825396825)
2850        assert allclose(domain.max_timestep, 0.0396825396825)
2851
2852        assert allclose(domain.quantities['stage'].centroid_values[:12],
2853                        [0.00171396, 0.02656103, 0.00241523, 0.02656103,
2854                        0.00241523, 0.02656103, 0.00241523, 0.02656103,
2855                        0.00241523, 0.02656103, 0.00241523, 0.0272623])
2856
2857        domain.distribute_to_vertices_and_edges()
2858        assert allclose(domain.quantities['stage'].vertex_values[:12,0],
2859                        [0.0001714, 0.02656103, 0.00024152,
2860                        0.02656103, 0.00024152, 0.02656103,
2861                        0.00024152, 0.02656103, 0.00024152,
2862                        0.02656103, 0.00024152, 0.0272623])
2863
2864        assert allclose(domain.quantities['stage'].vertex_values[:12,1],
2865                        [0.00315012, 0.02656103, 0.00024152, 0.02656103,
2866                         0.00024152, 0.02656103, 0.00024152, 0.02656103,
2867                         0.00024152, 0.02656103, 0.00040506, 0.0272623])
2868
2869        assert allclose(domain.quantities['stage'].vertex_values[:12,2],
2870                        [0.00182037, 0.02656103, 0.00676264,
2871                         0.02656103, 0.00676264, 0.02656103,
2872                         0.00676264, 0.02656103, 0.00676264,
2873                         0.02656103, 0.0065991, 0.0272623])
2874
2875        assert allclose(domain.quantities['xmomentum'].centroid_values[:12],
2876                        [0.00113961, 0.01302432, 0.00148672,
2877                         0.01302432, 0.00148672, 0.01302432,
2878                         0.00148672, 0.01302432, 0.00148672 ,
2879                         0.01302432, 0.00148672, 0.01337143])
2880
2881        assert allclose(domain.quantities['ymomentum'].centroid_values[:12],
2882                        [-2.91240050e-004 , 1.22721531e-004,
2883                         -1.22721531e-004, 1.22721531e-004 ,
2884                         -1.22721531e-004, 1.22721531e-004,
2885                         -1.22721531e-004 , 1.22721531e-004,
2886                         -1.22721531e-004, 1.22721531e-004,
2887                         -1.22721531e-004, -4.57969873e-005])
2888
2889        os.remove(domain.get_name() + '.sww')
2890
2891
2892    def test_second_order_flat_bed_moresteps(self):
2893
2894        from mesh_factory import rectangular
2895        from Numeric import array
2896
2897        #Create basic mesh
2898        points, vertices, boundary = rectangular(6, 6)
2899
2900        #Create shallow water domain
2901        domain = Domain(points, vertices, boundary)
2902        domain.smooth = False
2903        domain.default_order=2
2904
2905        # Boundary conditions
2906        Br = Reflective_boundary(domain)
2907        Bd = Dirichlet_boundary([0.1, 0., 0.])
2908        domain.set_boundary({'left': Bd, 'right': Br, 'top': Br, 'bottom': Br})
2909
2910        domain.check_integrity()
2911
2912        #Evolution
2913        for t in domain.evolve(yieldstep = 0.05, finaltime = 0.1):
2914            pass
2915
2916        #Data from earlier version of abstract_2d_finite_volumes
2917        #assert allclose(domain.min_timestep, 0.0396825396825)
2918        #assert allclose(domain.max_timestep, 0.0396825396825)
2919        #print domain.quantities['stage'].centroid_values
2920
2921        os.remove(domain.get_name() + '.sww')       
2922
2923
2924    def test_flatbed_first_order(self):
2925        from mesh_factory import rectangular
2926        from Numeric import array
2927
2928        #Create basic mesh
2929        N = 8
2930        points, vertices, boundary = rectangular(N, N)
2931
2932        #Create shallow water domain
2933        domain = Domain(points, vertices, boundary)
2934        domain.smooth = False
2935        domain.default_order=1
2936        domain.H0 = 0 # Backwards compatibility (6/2/7)       
2937
2938        # Boundary conditions
2939        Br = Reflective_boundary(domain)
2940        Bd = Dirichlet_boundary([0.2,0.,0.])
2941
2942        domain.set_boundary({'left': Bd, 'right': Br, 'top': Br, 'bottom': Br})
2943        domain.check_integrity()
2944
2945
2946        #Evolution
2947        for t in domain.evolve(yieldstep = 0.02, finaltime = 0.5):
2948            pass
2949            #domain.write_time()
2950
2951        #FIXME: These numbers were from version before 25/10
2952        #assert allclose(domain.min_timestep, 0.0140413643926)
2953        #assert allclose(domain.max_timestep, 0.0140947355753)
2954
2955        for i in range(3):
2956            #assert allclose(domain.quantities['stage'].edge_values[:4,i],
2957            #                [0.10730244,0.12337617,0.11200126,0.12605666])
2958
2959            assert allclose(domain.quantities['xmomentum'].edge_values[:4,i],
2960                            [0.07610894,0.06901572,0.07284461,0.06819712])
2961
2962            #assert allclose(domain.quantities['ymomentum'].edge_values[:4,i],
2963            #                [-0.0060238, -0.00157404, -0.00309633, -0.0001637])
2964
2965
2966        os.remove(domain.get_name() + '.sww')           
2967
2968    def test_flatbed_second_order(self):
2969        from mesh_factory import rectangular
2970        from Numeric import array
2971
2972        #Create basic mesh
2973        N = 8
2974        points, vertices, boundary = rectangular(N, N)
2975
2976        #Create shallow water domain
2977        domain = Domain(points, vertices, boundary)
2978        domain.smooth = False
2979        domain.default_order=2
2980        domain.beta_w      = 0.9
2981        domain.beta_w_dry  = 0.9
2982        domain.beta_uh     = 0.9
2983        domain.beta_uh_dry = 0.9
2984        domain.beta_vh     = 0.9
2985        domain.beta_vh_dry = 0.9       
2986        #domain.minimum_allowed_height = 0.0 #Makes it like the 'oldstyle' balance
2987        domain.H0 = 0 # Backwards compatibility (6/2/7)       
2988
2989        # Boundary conditions
2990        Br = Reflective_boundary(domain)
2991        Bd = Dirichlet_boundary([0.2,0.,0.])
2992
2993        domain.set_boundary({'left': Bd, 'right': Br, 'top': Br, 'bottom': Br})
2994        domain.check_integrity()
2995
2996        #Evolution
2997        for t in domain.evolve(yieldstep = 0.01, finaltime = 0.03):
2998            pass
2999
3000
3001        assert allclose(domain.min_timestep, 0.0210448446782)
3002        assert allclose(domain.max_timestep, 0.0210448446782)
3003
3004        #print domain.quantities['stage'].vertex_values[:4,0]
3005        #print domain.quantities['xmomentum'].vertex_values[:4,0]
3006        #print domain.quantities['ymomentum'].vertex_values[:4,0]
3007
3008        #FIXME: These numbers were from version before 25/10
3009        #assert allclose(domain.quantities['stage'].vertex_values[:4,0],
3010        #                [0.00101913,0.05352143,0.00104852,0.05354394])
3011
3012        #FIXME: These numbers were from version before 21/3/6 -
3013        #could be recreated by setting maximum_allowed_speed to 0 maybe
3014        #assert allclose(domain.quantities['xmomentum'].vertex_values[:4,0],
3015        #                [ 0.00064835, 0.03685719, 0.00085073, 0.03687313])
3016       
3017        assert allclose(domain.quantities['xmomentum'].vertex_values[:4,0],
3018                        [ 0.00090581, 0.03685719, 0.00088303, 0.03687313])
3019                       
3020                       
3021
3022        #assert allclose(domain.quantities['xmomentum'].vertex_values[:4,0],
3023        #                [0.00090581,0.03685719,0.00088303,0.03687313])
3024
3025        assert allclose(domain.quantities['ymomentum'].vertex_values[:4,0],
3026                        [-0.00139463,0.0006156,-0.00060364,0.00061827])
3027
3028
3029        os.remove(domain.get_name() + '.sww')
3030
3031       
3032    def test_flatbed_second_order_vmax_0(self):
3033        from mesh_factory import rectangular
3034        from Numeric import array
3035
3036        #Create basic mesh
3037        N = 8
3038        points, vertices, boundary = rectangular(N, N)
3039
3040        #Create shallow water domain
3041        domain = Domain(points, vertices, boundary)
3042        domain.smooth = False
3043        domain.default_order=2
3044        domain.beta_w      = 0.9
3045        domain.beta_w_dry  = 0.9
3046        domain.beta_uh     = 0.9
3047        domain.beta_uh_dry = 0.9
3048        domain.beta_vh     = 0.9
3049        domain.beta_vh_dry = 0.9       
3050        domain.maximum_allowed_speed = 0.0 #Makes it like the 'oldstyle'
3051        domain.H0 = 0 # Backwards compatibility (6/2/7)       
3052
3053        # Boundary conditions
3054        Br = Reflective_boundary(domain)
3055        Bd = Dirichlet_boundary([0.2,0.,0.])
3056
3057        domain.set_boundary({'left': Bd, 'right': Br, 'top': Br, 'bottom': Br})
3058        domain.check_integrity()
3059
3060        #Evolution
3061        for t in domain.evolve(yieldstep = 0.01, finaltime = 0.03):
3062            pass
3063
3064
3065        assert allclose(domain.min_timestep, 0.0210448446782)
3066        assert allclose(domain.max_timestep, 0.0210448446782)
3067
3068        #FIXME: These numbers were from version before 21/3/6 -
3069        #could be recreated by setting maximum_allowed_speed to 0 maybe
3070        assert allclose(domain.quantities['xmomentum'].vertex_values[:4,0],
3071                        [ 0.00064835, 0.03685719, 0.00085073, 0.03687313]) 
3072       
3073
3074        assert allclose(domain.quantities['ymomentum'].vertex_values[:4,0],
3075                        [-0.00139463,0.0006156,-0.00060364,0.00061827])
3076
3077
3078        os.remove(domain.get_name() + '.sww')
3079
3080       
3081
3082    def test_flatbed_second_order_distribute(self):
3083        #Use real data from anuga.abstract_2d_finite_volumes 2
3084        #painfully setup and extracted.
3085        from mesh_factory import rectangular
3086        from Numeric import array
3087
3088        #Create basic mesh
3089        N = 8
3090        points, vertices, boundary = rectangular(N, N)
3091
3092        #Create shallow water domain
3093        domain = Domain(points, vertices, boundary)
3094        domain.smooth = False
3095        domain.default_order=domain._order_=2
3096        domain.beta_w      = 0.9
3097        domain.beta_w_dry  = 0.9
3098        domain.beta_uh     = 0.9
3099        domain.beta_uh_dry = 0.9
3100        domain.beta_vh     = 0.9
3101        domain.beta_vh_dry = 0.9
3102        domain.H0 = 0 # Backwards compatibility (6/2/7)       
3103
3104        # Boundary conditions
3105        Br = Reflective_boundary(domain)
3106        Bd = Dirichlet_boundary([0.2,0.,0.])
3107
3108        domain.set_boundary({'left': Bd, 'right': Br, 'top': Br, 'bottom': Br})
3109        domain.check_integrity()
3110
3111
3112
3113        for V in [False, True]:
3114            if V:
3115                #Set centroids as if system had been evolved
3116                L = zeros(2*N*N, Float)
3117                L[:32] = [7.21205592e-003, 5.35214298e-002, 1.00910824e-002,
3118                          5.35439433e-002, 1.00910824e-002, 5.35439433e-002,
3119                          1.00910824e-002, 5.35439433e-002, 1.00910824e-002,
3120                          5.35439433e-002, 1.00910824e-002, 5.35439433e-002,
3121                          1.00910824e-002, 5.35393928e-002, 1.02344264e-002,
3122                          5.59605058e-002, 0.00000000e+000, 3.31027800e-004,
3123                          0.00000000e+000, 4.37962142e-005, 0.00000000e+000,
3124                          4.37962142e-005, 0.00000000e+000, 4.37962142e-005,
3125                          0.00000000e+000, 4.37962142e-005, 0.00000000e+000,
3126                          4.37962142e-005, 0.00000000e+000, 4.37962142e-005,
3127                          0.00000000e+000, 5.57305948e-005]
3128
3129                X = zeros(2*N*N, Float)
3130                X[:32] = [6.48351607e-003, 3.68571894e-002, 8.50733285e-003,
3131                          3.68731327e-002, 8.50733285e-003, 3.68731327e-002,
3132                          8.50733285e-003, 3.68731327e-002, 8.50733285e-003,
3133                          3.68731327e-002, 8.50733285e-003, 3.68731327e-002,
3134                          8.50733285e-003, 3.68693861e-002, 8.65220973e-003,
3135                          3.85055387e-002, 0.00000000e+000, 2.86060840e-004,
3136                          0.00000000e+000, 3.58905503e-005, 0.00000000e+000,
3137                          3.58905503e-005, 0.00000000e+000, 3.58905503e-005,
3138                          0.00000000e+000, 3.58905503e-005, 0.00000000e+000,
3139                          3.58905503e-005, 0.00000000e+000, 3.58905503e-005,
3140                          0.00000000e+000, 4.57662812e-005]
3141
3142                Y = zeros(2*N*N, Float)
3143                Y[:32]=[-1.39463104e-003, 6.15600298e-004, -6.03637382e-004,
3144                        6.18272251e-004, -6.03637382e-004, 6.18272251e-004,
3145                        -6.03637382e-004, 6.18272251e-004, -6.03637382e-004,
3146                        6.18272251e-004, -6.03637382e-004, 6.18272251e-004,
3147                        -6.03637382e-004, 6.18599320e-004, -6.74622797e-004,
3148                        -1.48934756e-004, 0.00000000e+000, -5.35079969e-005,
3149                        0.00000000e+000, -2.57264987e-005, 0.00000000e+000,
3150                        -2.57264987e-005, 0.00000000e+000, -2.57264987e-005,
3151                        0.00000000e+000, -2.57264987e-005, 0.00000000e+000,
3152                        -2.57264987e-005, 0.00000000e+000, -2.57264987e-005,
3153                        0.00000000e+000, -2.57635178e-005]
3154
3155
3156                domain.set_quantity('stage', L, location='centroids')
3157                domain.set_quantity('xmomentum', X, location='centroids')
3158                domain.set_quantity('ymomentum', Y, location='centroids')
3159
3160                domain.check_integrity()
3161            else:
3162                #Evolution
3163                for t in domain.evolve(yieldstep = 0.01, finaltime = 0.03):
3164                    pass
3165                assert allclose(domain.min_timestep, 0.0210448446782)
3166                assert allclose(domain.max_timestep, 0.0210448446782)
3167
3168
3169            #Centroids were correct but not vertices.
3170            #Hence the check of distribute below.
3171            assert allclose(domain.quantities['stage'].centroid_values[:4],
3172                            [0.00721206,0.05352143,0.01009108,0.05354394])
3173
3174            assert allclose(domain.quantities['xmomentum'].centroid_values[:4],
3175                            [0.00648352,0.03685719,0.00850733,0.03687313])
3176
3177            assert allclose(domain.quantities['ymomentum'].centroid_values[:4],
3178                            [-0.00139463,0.0006156,-0.00060364,0.00061827])
3179
3180            #print 'C17=', domain.quantities['xmomentum'].centroid_values[17]
3181            #print 'C19=', domain.quantities['xmomentum'].centroid_values[19]
3182
3183            #assert allclose(domain.quantities['xmomentum'].centroid_values[17],0.00028606084)
3184            ##print domain.quantities['xmomentum'].centroid_values[17], V
3185            ##print
3186            if not V:
3187                #FIXME: These numbers were from version before 21/3/6 -
3188                #could be recreated by setting maximum_allowed_speed to 0 maybe           
3189               
3190                #assert allclose(domain.quantities['xmomentum'].centroid_values[17], 0.0)               
3191                assert allclose(domain.quantities['xmomentum'].centroid_values[17], 0.000286060839592)                         
3192
3193            else:
3194                assert allclose(domain.quantities['xmomentum'].centroid_values[17], 0.00028606084)
3195
3196            import copy
3197            XX = copy.copy(domain.quantities['xmomentum'].centroid_values)
3198            assert allclose(domain.quantities['xmomentum'].centroid_values, XX)
3199
3200            domain.distribute_to_vertices_and_edges()
3201
3202            #assert allclose(domain.quantities['xmomentum'].centroid_values, XX)
3203
3204            #assert allclose(domain.quantities['xmomentum'].centroid_values[17],
3205            #                0.0)
3206            assert allclose(domain.quantities['xmomentum'].centroid_values[17], 0.000286060839592)                                                 
3207
3208
3209            #FIXME: These numbers were from version before 25/10
3210            #assert allclose(domain.quantities['stage'].vertex_values[:4,0],
3211            #                [0.00101913,0.05352143,0.00104852,0.05354394])
3212
3213            assert allclose(domain.quantities['ymomentum'].vertex_values[:4,0],
3214                            [-0.00139463,0.0006156,-0.00060364,0.00061827])
3215
3216
3217            assert allclose(domain.quantities['xmomentum'].vertex_values[:4,0],
3218                            [0.00090581,0.03685719,0.00088303,0.03687313])
3219
3220
3221            #NB NO longer relvant:
3222
3223            #This was the culprit. First triangles vertex 0 had an
3224            #x-momentum of 0.0064835 instead of 0.00090581 and
3225            #third triangle had 0.00850733 instead of 0.00088303
3226            #print domain.quantities['xmomentum'].vertex_values[:4,0]
3227
3228            #print domain.quantities['xmomentum'].vertex_values[:4,0]
3229            #assert allclose(domain.quantities['xmomentum'].vertex_values[:4,0],
3230            #                [0.00090581,0.03685719,0.00088303,0.03687313])
3231
3232        os.remove(domain.get_name() + '.sww')
3233
3234
3235
3236    def test_bedslope_problem_first_order(self):
3237
3238        from mesh_factory import rectangular
3239        from Numeric import array
3240
3241        #Create basic mesh
3242        points, vertices, boundary = rectangular(6, 6)
3243
3244        #Create shallow water domain
3245        domain = Domain(points, vertices, boundary)
3246        domain.smooth = False
3247        domain.default_order = 1
3248
3249        #Bed-slope and friction
3250        def x_slope(x, y):
3251            return -x/3
3252
3253        domain.set_quantity('elevation', x_slope)
3254
3255        # Boundary conditions
3256        Br = Reflective_boundary(domain)
3257        domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
3258
3259        #Initial condition
3260        domain.set_quantity('stage', Constant_height(x_slope, 0.05))
3261        domain.check_integrity()
3262
3263        #Evolution
3264        for t in domain.evolve(yieldstep = 0.05, finaltime = 0.05):
3265            pass# domain.write_time()
3266
3267        # FIXME (Ole): Need some other assertion here!
3268        #print domain.min_timestep, domain.max_timestep   
3269        #assert allclose(domain.min_timestep, 0.050010003001)
3270        #assert allclose(domain.max_timestep, 0.050010003001)
3271
3272
3273        os.remove(domain.get_name() + '.sww')
3274       
3275    def test_bedslope_problem_first_order_moresteps(self):
3276
3277        from mesh_factory import rectangular
3278        from Numeric import array
3279
3280        #Create basic mesh
3281        points, vertices, boundary = rectangular(6, 6)
3282
3283        #Create shallow water domain
3284        domain = Domain(points, vertices, boundary)
3285        domain.smooth = False
3286        domain.default_order = 1
3287        domain.beta_h = 0.0 # Use first order in h-limiter
3288       
3289        # FIXME (Ole): Need tests where these two are commented out
3290        domain.H0 = 0        # Backwards compatibility (6/2/7)       
3291        domain.tight_slope_limiters = 0 # Backwards compatibility (14/4/7)         
3292
3293        #Bed-slope and friction
3294        def x_slope(x, y):
3295            return -x/3
3296
3297        domain.set_quantity('elevation', x_slope)
3298
3299        # Boundary conditions
3300        Br = Reflective_boundary(domain)
3301        domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
3302
3303        #Initial condition
3304        domain.set_quantity('stage', Constant_height(x_slope, 0.05))
3305        domain.check_integrity()
3306
3307        #Evolution
3308        for t in domain.evolve(yieldstep = 0.05, finaltime = 0.5):
3309            pass# domain.write_time()
3310
3311        #Data from earlier version of abstract_2d_finite_volumes
3312        #print domain.quantities['stage'].centroid_values
3313
3314        assert allclose(domain.quantities['stage'].centroid_values,
3315                        [-0.02998628, -0.01520652, -0.03043492,
3316                        -0.0149132, -0.03004706, -0.01476251,
3317                        -0.0298215, -0.01467976, -0.02988158,
3318                        -0.01474662, -0.03036161, -0.01442995,
3319                        -0.07624583, -0.06297061, -0.07733792,
3320                        -0.06342237, -0.07695439, -0.06289595,
3321                        -0.07635559, -0.0626065, -0.07633628,
3322                        -0.06280072, -0.07739632, -0.06386738,
3323                        -0.12161738, -0.11028239, -0.1223796,
3324                        -0.11095953, -0.12189744, -0.11048616,
3325                        -0.12074535, -0.10987605, -0.12014311,
3326                        -0.10976691, -0.12096859, -0.11087692,
3327                        -0.16868259, -0.15868061, -0.16801135,
3328                        -0.1588003, -0.16674343, -0.15813323,
3329                        -0.16457595, -0.15693826, -0.16281096,
3330                        -0.15585154, -0.16283873, -0.15540068,
3331                        -0.17450362, -0.19919913, -0.18062882,
3332                        -0.19764131, -0.17783111, -0.19407213,
3333                        -0.1736915, -0.19053624, -0.17228678,
3334                        -0.19105634, -0.17920133, -0.1968828,
3335                        -0.14244395, -0.14604641, -0.14473537,
3336                        -0.1506107, -0.14510055, -0.14919522,
3337                        -0.14175896, -0.14560798, -0.13911658,
3338                        -0.14439383, -0.13924047, -0.14829043])
3339
3340        os.remove(domain.get_name() + '.sww')
3341       
3342    def test_bedslope_problem_second_order_one_step(self):
3343
3344        from mesh_factory import rectangular
3345        from Numeric import array
3346
3347        #Create basic mesh
3348        points, vertices, boundary = rectangular(6, 6)
3349
3350        #Create shallow water domain
3351        domain = Domain(points, vertices, boundary)
3352        domain.smooth = False
3353        domain.default_order=2
3354        domain.beta_w      = 0.9
3355        domain.beta_w_dry  = 0.9
3356        domain.beta_uh     = 0.9
3357        domain.beta_uh_dry = 0.9
3358        domain.beta_vh     = 0.9
3359        domain.beta_vh_dry = 0.9
3360
3361       
3362        # FIXME (Ole): Need tests where this is commented out
3363        domain.tight_slope_limiters = 0 # Backwards compatibility (14/4/7)         
3364       
3365        #Bed-slope and friction at vertices (and interpolated elsewhere)
3366        def x_slope(x, y):
3367            return -x/3
3368
3369        domain.set_quantity('elevation', x_slope)
3370
3371        # Boundary conditions
3372        Br = Reflective_boundary(domain)
3373        domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
3374
3375        #Initial condition
3376        domain.set_quantity('stage', Constant_height(x_slope, 0.05))
3377        domain.check_integrity()
3378
3379        assert allclose(domain.quantities['stage'].centroid_values,
3380                        [0.01296296, 0.03148148, 0.01296296,
3381                        0.03148148, 0.01296296, 0.03148148,
3382                        0.01296296, 0.03148148, 0.01296296,
3383                        0.03148148, 0.01296296, 0.03148148,
3384                        -0.04259259, -0.02407407, -0.04259259,
3385                        -0.02407407, -0.04259259, -0.02407407,
3386                        -0.04259259, -0.02407407, -0.04259259,
3387                        -0.02407407, -0.04259259, -0.02407407,
3388                        -0.09814815, -0.07962963, -0.09814815,
3389                        -0.07962963, -0.09814815, -0.07962963,
3390                        -0.09814815, -0.07962963, -0.09814815,
3391                        -0.07962963, -0.09814815, -0.07962963,
3392                        -0.1537037 , -0.13518519, -0.1537037,
3393                        -0.13518519, -0.1537037, -0.13518519,
3394                        -0.1537037 , -0.13518519, -0.1537037,
3395                        -0.13518519, -0.1537037, -0.13518519,
3396                        -0.20925926, -0.19074074, -0.20925926,
3397                        -0.19074074, -0.20925926, -0.19074074,
3398                        -0.20925926, -0.19074074, -0.20925926,
3399                        -0.19074074, -0.20925926, -0.19074074,
3400                        -0.26481481, -0.2462963, -0.26481481,
3401                        -0.2462963, -0.26481481, -0.2462963,
3402                        -0.26481481, -0.2462963, -0.26481481,
3403                        -0.2462963, -0.26481481, -0.2462963])
3404
3405
3406        #print domain.quantities['stage'].extrapolate_second_order()
3407        #domain.distribute_to_vertices_and_edges()
3408        #print domain.quantities['stage'].vertex_values[:,0]
3409
3410        #Evolution
3411        for t in domain.evolve(yieldstep = 0.05, finaltime = 0.05):
3412            #domain.write_time()
3413            pass
3414
3415
3416        #print domain.quantities['stage'].centroid_values
3417        assert allclose(domain.quantities['stage'].centroid_values,
3418                 [0.01290985, 0.02356019, 0.01619096, 0.02356019, 0.01619096,
3419                  0.02356019, 0.01619096, 0.02356019, 0.01619096, 0.02356019,
3420                  0.01619096, 0.0268413, -0.04411074, -0.0248011, -0.04186556,
3421                  -0.0248011, -0.04186556, -0.0248011, -0.04186556, -0.0248011,
3422                  -0.04186556, -0.0248011, -0.04186556, -0.02255593,
3423                  -0.09966629, -0.08035666, -0.09742112, -0.08035666,
3424                  -0.09742112, -0.08035666, -0.09742112, -0.08035666,
3425                  -0.09742112, -0.08035666, -0.09742112, -0.07811149,
3426                  -0.15522185, -0.13591222, -0.15297667, -0.13591222,
3427                  -0.15297667, -0.13591222, -0.15297667, -0.13591222,
3428                  -0.15297667, -0.13591222, -0.15297667, -0.13366704,
3429                  -0.2107774, -0.19146777, -0.20853223, -0.19146777,
3430                  -0.20853223, -0.19146777, -0.20853223, -0.19146777,
3431                  -0.20853223, -0.19146777, -0.20853223, -0.1892226,
3432                  -0.26120669, -0.24776246, -0.25865535, -0.24776246,
3433                  -0.25865535, -0.24776246, -0.25865535, -0.24776246,
3434                  -0.25865535, -0.24776246, -0.25865535, -0.24521113])
3435
3436        os.remove(domain.get_name() + '.sww')
3437
3438    def test_bedslope_problem_second_order_two_steps(self):
3439
3440        from mesh_factory import rectangular
3441        from Numeric import array
3442
3443        #Create basic mesh
3444        points, vertices, boundary = rectangular(6, 6)
3445
3446        #Create shallow water domain
3447        domain = Domain(points, vertices, boundary)
3448        domain.smooth = False
3449        domain.default_order=2
3450        domain.beta_w      = 0.9
3451        domain.beta_w_dry  = 0.9
3452        domain.beta_uh     = 0.9
3453        domain.beta_uh_dry = 0.9
3454        domain.beta_vh     = 0.9
3455        domain.beta_vh_dry = 0.9
3456        domain.beta_h = 0.0 #Use first order in h-limiter
3457       
3458        # FIXME (Ole): Need tests where this is commented out
3459        domain.tight_slope_limiters = 0 # Backwards compatibility (14/4/7)                 
3460        domain.H0 = 0 # Backwards compatibility (6/2/7)       
3461
3462        #Bed-slope and friction at vertices (and interpolated elsewhere)
3463        def x_slope(x, y):
3464            return -x/3
3465
3466        domain.set_quantity('elevation', x_slope)
3467
3468        # Boundary conditions
3469        Br = Reflective_boundary(domain)
3470        domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
3471
3472        #Initial condition
3473        domain.set_quantity('stage', Constant_height(x_slope, 0.05))
3474        domain.check_integrity()
3475
3476        assert allclose(domain.quantities['stage'].centroid_values,
3477                        [0.01296296, 0.03148148, 0.01296296,
3478                        0.03148148, 0.01296296, 0.03148148,
3479                        0.01296296, 0.03148148, 0.01296296,
3480                        0.03148148, 0.01296296, 0.03148148,
3481                        -0.04259259, -0.02407407, -0.04259259,
3482                        -0.02407407, -0.04259259, -0.02407407,
3483                        -0.04259259, -0.02407407, -0.04259259,
3484                        -0.02407407, -0.04259259, -0.02407407,
3485                        -0.09814815, -0.07962963, -0.09814815,
3486                        -0.07962963, -0.09814815, -0.07962963,
3487                        -0.09814815, -0.07962963, -0.09814815,
3488                        -0.07962963, -0.09814815, -0.07962963,
3489                        -0.1537037 , -0.13518519, -0.1537037,
3490                        -0.13518519, -0.1537037, -0.13518519,
3491                        -0.1537037 , -0.13518519, -0.1537037,
3492                        -0.13518519, -0.1537037, -0.13518519,
3493                        -0.20925926, -0.19074074, -0.20925926,
3494                        -0.19074074, -0.20925926, -0.19074074,
3495                        -0.20925926, -0.19074074, -0.20925926,
3496                        -0.19074074, -0.20925926, -0.19074074,
3497                        -0.26481481, -0.2462963, -0.26481481,
3498                        -0.2462963, -0.26481481, -0.2462963,
3499                        -0.26481481, -0.2462963, -0.26481481,
3500                        -0.2462963, -0.26481481, -0.2462963])
3501
3502
3503        #print domain.quantities['stage'].extrapolate_second_order()
3504        #domain.distribute_to_vertices_and_edges()
3505        #print domain.quantities['stage'].vertex_values[:,0]
3506
3507        #Evolution
3508        for t in domain.evolve(yieldstep = 0.05, finaltime = 0.1):
3509            pass
3510
3511
3512        #Data from earlier version of abstract_2d_finite_volumes ft=0.1
3513        assert allclose(domain.min_timestep, 0.0376895634803)
3514        assert allclose(domain.max_timestep, 0.0415635655309)
3515
3516
3517        assert allclose(domain.quantities['stage'].centroid_values,
3518                        [0.00855788, 0.01575204, 0.00994606, 0.01717072,
3519                         0.01005985, 0.01716362, 0.01005985, 0.01716299,
3520                         0.01007098, 0.01736248, 0.01216452, 0.02026776,
3521                         -0.04462374, -0.02479045, -0.04199789, -0.0229465,
3522                         -0.04184033, -0.02295693, -0.04184013, -0.02295675,
3523                         -0.04184486, -0.0228168, -0.04028876, -0.02036486,
3524                         -0.10029444, -0.08170809, -0.09772846, -0.08021704,
3525                         -0.09760006, -0.08022143, -0.09759984, -0.08022124,
3526                         -0.09760261, -0.08008893, -0.09603914, -0.07758209,
3527                         -0.15584152, -0.13723138, -0.15327266, -0.13572906,
3528                         -0.15314427, -0.13573349, -0.15314405, -0.13573331,
3529                         -0.15314679, -0.13560104, -0.15158523, -0.13310701,
3530                         -0.21208605, -0.19283913, -0.20955631, -0.19134189,
3531                         -0.20942821, -0.19134598, -0.20942799, -0.1913458,
3532                         -0.20943005, -0.19120952, -0.20781177, -0.18869401,
3533                         -0.25384082, -0.2463294, -0.25047649, -0.24464654,
3534                         -0.25031159, -0.24464253, -0.25031112, -0.24464253,
3535                         -0.25031463, -0.24454764, -0.24885323, -0.24286438])
3536
3537
3538        os.remove(domain.get_name() + '.sww')
3539
3540    def test_bedslope_problem_second_order_two_yieldsteps(self):
3541
3542        from mesh_factory import rectangular
3543        from Numeric import array
3544
3545        #Create basic mesh
3546        points, vertices, boundary = rectangular(6, 6)
3547
3548        #Create shallow water domain
3549        domain = Domain(points, vertices, boundary)
3550        domain.smooth = False
3551        domain.default_order=2
3552        domain.beta_w      = 0.9
3553        domain.beta_w_dry  = 0.9
3554        domain.beta_uh     = 0.9
3555        domain.beta_uh_dry = 0.9
3556        domain.beta_vh     = 0.9
3557        domain.beta_vh_dry = 0.9
3558        domain.beta_h = 0.0 #Use first order in h-limiter
3559       
3560        # FIXME (Ole): Need tests where this is commented out
3561        domain.tight_slope_limiters = 0 # Backwards compatibility (14/4/7)                 
3562        domain.H0 = 0 # Backwards compatibility (6/2/7)
3563
3564        #Bed-slope and friction at vertices (and interpolated elsewhere)
3565        def x_slope(x, y):
3566            return -x/3
3567
3568        domain.set_quantity('elevation', x_slope)
3569
3570        # Boundary conditions
3571        Br = Reflective_boundary(domain)
3572        domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
3573
3574        #Initial condition
3575        domain.set_quantity('stage', Constant_height(x_slope, 0.05))
3576        domain.check_integrity()
3577
3578        assert allclose(domain.quantities['stage'].centroid_values,
3579                        [0.01296296, 0.03148148, 0.01296296,
3580                        0.03148148, 0.01296296, 0.03148148,
3581                        0.01296296, 0.03148148, 0.01296296,
3582                        0.03148148, 0.01296296, 0.03148148,
3583                        -0.04259259, -0.02407407, -0.04259259,
3584                        -0.02407407, -0.04259259, -0.02407407,
3585                        -0.04259259, -0.02407407, -0.04259259,
3586                        -0.02407407, -0.04259259, -0.02407407,
3587                        -0.09814815, -0.07962963, -0.09814815,
3588                        -0.07962963, -0.09814815, -0.07962963,
3589                        -0.09814815, -0.07962963, -0.09814815,
3590                        -0.07962963, -0.09814815, -0.07962963,
3591                        -0.1537037 , -0.13518519, -0.1537037,
3592                        -0.13518519, -0.1537037, -0.13518519,
3593                        -0.1537037 , -0.13518519, -0.1537037,
3594                        -0.13518519, -0.1537037, -0.13518519,
3595                        -0.20925926, -0.19074074, -0.20925926,
3596                        -0.19074074, -0.20925926, -0.19074074,
3597                        -0.20925926, -0.19074074, -0.20925926,
3598                        -0.19074074, -0.20925926, -0.19074074,
3599                        -0.26481481, -0.2462963, -0.26481481,
3600                        -0.2462963, -0.26481481, -0.2462963,
3601                        -0.26481481, -0.2462963, -0.26481481,
3602                        -0.2462963, -0.26481481, -0.2462963])
3603
3604
3605        #print domain.quantities['stage'].extrapolate_second_order()
3606        #domain.distribute_to_vertices_and_edges()
3607        #print domain.quantities['stage'].vertex_values[:,0]
3608
3609        #Evolution
3610        for t in domain.evolve(yieldstep = 0.05, finaltime = 0.1):   #0.05??
3611            #domain.write_time()
3612            pass
3613
3614
3615
3616        assert allclose(domain.quantities['stage'].centroid_values,
3617                 [0.00855788, 0.01575204, 0.00994606, 0.01717072, 0.01005985,
3618                  0.01716362, 0.01005985, 0.01716299, 0.01007098, 0.01736248,
3619                  0.01216452, 0.02026776, -0.04462374, -0.02479045, -0.04199789,
3620                  -0.0229465, -0.04184033, -0.02295693, -0.04184013,
3621                  -0.02295675, -0.04184486, -0.0228168, -0.04028876,
3622                  -0.02036486, -0.10029444, -0.08170809, -0.09772846,
3623                  -0.08021704, -0.09760006, -0.08022143, -0.09759984,
3624                  -0.08022124, -0.09760261, -0.08008893, -0.09603914,
3625                  -0.07758209, -0.15584152, -0.13723138, -0.15327266,
3626                  -0.13572906, -0.15314427, -0.13573349, -0.15314405,
3627                  -0.13573331, -0.15314679, -0.13560104, -0.15158523,
3628                  -0.13310701, -0.21208605, -0.19283913, -0.20955631,
3629                  -0.19134189, -0.20942821, -0.19134598, -0.20942799,
3630                  -0.1913458, -0.20943005, -0.19120952, -0.20781177,
3631                  -0.18869401, -0.25384082, -0.2463294, -0.25047649,
3632                  -0.24464654, -0.25031159, -0.24464253, -0.25031112,
3633                  -0.24464253, -0.25031463, -0.24454764, -0.24885323,
3634                  -0.24286438])
3635
3636        os.remove(domain.get_name() + '.sww')
3637
3638    def test_bedslope_problem_second_order_more_steps(self):
3639
3640        from mesh_factory import rectangular
3641        from Numeric import array
3642
3643        #Create basic mesh
3644        points, vertices, boundary = rectangular(6, 6)
3645
3646        #Create shallow water domain
3647        domain = Domain(points, vertices, boundary)
3648        domain.smooth = False
3649        domain.default_order=2
3650        domain.beta_w      = 0.9
3651        domain.beta_w_dry  = 0.9
3652        domain.beta_uh     = 0.9
3653        domain.beta_uh_dry = 0.9
3654        domain.beta_vh     = 0.9
3655        domain.beta_vh_dry = 0.9
3656        domain.beta_h = 0.0 #Use first order in h-limiter
3657       
3658       
3659        # FIXME (Ole): Need tests where these two are commented out
3660        domain.H0 = 0        # Backwards compatibility (6/2/7)       
3661        domain.tight_slope_limiters = 0 # Backwards compatibility (14/4/7)                 
3662       
3663               
3664
3665        #Bed-slope and friction at vertices (and interpolated elsewhere)
3666        def x_slope(x, y):
3667            return -x/3
3668
3669        domain.set_quantity('elevation', x_slope)
3670
3671        # Boundary conditions
3672        Br = Reflective_boundary(domain)
3673        domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
3674
3675        #Initial condition
3676        domain.set_quantity('stage', expression = 'elevation + 0.05')
3677        domain.check_integrity()
3678
3679        assert allclose(domain.quantities['stage'].centroid_values,
3680                        [0.01296296, 0.03148148, 0.01296296,
3681                        0.03148148, 0.01296296, 0.03148148,
3682                        0.01296296, 0.03148148, 0.01296296,
3683                        0.03148148, 0.01296296, 0.03148148,
3684                        -0.04259259, -0.02407407, -0.04259259,
3685                        -0.02407407, -0.04259259, -0.02407407,
3686                        -0.04259259, -0.02407407, -0.04259259,
3687                        -0.02407407, -0.04259259, -0.02407407,
3688                        -0.09814815, -0.07962963, -0.09814815,
3689                        -0.07962963, -0.09814815, -0.07962963,
3690                        -0.09814815, -0.07962963, -0.09814815,
3691                        -0.07962963, -0.09814815, -0.07962963,
3692                        -0.1537037 , -0.13518519, -0.1537037,
3693                        -0.13518519, -0.1537037, -0.13518519,
3694                        -0.1537037 , -0.13518519, -0.1537037,
3695                        -0.13518519, -0.1537037, -0.13518519,
3696                        -0.20925926, -0.19074074, -0.20925926,
3697                        -0.19074074, -0.20925926, -0.19074074,
3698                        -0.20925926, -0.19074074, -0.20925926,
3699                        -0.19074074, -0.20925926, -0.19074074,
3700                        -0.26481481, -0.2462963, -0.26481481,
3701                        -0.2462963, -0.26481481, -0.2462963,
3702                        -0.26481481, -0.2462963, -0.26481481,
3703                        -0.2462963, -0.26481481, -0.2462963])
3704
3705
3706        #print domain.quantities['stage'].extrapolate_second_order()
3707        #domain.distribute_to_vertices_and_edges()
3708        #print domain.quantities['stage'].vertex_values[:,0]
3709
3710        #Evolution
3711        for t in domain.evolve(yieldstep = 0.05, finaltime = 0.5):
3712            pass
3713
3714
3715        assert allclose(domain.quantities['stage'].centroid_values,
3716     [-0.02907028, -0.01475478, -0.02973417, -0.01447186, -0.02932665, -0.01428285,
3717      -0.02901975, -0.0141361,  -0.02898816, -0.01418135, -0.02961409, -0.01403487,
3718      -0.07597998, -0.06252591, -0.07664854, -0.06312532, -0.07638287, -0.06265139,
3719      -0.07571145, -0.06235231, -0.0756817,  -0.06245309, -0.07652292, -0.06289946,
3720      -0.12367464, -0.11088981, -0.12237277, -0.11115338, -0.1218934,  -0.1107174,
3721      -0.12081485, -0.11000491, -0.12038451, -0.11010335, -0.12102113, -0.11012105,
3722      -0.16909116, -0.15831543, -0.16730214, -0.15786249, -0.1665493,  -0.15697919,
3723      -0.16496618, -0.15559852, -0.16338679, -0.15509088, -0.16364092, -0.15424423,
3724      -0.18771107, -0.19903904, -0.18903759, -0.19858437, -0.18701552, -0.19697797,
3725      -0.1833593,  -0.19505871, -0.1818806,  -0.19418042, -0.18586159, -0.19576946,
3726      -0.13986873, -0.14170053, -0.14132188, -0.14560674, -0.14095617, -0.14373292,
3727      -0.13785933, -0.14033364, -0.13592955, -0.13936356, -0.13596008, -0.14216296])
3728
3729        assert allclose(domain.quantities['xmomentum'].centroid_values,
3730     [ 0.00831121,  0.00317948,  0.00731797,  0.00334939,  0.00764717,  0.00348053,
3731       0.00788729,  0.00356522,  0.00780649,  0.00341919,  0.00693525,  0.00310375,
3732       0.02166196,  0.01421475,  0.02017737,  0.01316839,  0.02037015,  0.01368659,
3733       0.02106,     0.01399161,  0.02074514,  0.01354935,  0.01887407,  0.0123113,
3734       0.03775083,  0.02855197,  0.03689337,  0.02759782,  0.03732848,  0.02812072,
3735       0.03872545,  0.02913348,  0.03880939,  0.02803804,  0.03546499,  0.0260039,
3736       0.0632131,   0.04730634,  0.0576324,   0.04592336,  0.05790921,  0.04690514,
3737       0.05986467,  0.04871165,  0.06170068,  0.04811572,  0.05657041,  0.04416292,
3738       0.08489642,  0.07188097,  0.07835261,  0.06843406,  0.07986412,  0.0698247,
3739       0.08201071,  0.07216756,  0.08378418,  0.07273624,  0.080399,    0.06645841,
3740       0.01631548,  0.04691608,  0.0206632,   0.044409,    0.02115518,  0.04560305,
3741       0.02160608,  0.04663725,  0.02174734,  0.04795559,  0.02281427,  0.05667111])
3742
3743
3744        assert allclose(domain.quantities['ymomentum'].centroid_values,
3745     [ 1.45876601e-004, -3.24627393e-004, -1.57572719e-004, -2.92790187e-004,
3746      -9.90988382e-005, -3.06677335e-004, -1.62493106e-004, -3.71310004e-004,
3747      -1.99445058e-004, -3.28493467e-004,  6.68217349e-005, -8.42042805e-006,
3748       5.05093371e-004, -1.42842214e-004, -6.81454718e-005, -5.02084057e-004,
3749      -8.50583861e-005, -4.65443981e-004, -1.96406564e-004, -5.88889562e-004,
3750      -2.70160173e-004, -5.35485454e-004,  2.60780997e-004,  3.12145471e-005,
3751       5.16189608e-004,  1.07069062e-004,  9.29989252e-005, -3.71211119e-004,
3752       1.16350246e-004, -3.82407830e-004, -1.62077969e-004, -6.30906636e-004,
3753      -4.74025708e-004, -6.94463009e-004,  6.15092843e-005,  2.22106820e-004,
3754      -6.29589294e-004,  2.43611937e-004, -5.88125094e-004, -6.94293192e-005,
3755      -4.17914641e-004,  6.64609019e-005, -7.68334577e-004, -3.40232101e-004,
3756      -1.67424308e-003, -7.39485066e-004, -1.59966988e-003,  5.68262838e-005,
3757      -1.48470633e-003, -1.84554882e-003, -2.27200099e-003, -1.67506848e-003,
3758      -1.95610258e-003, -1.47638801e-003, -1.73779477e-003, -1.85498791e-003,
3759      -2.01357843e-003, -2.17675471e-003, -1.65783870e-003, -1.15818681e-003,
3760      -1.18663036e-003, -2.94229849e-003, -3.59309018e-003, -5.13496584e-003,
3761      -6.17359400e-003, -5.98761937e-003, -6.00540116e-003, -5.01121966e-003,
3762      -4.50964850e-003, -3.06319963e-003,  6.08950810e-004, -4.79537921e-004])
3763
3764        os.remove(domain.get_name() + '.sww')
3765
3766
3767
3768    def NOtest_bedslope_problem_second_order_more_steps_feb_2007(self):
3769        """test_bedslope_problem_second_order_more_steps_feb_2007
3770
3771        Test shallow water finite volumes, using parameters from
3772        feb 2007 rather than backward compatibility ad infinitum
3773       
3774        """
3775        from mesh_factory import rectangular
3776        from Numeric import array
3777
3778        #Create basic mesh
3779        points, vertices, boundary = rectangular(6, 6)
3780
3781        #Create shallow water domain
3782        domain = Domain(points, vertices, boundary)
3783        domain.smooth = False
3784        domain.default_order = 2
3785        domain.beta_w      = 0.9
3786        domain.beta_w_dry  = 0.9
3787        domain.beta_uh     = 0.9
3788        domain.beta_uh_dry = 0.9
3789        domain.beta_vh     = 0.9
3790        domain.beta_vh_dry = 0.9
3791        domain.beta_h = 0.0 #Use first order in h-limiter
3792        domain.H0 = 0.001
3793        domain.tight_slope_limiters = 1
3794
3795        #Bed-slope and friction at vertices (and interpolated elsewhere)
3796        def x_slope(x, y):
3797            return -x/3
3798
3799        domain.set_quantity('elevation', x_slope)
3800
3801        # Boundary conditions
3802        Br = Reflective_boundary(domain)
3803        domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
3804
3805        #Initial condition
3806        domain.set_quantity('stage', expression = 'elevation + 0.05')
3807        domain.check_integrity()
3808
3809        assert allclose(domain.quantities['stage'].centroid_values,
3810                        [0.01296296, 0.03148148, 0.01296296,
3811                        0.03148148, 0.01296296, 0.03148148,
3812                        0.01296296, 0.03148148, 0.01296296,
3813                        0.03148148, 0.01296296, 0.03148148,
3814                        -0.04259259, -0.02407407, -0.04259259,
3815                        -0.02407407, -0.04259259, -0.02407407,
3816                        -0.04259259, -0.02407407, -0.04259259,
3817                        -0.02407407, -0.04259259, -0.02407407,
3818                        -0.09814815, -0.07962963, -0.09814815,
3819                        -0.07962963, -0.09814815, -0.07962963,
3820                        -0.09814815, -0.07962963, -0.09814815,
3821                        -0.07962963, -0.09814815, -0.07962963,
3822                        -0.1537037 , -0.13518519, -0.1537037,
3823                        -0.13518519, -0.1537037, -0.13518519,
3824                        -0.1537037 , -0.13518519, -0.1537037,
3825                        -0.13518519, -0.1537037, -0.13518519,
3826                        -0.20925926, -0.19074074, -0.20925926,
3827                        -0.19074074, -0.20925926, -0.19074074,
3828                        -0.20925926, -0.19074074, -0.20925926,
3829                        -0.19074074, -0.20925926, -0.19074074,
3830                        -0.26481481, -0.2462963, -0.26481481,
3831                        -0.2462963, -0.26481481, -0.2462963,
3832                        -0.26481481, -0.2462963, -0.26481481,
3833                        -0.2462963, -0.26481481, -0.2462963])
3834
3835
3836        #print domain.quantities['stage'].extrapolate_second_order()
3837        #domain.distribute_to_vertices_and_edges()
3838        #print domain.quantities['stage'].vertex_values[:,0]
3839
3840        #Evolution
3841        for t in domain.evolve(yieldstep = 0.05, finaltime = 0.5):
3842            pass
3843
3844
3845        #print domain.quantities['stage'].centroid_values
3846           
3847        assert allclose(domain.quantities['stage'].centroid_values,
3848         [-0.03348416, -0.01749303, -0.03299091, -0.01739241, -0.03246447, -0.01732016,
3849          -0.03205390, -0.01717833, -0.03146383, -0.01699831, -0.03076577, -0.01671795,
3850          -0.07952656, -0.06684763, -0.07721455, -0.06668388, -0.07632976, -0.06600113,
3851          -0.07523678, -0.06546373, -0.07447040, -0.06508861, -0.07438723, -0.06359288,
3852          -0.12526729, -0.11205668, -0.12179433, -0.11068104, -0.12048395, -0.10968948,
3853          -0.11912023, -0.10862628, -0.11784090, -0.10803744, -0.11790629, -0.10742354,
3854          -0.16859613, -0.15427413, -0.16664444, -0.15464452, -0.16570816, -0.15327556,
3855          -0.16409162, -0.15204092, -0.16264608, -0.15102139, -0.16162736, -0.14969205,
3856          -0.18736511, -0.19874036, -0.18811230, -0.19758289, -0.18590182, -0.19580301,
3857          -0.18234588, -0.19423215, -0.18100376, -0.19380116, -0.18509710, -0.19501636,
3858          -0.13982382, -0.14166819, -0.14132775, -0.14528694, -0.14096905, -0.14351126,
3859          -0.13800356, -0.14027920, -0.13613538, -0.13936795, -0.13621902, -0.14204982])
3860
3861                     
3862        assert allclose(domain.quantities['xmomentum'].centroid_values,
3863      [0.00600290,  0.00175780,  0.00591905,  0.00190903,  0.00644462,  0.00203095,
3864       0.00684561,  0.00225089,  0.00708208,  0.00236235,  0.00649095,  0.00222343,
3865       0.02068693,  0.01164034,  0.01983343,  0.01159526,  0.02044611,  0.01233252,
3866       0.02135685,  0.01301289,  0.02161290,  0.01260280,  0.01867612,  0.01133078,
3867       0.04091313,  0.02668283,  0.03634781,  0.02733469,  0.03767692,  0.02836840,
3868       0.03906338,  0.02958073,  0.04025669,  0.02953292,  0.03665616,  0.02583565,
3869       0.06314558,  0.04830935,  0.05663609,  0.04564362,  0.05756200,  0.04739673,
3870       0.05967379,  0.04919083,  0.06124330,  0.04965808,  0.05879240,  0.04629319,
3871       0.08220739,  0.06924725,  0.07713556,  0.06782640,  0.07909499,  0.06992544,
3872       0.08116621,  0.07210181,  0.08281548,  0.07222669,  0.07941059,  0.06755612,
3873       0.01581588,  0.04533609,  0.02017939,  0.04342565,  0.02073232,  0.04476108,
3874       0.02117439,  0.04573358,  0.02129473,  0.04694267,  0.02220398,  0.05533458])
3875
3876       
3877        assert allclose(domain.quantities['ymomentum'].centroid_values,
3878     [-7.65882069e-005, -1.46087080e-004, -1.09630102e-004, -7.80950424e-005,
3879      -1.15922807e-005, -9.09134899e-005, -1.35994542e-004, -1.95673476e-004,
3880      -4.25779199e-004, -2.95890312e-004, -4.00060341e-004, -9.42021290e-005,
3881      -3.41372596e-004, -1.54560195e-004, -2.94810038e-004, -1.08844546e-004,
3882      -6.97240892e-005,  3.50299623e-005, -2.40159184e-004, -2.01805883e-004,
3883      -7.60732405e-004, -5.10897642e-004, -1.00940001e-003, -1.38037759e-004,
3884      -1.06169131e-003, -3.12307760e-004, -9.90602307e-004, -4.21634250e-005,
3885      -6.02424239e-004,  1.52230578e-004, -7.63833035e-004, -1.10273481e-004,
3886      -1.40187071e-003, -5.57831837e-004, -1.63988285e-003, -2.48018092e-004,
3887      -1.83309840e-003, -6.19360836e-004, -1.29955242e-003, -3.76237145e-004,
3888      -1.00613007e-003, -8.63641918e-005, -1.13604124e-003, -3.90589728e-004,
3889      -1.91457355e-003, -9.43783961e-004, -2.28090840e-003, -5.79107025e-004,
3890      -1.54091533e-003, -2.39785792e-003, -2.47947427e-003, -2.02694009e-003,
3891      -2.10441194e-003, -1.82082650e-003, -1.80229336e-003, -2.10418336e-003,
3892      -1.93104408e-003, -2.23200334e-003, -1.57239706e-003, -1.31486358e-003,
3893      -1.17564993e-003, -2.85846494e-003, -3.52956754e-003, -5.12658193e-003,
3894      -6.24238960e-003, -6.01820113e-003, -6.09602201e-003, -5.04787190e-003,
3895      -4.59373845e-003, -3.01393146e-003,  5.08550095e-004, -4.35896549e-004])
3896
3897        os.remove(domain.get_name() + '.sww')
3898
3899
3900    def test_temp_play(self):
3901
3902        from mesh_factory import rectangular
3903        from Numeric import array
3904
3905        #Create basic mesh
3906        points, vertices, boundary = rectangular(5, 5)
3907
3908        #Create shallow water domain
3909        domain = Domain(points, vertices, boundary)
3910        domain.smooth = False
3911        domain.default_order=2
3912        domain.beta_w      = 0.9
3913        domain.beta_w_dry  = 0.9
3914        domain.beta_uh     = 0.9
3915        domain.beta_uh_dry = 0.9
3916        domain.beta_vh     = 0.9
3917        domain.beta_vh_dry = 0.9
3918        domain.beta_h = 0.0 #Use first order in h-limiter
3919       
3920        # FIXME (Ole): Need tests where these two are commented out
3921        domain.H0 = 0        # Backwards compatibility (6/2/7)       
3922        domain.tight_slope_limiters = 0 # Backwards compatibility (14/4/7)                         
3923       
3924
3925        #Bed-slope and friction at vertices (and interpolated elsewhere)
3926        def x_slope(x, y):
3927            return -x/3
3928
3929        domain.set_quantity('elevation', x_slope)
3930
3931        # Boundary conditions
3932        Br = Reflective_boundary(domain)
3933        domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
3934
3935        #Initial condition
3936        domain.set_quantity('stage', Constant_height(x_slope, 0.05))
3937        domain.check_integrity()
3938
3939        #Evolution
3940        for t in domain.evolve(yieldstep = 0.05, finaltime = 0.1):
3941            pass
3942
3943        assert allclose(domain.quantities['stage'].centroid_values[:4],
3944                        [0.00206836, 0.01296714, 0.00363415, 0.01438924])
3945        #print domain.quantities['xmomentum'].centroid_values[:4]
3946        assert allclose(domain.quantities['xmomentum'].centroid_values[:4],
3947                        [0.01360154, 0.00671133, 0.01264578, 0.00648503])
3948        assert allclose(domain.quantities['ymomentum'].centroid_values[:4],
3949                        [-1.19201077e-003, -7.23647546e-004,
3950                        -6.39083123e-005, 6.29815168e-005])
3951
3952        os.remove(domain.get_name() + '.sww')
3953
3954    def test_complex_bed(self):
3955        #No friction is tested here
3956
3957        from mesh_factory import rectangular
3958        from Numeric import array
3959
3960        N = 12
3961        points, vertices, boundary = rectangular(N, N/2, len1=1.2,len2=0.6,
3962                                                 origin=(-0.07, 0))
3963
3964
3965        domain = Domain(points, vertices, boundary)
3966        domain.smooth = False
3967        domain.default_order=2
3968
3969
3970        inflow_stage = 0.1
3971        Z = Weir(inflow_stage)
3972        domain.set_quantity('elevation', Z)
3973
3974        Br = Reflective_boundary(domain)
3975        Bd = Dirichlet_boundary([inflow_stage, 0.0, 0.0])
3976        domain.set_boundary({'left': Bd, 'right': Br, 'bottom': Br, 'top': Br})
3977
3978        domain.set_quantity('stage', Constant_height(Z, 0.))
3979
3980        for t in domain.evolve(yieldstep = 0.02, finaltime = 0.2):
3981            pass
3982
3983
3984        #print domain.quantities['stage'].centroid_values
3985
3986        #FIXME: These numbers were from version before 25/10
3987        #assert allclose(domain.quantities['stage'].centroid_values,
3988# [3.95822638e-002,  5.61022588e-002,  4.66437868e-002,  5.73081011e-002,
3989#  4.72394613e-002,  5.74684939e-002,  4.74309483e-002,  5.77458084e-002,
3990#  4.80628177e-002,  5.85656225e-002,  4.90498542e-002,  6.02609831e-002,
3991#  1.18470315e-002,  1.75136443e-002,  1.18035266e-002,  2.15565695e-002,
3992#  1.31620268e-002,  2.14351640e-002,  1.32351076e-002,  2.15450687e-002,
3993#  1.36414028e-002,  2.24274619e-002,  1.51689511e-002,  2.21789655e-002,
3994# -7.54337535e-003, -6.86362021e-004, -7.74146760e-003, -1.83756530e-003,
3995# -8.16773628e-003, -4.49916813e-004, -8.08202599e-003, -3.91118720e-004,
3996# -8.10292716e-003, -3.88584984e-004, -7.35226124e-003,  2.73985295e-004,
3997#  1.86166683e-001,  8.74070369e-002,  1.86166712e-001,  8.74035875e-002,
3998#  6.11666935e-002, -3.76173225e-002, -6.38333276e-002, -3.76147365e-002,
3999#  6.11666725e-002,  8.73846774e-002,  1.86166697e-001,  8.74171550e-002,
4000# -4.83333333e-002,  1.18333333e-001, -4.83333333e-002,  1.18333333e-001,
4001# -4.83333333e-002, -6.66666667e-003, -1.73333333e-001, -1.31666667e-001,
4002# -1.73333333e-001, -6.66666667e-003, -4.83333333e-002,  1.18333333e-001,
4003# -2.48333333e-001, -2.31666667e-001, -2.48333333e-001, -2.31666667e-001,
4004# -2.48333333e-001, -2.31666667e-001, -2.48333333e-001, -2.31666667e-001,
4005# -2.48333333e-001, -2.31666667e-001, -2.48333333e-001, -2.31666667e-001,
4006# -4.65000000e-001, -3.65000000e-001, -4.65000000e-001, -3.65000000e-001,
4007# -4.65000000e-001, -3.65000000e-001, -4.65000000e-001, -3.65000000e-001,
4008# -4.65000000e-001, -3.65000000e-001, -4.65000000e-001, -3.65000000e-001,
4009# -5.98333333e-001, -5.81666667e-001, -5.98333333e-001, -5.81666667e-001,
4010# -5.98333333e-001, -5.81666667e-001, -5.98333333e-001, -5.81666667e-001,
4011# -5.98333333e-001, -5.81666667e-001, -5.98333333e-001, -5.81666667e-001,
4012# -6.48333333e-001, -6.31666667e-001, -6.48333333e-001, -6.31666667e-001,
4013# -6.48333333e-001, -6.31666667e-001, -6.48333333e-001, -6.31666667e-001,
4014# -6.48333333e-001, -6.31666667e-001, -6.48333333e-001, -6.31666667e-001,
4015# -5.31666667e-001, -5.98333333e-001, -5.31666667e-001, -5.98333333e-001,
4016# -5.31666667e-001, -5.98333333e-001, -5.31666667e-001, -5.98333333e-001,
4017# -5.31666667e-001, -5.98333333e-001, -5.31666667e-001, -5.98333333e-001,
4018# -4.98333333e-001, -4.81666667e-001, -4.98333333e-001, -4.81666667e-001,
4019# -4.98333333e-001, -4.81666667e-001, -4.98333333e-001, -4.81666667e-001,
4020# -4.98333333e-001, -4.81666667e-001, -4.98333333e-001, -4.81666667e-001,
4021# -5.48333333e-001, -5.31666667e-001, -5.48333333e-001, -5.31666667e-001,
4022# -5.48333333e-001, -5.31666667e-001, -5.48333333e-001, -5.31666667e-001,
4023# -5.48333333e-001, -5.31666667e-001, -5.48333333e-001, -5.31666667e-001])
4024
4025        os.remove(domain.get_name() + '.sww')
4026
4027    def test_spatio_temporal_boundary_1(self):
4028        """Test that boundary values can be read from file and interpolated
4029        in both time and space.
4030
4031        Verify that the same steady state solution is arrived at and that
4032        time interpolation works.
4033
4034        The full solution history is not exactly the same as
4035        file boundary must read and interpolate from *smoothed* version
4036        as stored in sww.
4037        """
4038        import time
4039
4040        #Create sww file of simple propagation from left to right
4041        #through rectangular domain
4042
4043        from mesh_factory import rectangular
4044
4045        #Create basic mesh
4046        points, vertices, boundary = rectangular(3, 3)
4047
4048        #Create shallow water domain
4049        domain1 = Domain(points, vertices, boundary)
4050
4051        domain1.reduction = mean
4052        domain1.smooth = False  #Exact result
4053
4054        domain1.default_order = 2
4055        domain1.store = True
4056        domain1.set_datadir('.')
4057        domain1.set_name('spatio_temporal_boundary_source' + str(time.time()))
4058
4059        #FIXME: This is extremely important!
4060        #How can we test if they weren't stored?
4061        domain1.quantities_to_be_stored = ['stage', 'xmomentum', 'ymomentum']
4062
4063
4064        #Bed-slope and friction at vertices (and interpolated elsewhere)
4065        domain1.set_quantity('elevation', 0)
4066        domain1.set_quantity('friction', 0)
4067
4068        # Boundary conditions
4069        Br = Reflective_boundary(domain1)
4070        Bd = Dirichlet_boundary([0.3,0,0])
4071        domain1.set_boundary({'left': Bd, 'top': Bd, 'right': Br, 'bottom': Br})
4072        #Initial condition
4073        domain1.set_quantity('stage', 0)
4074        domain1.check_integrity()
4075
4076        finaltime = 5
4077        #Evolution  (full domain - large steps)
4078        for t in domain1.evolve(yieldstep = 0.671, finaltime = finaltime):
4079            pass
4080            #domain1.write_time()
4081
4082        cv1 = domain1.quantities['stage'].centroid_values
4083
4084
4085        #Create a triangle shaped domain (reusing coordinates from domain 1),
4086        #formed from the lower and right hand  boundaries and
4087        #the sw-ne diagonal
4088        #from domain 1. Call it domain2
4089
4090        points = [ [0,0], [1.0/3,0], [1.0/3,1.0/3],
4091                   [2.0/3,0], [2.0/3,1.0/3], [2.0/3,2.0/3],
4092                   [1,0], [1,1.0/3], [1,2.0/3], [1,1]]
4093
4094        vertices = [ [1,2,0], [3,4,1], [2,1,4], [4,5,2],
4095                     [6,7,3], [4,3,7], [7,8,4], [5,4,8], [8,9,5]]
4096
4097        boundary = { (0,1):'bottom', (1,1):'bottom', (4,1): 'bottom',
4098                     (4,2):'right', (6,2):'right', (8,2):'right',
4099                     (0,0):'diagonal', (3,0):'diagonal', (8,0):'diagonal'}
4100
4101        domain2 = Domain(points, vertices, boundary)
4102
4103        domain2.reduction = domain1.reduction
4104        domain2.smooth = False
4105        domain2.default_order = 2
4106
4107        #Bed-slope and friction at vertices (and interpolated elsewhere)
4108        domain2.set_quantity('elevation', 0)
4109        domain2.set_quantity('friction', 0)
4110        domain2.set_quantity('stage', 0)
4111
4112        # Boundary conditions
4113        Br = Reflective_boundary(domain2)
4114        #Bf = Spatio_temporal_boundary(domain1.get_name() + '.' +\
4115        #                              domain1.format, domain2)
4116        Bf = Field_boundary(domain1.get_name() + '.' +\
4117                            domain1.format, domain2)       
4118        domain2.set_boundary({'right':Br, 'bottom':Br, 'diagonal':Bf})
4119        domain2.check_integrity()
4120
4121
4122
4123        #Evolution (small steps)
4124        for t in domain2.evolve(yieldstep = 0.0711, finaltime = finaltime):
4125            pass
4126
4127
4128        #Use output from domain1 as spatio-temporal boundary for domain2
4129        #and verify that results at right hand side are close.
4130
4131        cv2 = domain2.quantities['stage'].centroid_values
4132
4133        #print take(cv1, (12,14,16))  #Right
4134        #print take(cv2, (4,6,8))
4135        #print take(cv1, (0,6,12))  #Bottom
4136        #print take(cv2, (0,1,4))
4137        #print take(cv1, (0,8,16))  #Diag
4138        #print take(cv2, (0,3,8))
4139
4140        assert allclose( take(cv1, (0,8,16)), take(cv2, (0,3,8))) #Diag
4141        assert allclose( take(cv1, (0,6,12)), take(cv2, (0,1,4))) #Bottom
4142        assert allclose( take(cv1, (12,14,16)), take(cv2, (4,6,8))) #RHS
4143
4144        #Cleanup
4145        os.remove(domain1.get_name() + '.' + domain1.format)
4146        os.remove(domain2.get_name() + '.' + domain2.format)       
4147
4148
4149
4150    def test_spatio_temporal_boundary_2(self):
4151        """Test that boundary values can be read from file and interpolated
4152        in both time and space.
4153        This is a more basic test, verifying that boundary object
4154        produces the expected results
4155
4156
4157        """
4158        import time
4159
4160        #Create sww file of simple propagation from left to right
4161        #through rectangular domain
4162
4163        from mesh_factory import rectangular
4164
4165        #Create basic mesh
4166        points, vertices, boundary = rectangular(3, 3)
4167
4168        #Create shallow water domain
4169        domain1 = Domain(points, vertices, boundary)
4170
4171        domain1.reduction = mean
4172        domain1.smooth = True #To mimic MOST output
4173
4174        domain1.default_order = 2
4175        domain1.store = True
4176        domain1.set_datadir('.')
4177        domain1.set_name('spatio_temporal_boundary_source' + str(time.time()))
4178
4179        #FIXME: This is extremely important!
4180        #How can we test if they weren't stored?
4181        domain1.quantities_to_be_stored = ['stage', 'xmomentum', 'ymomentum']
4182
4183
4184        #Bed-slope and friction at vertices (and interpolated elsewhere)
4185        domain1.set_quantity('elevation', 0)
4186        domain1.set_quantity('friction', 0)
4187
4188        # Boundary conditions
4189        Br = Reflective_boundary(domain1)
4190        Bd = Dirichlet_boundary([0.3,0,0])
4191        domain1.set_boundary({'left': Bd, 'top': Bd, 'right': Br, 'bottom': Br})
4192        #Initial condition
4193        domain1.set_quantity('stage', 0)
4194        domain1.check_integrity()
4195
4196        finaltime = 5
4197        #Evolution  (full domain - large steps)
4198        for t in domain1.evolve(yieldstep = 1, finaltime = finaltime):
4199            pass
4200            #domain1.write_time()
4201
4202
4203        #Create an triangle shaped domain (coinciding with some
4204        #coordinates from domain 1),
4205        #formed from the lower and right hand  boundaries and
4206        #the sw-ne diagonal
4207        #from domain 1. Call it domain2
4208
4209        points = [ [0,0], [1.0/3,0], [1.0/3,1.0/3],
4210                   [2.0/3,0], [2.0/3,1.0/3], [2.0/3,2.0/3],
4211                   [1,0], [1,1.0/3], [1,2.0/3], [1,1]]
4212
4213        vertices = [ [1,2,0],
4214                     [3,4,1], [2,1,4], [4,5,2],
4215                     [6,7,3], [4,3,7], [7,8,4], [5,4,8], [8,9,5]]
4216
4217        boundary = { (0,1):'bottom', (1,1):'bottom', (4,1): 'bottom',
4218                     (4,2):'right', (6,2):'right', (8,2):'right',
4219                     (0,0):'diagonal', (3,0):'diagonal', (8,0):'diagonal'}
4220
4221        domain2 = Domain(points, vertices, boundary)
4222
4223        domain2.reduction = domain1.reduction
4224        domain2.smooth = False
4225        domain2.default_order = 2
4226
4227        #Bed-slope and friction at vertices (and interpolated elsewhere)
4228        domain2.set_quantity('elevation', 0)
4229        domain2.set_quantity('friction', 0)
4230        domain2.set_quantity('stage', 0)
4231
4232
4233        #Read results for specific timesteps t=1 and t=2
4234        from Scientific.IO.NetCDF import NetCDFFile
4235        fid = NetCDFFile(domain1.get_name() + '.' + domain1.format)
4236
4237        x = fid.variables['x'][:]
4238        y = fid.variables['y'][:]
4239        s1 = fid.variables['stage'][1,:]
4240        s2 = fid.variables['stage'][2,:]
4241        fid.close()
4242
4243        from Numeric import take, reshape, concatenate
4244        shp = (len(x), 1)
4245        points = concatenate( (reshape(x, shp), reshape(y, shp)), axis=1)
4246        #The diagonal points of domain 1 are 0, 5, 10, 15
4247
4248        #print points[0], points[5], points[10], points[15]
4249        assert allclose( take(points, [0,5,10,15]),
4250                         [[0,0], [1.0/3, 1.0/3], [2.0/3, 2.0/3], [1,1]])
4251
4252
4253        # Boundary conditions
4254        Br = Reflective_boundary(domain2)
4255        #Bf = Spatio_temporal_boundary(domain1.get_name() + '.' + domain1.format,
4256        #                              domain2)
4257        Bf = Field_boundary(domain1.get_name() + '.' + domain1.format,
4258                            domain2)       
4259        domain2.set_boundary({'right':Br, 'bottom':Br, 'diagonal':Bf})
4260        domain2.check_integrity()
4261
4262        #Test that interpolation points are the mid points of the all boundary
4263        #segments
4264
4265        boundary_midpoints = [[1.0/6, 0], [1.0/2, 0], [5.0/6,0],
4266                              [1.0, 1.0/6], [1.0, 1.0/2], [1.0, 5.0/6],
4267                              [1.0/6, 1.0/6], [0.5, 0.5], [5.0/6, 5.0/6]]
4268
4269        boundary_midpoints.sort()
4270        R = Bf.F.interpolation_points.tolist()
4271        R.sort()
4272        assert allclose(boundary_midpoints, R)
4273
4274        #Check spatially interpolated output at time == 1
4275        domain2.time = 1
4276
4277        #First diagonal midpoint
4278        R0 = Bf.evaluate(0,0)
4279        assert allclose(R0[0], (s1[0] + s1[5])/2)
4280
4281        #Second diagonal midpoint
4282        R0 = Bf.evaluate(3,0)
4283        assert allclose(R0[0], (s1[5] + s1[10])/2)
4284
4285        #First diagonal midpoint
4286        R0 = Bf.evaluate(8,0)
4287        assert allclose(R0[0], (s1[10] + s1[15])/2)
4288
4289        #Check spatially interpolated output at time == 2
4290        domain2.time = 2
4291
4292        #First diagonal midpoint
4293        R0 = Bf.evaluate(0,0)
4294        assert allclose(R0[0], (s2[0] + s2[5])/2)
4295
4296        #Second diagonal midpoint
4297        R0 = Bf.evaluate(3,0)
4298        assert allclose(R0[0], (s2[5] + s2[10])/2)
4299
4300        #First diagonal midpoint
4301        R0 = Bf.evaluate(8,0)
4302        assert allclose(R0[0], (s2[10] + s2[15])/2)
4303
4304
4305        #Now check temporal interpolation
4306
4307        domain2.time = 1 + 2.0/3
4308
4309        #First diagonal midpoint
4310        R0 = Bf.evaluate(0,0)
4311        assert allclose(R0[0], ((s1[0] + s1[5])/2 + 2.0*(s2[0] + s2[5])/2)/3)
4312
4313        #Second diagonal midpoint
4314        R0 = Bf.evaluate(3,0)
4315        assert allclose(R0[0], ((s1[5] + s1[10])/2 + 2.0*(s2[5] + s2[10])/2)/3)
4316
4317        #First diagonal midpoint
4318        R0 = Bf.evaluate(8,0)
4319        assert allclose(R0[0], ((s1[10] + s1[15])/2 + 2.0*(s2[10] + s2[15])/2)/3)
4320
4321
4322
4323        #Cleanup
4324        os.remove(domain1.get_name() + '.' + domain1.format)
4325
4326
4327    def test_spatio_temporal_boundary_3(self):
4328        """Test that boundary values can be read from file and interpolated
4329        in both time and space.
4330        This is a more basic test, verifying that boundary object
4331        produces the expected results
4332
4333        This tests adjusting using mean_stage
4334
4335        """
4336
4337        import time
4338
4339        mean_stage = 5.2 # Adjust stage by this amount in boundary
4340
4341        #Create sww file of simple propagation from left to right
4342        #through rectangular domain
4343
4344        from mesh_factory import rectangular
4345
4346        #Create basic mesh
4347        points, vertices, boundary = rectangular(3, 3)
4348
4349        #Create shallow water domain
4350        domain1 = Domain(points, vertices, boundary)
4351
4352        domain1.reduction = mean
4353        domain1.smooth = True #To mimic MOST output
4354
4355        domain1.default_order = 2
4356        domain1.store = True
4357        domain1.set_datadir('.')
4358        domain1.set_name('spatio_temporal_boundary_source' + str(time.time()))
4359
4360        #FIXME: This is extremely important!
4361        #How can we test if they weren't stored?
4362        domain1.quantities_to_be_stored = ['stage', 'xmomentum', 'ymomentum']
4363
4364
4365        #Bed-slope and friction at vertices (and interpolated elsewhere)
4366        domain1.set_quantity('elevation', 0)
4367        domain1.set_quantity('friction', 0)
4368
4369        # Boundary conditions
4370        Br = Reflective_boundary(domain1)
4371        Bd = Dirichlet_boundary([0.3,0,0])
4372        domain1.set_boundary({'left': Bd, 'top': Bd, 'right': Br, 'bottom': Br})
4373        #Initial condition
4374        domain1.set_quantity('stage', 0)
4375        domain1.check_integrity()
4376
4377        finaltime = 5
4378        #Evolution  (full domain - large steps)
4379        for t in domain1.evolve(yieldstep = 1, finaltime = finaltime):
4380            pass
4381            #domain1.write_time()
4382
4383
4384        #Create an triangle shaped domain (coinciding with some
4385        #coordinates from domain 1),
4386        #formed from the lower and right hand  boundaries and
4387        #the sw-ne diagonal
4388        #from domain 1. Call it domain2
4389
4390        points = [ [0,0],
4391                   [1.0/3,0], [1.0/3,1.0/3],
4392                   [2.0/3,0], [2.0/3,1.0/3], [2.0/3,2.0/3],
4393                   [1,0],     [1,1.0/3],     [1,2.0/3],     [1,1]] 
4394                   
4395        vertices = [ [1,2,0],
4396                     [3,4,1], [2,1,4], [4,5,2],
4397                     [6,7,3], [4,3,7], [7,8,4], [5,4,8], [8,9,5]]
4398
4399        boundary = { (0,1):'bottom', (1,1):'bottom', (4,1): 'bottom',
4400                     (4,2):'right', (6,2):'right', (8,2):'right',
4401                     (0,0):'diagonal', (3,0):'diagonal', (8,0):'diagonal'}
4402
4403        domain2 = Domain(points, vertices, boundary)
4404
4405        domain2.reduction = domain1.reduction
4406        domain2.smooth = False
4407        domain2.default_order = 2
4408
4409        #Bed-slope and friction at vertices (and interpolated elsewhere)
4410        domain2.set_quantity('elevation', 0)
4411        domain2.set_quantity('friction', 0)
4412        domain2.set_quantity('stage', 0)
4413
4414
4415        #Read results for specific timesteps t=1 and t=2
4416        from Scientific.IO.NetCDF import NetCDFFile
4417        fid = NetCDFFile(domain1.get_name() + '.' + domain1.format)
4418
4419        x = fid.variables['x'][:]
4420        y = fid.variables['y'][:]
4421        s1 = fid.variables['stage'][1,:]
4422        s2 = fid.variables['stage'][2,:]
4423        fid.close()
4424
4425        from Numeric import take, reshape, concatenate
4426        shp = (len(x), 1)
4427        points = concatenate( (reshape(x, shp), reshape(y, shp)), axis=1)
4428        #The diagonal points of domain 1 are 0, 5, 10, 15
4429
4430        #print points[0], points[5], points[10], points[15]
4431        assert allclose( take(points, [0,5,10,15]),
4432                         [[0,0], [1.0/3, 1.0/3], [2.0/3, 2.0/3], [1,1]])
4433
4434
4435        # Boundary conditions
4436        Br = Reflective_boundary(domain2)
4437        #Bf = Spatio_temporal_boundary(domain1.get_name() + '.' + domain1.format,
4438        #                              domain2)
4439        Bf = Field_boundary(domain1.get_name() + '.' + domain1.format,
4440                            domain2, mean_stage=mean_stage)
4441       
4442        domain2.set_boundary({'right':Br, 'bottom':Br, 'diagonal':Bf})
4443        domain2.check_integrity()
4444
4445        #Test that interpolation points are the mid points of the all boundary
4446        #segments
4447
4448        boundary_midpoints = [[1.0/6, 0], [1.0/2, 0], [5.0/6,0],
4449                              [1.0, 1.0/6], [1.0, 1.0/2], [1.0, 5.0/6],
4450                              [1.0/6, 1.0/6], [0.5, 0.5], [5.0/6, 5.0/6]]
4451
4452        boundary_midpoints.sort()
4453        R = Bf.F.interpolation_points.tolist()
4454        R.sort()
4455        assert allclose(boundary_midpoints, R)
4456
4457        #Check spatially interpolated output at time == 1
4458        domain2.time = 1
4459
4460        #First diagonal midpoint
4461        R0 = Bf.evaluate(0,0)
4462        assert allclose(R0[0], (s1[0] + s1[5])/2 + mean_stage)
4463
4464        #Second diagonal midpoint
4465        R0 = Bf.evaluate(3,0)
4466        assert allclose(R0[0], (s1[5] + s1[10])/2 + mean_stage)
4467
4468        #First diagonal midpoint
4469        R0 = Bf.evaluate(8,0)
4470        assert allclose(R0[0], (s1[10] + s1[15])/2 + mean_stage)
4471
4472        #Check spatially interpolated output at time == 2
4473        domain2.time = 2
4474
4475        #First diagonal midpoint
4476        R0 = Bf.evaluate(0,0)
4477        assert allclose(R0[0], (s2[0] + s2[5])/2 + mean_stage)
4478
4479        #Second diagonal midpoint
4480        R0 = Bf.evaluate(3,0)
4481        assert allclose(R0[0], (s2[5] + s2[10])/2 + mean_stage)
4482
4483        #First diagonal midpoint
4484        R0 = Bf.evaluate(8,0)
4485        assert allclose(R0[0], (s2[10] + s2[15])/2 + mean_stage)
4486
4487
4488        #Now check temporal interpolation
4489
4490        domain2.time = 1 + 2.0/3
4491
4492        #First diagonal midpoint
4493        R0 = Bf.evaluate(0,0)
4494        assert allclose(R0[0], ((s1[0] + s1[5])/2 + 2.0*(s2[0] + s2[5])/2)/3 + mean_stage)
4495
4496        #Second diagonal midpoint
4497        R0 = Bf.evaluate(3,0)
4498        assert allclose(R0[0], ((s1[5] + s1[10])/2 + 2.0*(s2[5] + s2[10])/2)/3 + mean_stage)
4499
4500        #First diagonal midpoint
4501        R0 = Bf.evaluate(8,0)
4502        assert allclose(R0[0], ((s1[10] + s1[15])/2 + 2.0*(s2[10] + s2[15])/2)/3 + mean_stage)
4503
4504
4505        #Cleanup
4506        os.remove(domain1.get_name() + '.' + domain1.format)
4507
4508
4509    def test_spatio_temporal_boundary_outside(self):
4510        """Test that field_boundary catches if a point is outside the sww that defines it
4511        """
4512
4513        import time
4514        #Create sww file of simple propagation from left to right
4515        #through rectangular domain
4516
4517        from mesh_factory import rectangular
4518
4519        #Create basic mesh
4520        points, vertices, boundary = rectangular(3, 3)
4521
4522        #Create shallow water domain
4523        domain1 = Domain(points, vertices, boundary)
4524
4525        domain1.reduction = mean
4526        domain1.smooth = True #To mimic MOST output
4527
4528        domain1.default_order = 2
4529        domain1.store = True
4530        domain1.set_datadir('.')
4531        domain1.set_name('spatio_temporal_boundary_source' + str(time.time()))
4532
4533        #FIXME: This is extremely important!
4534        #How can we test if they weren't stored?
4535        domain1.quantities_to_be_stored = ['stage', 'xmomentum', 'ymomentum']
4536
4537
4538        #Bed-slope and friction at vertices (and interpolated elsewhere)
4539        domain1.set_quantity('elevation', 0)
4540        domain1.set_quantity('friction', 0)
4541
4542        # Boundary conditions
4543        Br = Reflective_boundary(domain1)
4544        Bd = Dirichlet_boundary([0.3,0,0])
4545        domain1.set_boundary({'left': Bd, 'top': Bd, 'right': Br, 'bottom': Br})
4546        #Initial condition
4547        domain1.set_quantity('stage', 0)
4548        domain1.check_integrity()
4549
4550        finaltime = 5
4551        #Evolution  (full domain - large steps)
4552        for t in domain1.evolve(yieldstep = 1, finaltime = finaltime):
4553            pass
4554            #domain1.write_time()
4555
4556
4557        #Create an triangle shaped domain (coinciding with some
4558        #coordinates from domain 1, but one edge outside!),
4559        #formed from the lower and right hand  boundaries and
4560        #the sw-ne diagonal as in the previous test but scaled
4561        #in the x direction by a factor of 2
4562
4563        points = [ [0,0],
4564                   [2.0/3,0], [2.0/3,1.0/3],
4565                   [4.0/3,0], [4.0/3,1.0/3], [4.0/3,2.0/3],
4566                   [2,0],     [2,1.0/3],     [2,2.0/3],     [2,1] 
4567                   ]
4568
4569        vertices = [ [1,2,0],
4570                     [3,4,1], [2,1,4], [4,5,2],
4571                     [6,7,3], [4,3,7], [7,8,4], [5,4,8], [8,9,5]]
4572
4573        boundary = { (0,1):'bottom', (1,1):'bottom', (4,1): 'bottom',
4574                     (4,2):'right', (6,2):'right', (8,2):'right',
4575                     (0,0):'diagonal', (3,0):'diagonal', (8,0):'diagonal'}
4576
4577        domain2 = Domain(points, vertices, boundary)
4578
4579        domain2.reduction = domain1.reduction
4580        domain2.smooth = False
4581        domain2.default_order = 2
4582
4583        #Bed-slope and friction at vertices (and interpolated elsewhere)
4584        domain2.set_quantity('elevation', 0)
4585        domain2.set_quantity('friction', 0)
4586        domain2.set_quantity('stage', 0)
4587
4588
4589        #Read results for specific timesteps t=1 and t=2
4590        from Scientific.IO.NetCDF import NetCDFFile
4591        fid = NetCDFFile(domain1.get_name() + '.' + domain1.format)
4592
4593        x = fid.variables['x'][:]
4594        y = fid.variables['y'][:]
4595        s1 = fid.variables['stage'][1,:]
4596        s2 = fid.variables['stage'][2,:]
4597        fid.close()
4598
4599        from Numeric import take, reshape, concatenate
4600        shp = (len(x), 1)
4601        points = concatenate( (reshape(x, shp), reshape(y, shp)), axis=1)
4602        #The diagonal points of domain 1 are 0, 5, 10, 15
4603
4604        #print points[0], points[5], points[10], points[15]
4605        assert allclose( take(points, [0,5,10,15]),
4606                         [[0,0], [1.0/3, 1.0/3], [2.0/3, 2.0/3], [1,1]])
4607
4608
4609        # Boundary conditions
4610        Br = Reflective_boundary(domain2)
4611        #Bf = Spatio_temporal_boundary(domain1.get_name() + '.' + domain1.format,
4612        #                              domain2)
4613        Bf = Field_boundary(domain1.get_name() + '.' + domain1.format,
4614                            domain2, mean_stage=1)
4615       
4616        domain2.set_boundary({'right':Br, 'bottom':Br, 'diagonal':Bf})
4617        domain2.check_integrity()
4618
4619        try:
4620            for t in domain2.evolve(yieldstep = 1, finaltime = finaltime):
4621                pass
4622        except:
4623            pass
4624        else:
4625            msg = 'This should have caught NAN at boundary'
4626            raise Exception, msg
4627
4628
4629        #Cleanup
4630        os.remove(domain1.get_name() + '.' + domain1.format)
4631
4632
4633
4634    def test_tight_slope_limiters(self):
4635        """Test that new slope limiters (Feb 2007) don't induce extremely
4636        small timesteps. This test actually reveals the problem as it
4637        was in March-April 2007
4638        """
4639
4640        import time, os
4641        from Numeric import array, zeros, allclose, Float, concatenate
4642        from Scientific.IO.NetCDF import NetCDFFile
4643        from data_manager import get_dataobject, extent_sww
4644        from mesh_factory import rectangular
4645
4646       
4647        #Create basic mesh
4648        points, vertices, boundary = rectangular(2, 2)
4649
4650        #Create shallow water domain
4651        domain = Domain(points, vertices, boundary)
4652        domain.default_order = 2
4653
4654        # This will pass
4655        #domain.tight_slope_limiters = 1
4656        #domain.H0 = 0.01
4657       
4658        # This will fail
4659        #domain.tight_slope_limiters = 1
4660        #domain.H0 = 0.001
4661
4662        # This will pass provided C extension implements limiting of
4663        # momentum in _compute_speeds
4664        domain.tight_slope_limiters = 1
4665        domain.H0 = 0.001       
4666        domain.protect_against_isolated_degenerate_timesteps = True       
4667
4668        #Set some field values
4669        domain.set_quantity('elevation', lambda x,y: -x)
4670        domain.set_quantity('friction', 0.03)
4671
4672
4673        ######################
4674        # Boundary conditions
4675        B = Transmissive_boundary(domain)
4676        domain.set_boundary( {'left': B, 'right': B, 'top': B, 'bottom': B})
4677
4678
4679        ######################
4680        #Initial condition - with jumps
4681
4682
4683        bed = domain.quantities['elevation'].vertex_values
4684        stage = zeros(bed.shape, Float)
4685
4686        h = 0.3
4687        for i in range(stage.shape[0]):
4688            if i % 2 == 0:
4689                stage[i,:] = bed[i,:] + h
4690            else:
4691                stage[i,:] = bed[i,:]
4692
4693        domain.set_quantity('stage', stage)
4694
4695
4696        domain.distribute_to_vertices_and_edges()               
4697
4698       
4699
4700        domain.set_name('tight_limiters')
4701        domain.format = 'sww'
4702        domain.smooth = True
4703        domain.reduction = mean
4704        domain.set_datadir('.')
4705        domain.smooth = False
4706        domain.store = True
4707        domain.beta_h = 0
4708       
4709
4710        #Evolution
4711        for t in domain.evolve(yieldstep = 0.1, finaltime = 0.3):
4712           
4713            #domain.write_time(track_speeds=True)
4714            stage = domain.quantities['stage'].vertex_values
4715
4716            #Get NetCDF
4717            fid = NetCDFFile(domain.writer.filename, 'r')
4718            stage_file = fid.variables['stage']
4719           
4720            fid.close()
4721
4722        os.remove(domain.writer.filename)
4723
4724
4725    def test_pmesh2Domain(self):
4726         import os
4727         import tempfile
4728
4729         fileName = tempfile.mktemp(".tsh")
4730         file = open(fileName,"w")
4731         file.write("4 3 # <vertex #> <x> <y> [attributes]\n \
47320 0.0 0.0 0.0 0.0 0.01 \n \
47331 1.0 0.0 10.0 10.0 0.02  \n \
47342 0.0 1.0 0.0 10.0 0.03  \n \
47353 0.5 0.25 8.0 12.0 0.04  \n \
4736# Vert att title  \n \
4737elevation  \n \
4738stage  \n \
4739friction  \n \
47402 # <triangle #> [<vertex #>] [<neigbouring triangle #>]  \n\
47410 0 3 2 -1  -1  1 dsg\n\
47421 0 1 3 -1  0 -1   ole nielsen\n\
47434 # <segment #> <vertex #>  <vertex #> [boundary tag] \n\
47440 1 0 2 \n\
47451 0 2 3 \n\
47462 2 3 \n\
47473 3 1 1 \n\
47483 0 # <x> <y> [attributes] ...Mesh Vertices... \n \
47490 216.0 -86.0 \n \
47501 160.0 -167.0 \n \
47512 114.0 -91.0 \n \
47523 # <vertex #>  <vertex #> [boundary tag] ...Mesh Segments... \n \
47530 0 1 0 \n \
47541 1 2 0 \n \
47552 2 0 0 \n \
47560 # <x> <y> ...Mesh Holes... \n \
47570 # <x> <y> <attribute>...Mesh Regions... \n \
47580 # <x> <y> <attribute>...Mesh Regions, area... \n\
4759#Geo reference \n \
476056 \n \
4761140 \n \
4762120 \n")
4763         file.close()
4764
4765         tags = {}
4766         b1 =  Dirichlet_boundary(conserved_quantities = array([0.0]))
4767         b2 =  Dirichlet_boundary(conserved_quantities = array([1.0]))
4768         b3 =  Dirichlet_boundary(conserved_quantities = array([2.0]))
4769         tags["1"] = b1
4770         tags["2"] = b2
4771         tags["3"] = b3
4772
4773         #from anuga.abstract_2d_finite_volumes.pmesh2domain import pmesh_to_domain_instance
4774         #domain = pmesh_to_domain_instance(fileName, Domain)
4775
4776         domain = Domain(mesh_filename=fileName)
4777                         #verbose=True, use_cache=True)
4778         
4779         #print "domain.tagged_elements", domain.tagged_elements
4780         ## check the quantities
4781         #print domain.quantities['elevation'].vertex_values
4782         answer = [[0., 8., 0.],
4783                   [0., 10., 8.]]
4784         assert allclose(domain.quantities['elevation'].vertex_values,
4785                        answer)
4786
4787         #print domain.quantities['stage'].vertex_values
4788         answer = [[0., 12., 10.],
4789                   [0., 10., 12.]]
4790         assert allclose(domain.quantities['stage'].vertex_values,
4791                        answer)
4792
4793         #print domain.quantities['friction'].vertex_values
4794         answer = [[0.01, 0.04, 0.03],
4795                   [0.01, 0.02, 0.04]]
4796         assert allclose(domain.quantities['friction'].vertex_values,
4797                        answer)
4798
4799         #print domain.quantities['friction'].vertex_values
4800         assert allclose(domain.tagged_elements['dsg'][0],0)
4801         assert allclose(domain.tagged_elements['ole nielsen'][0],1)
4802
4803         self.failUnless( domain.boundary[(1, 0)]  == '1',
4804                          "test_tags_to_boundaries  failed. Single boundary wasn't added.")
4805         self.failUnless( domain.boundary[(1, 2)]  == '2',
4806                          "test_tags_to_boundaries  failed. Single boundary wasn't added.")
4807         self.failUnless( domain.boundary[(0, 1)]  == '3',
4808                          "test_tags_to_boundaries  failed. Single boundary wasn't added.")
4809         self.failUnless( domain.boundary[(0, 0)]  == 'exterior',
4810                          "test_tags_to_boundaries  failed. Single boundary wasn't added.")
4811         #print "domain.boundary",domain.boundary
4812         self.failUnless( len(domain.boundary)  == 4,
4813                          "test_pmesh2Domain Too many boundaries")
4814         #FIXME change to use get_xllcorner
4815         #print "d.geo_reference.xllcorner",domain.geo_reference.xllcorner
4816         self.failUnless(domain.geo_reference.xllcorner  == 140.0,
4817                          "bad geo_referece")
4818
4819
4820         #************
4821
4822   
4823         domain = Domain(fileName)
4824         
4825         #print "domain.tagged_elements", domain.tagged_elements
4826         ## check the quantities
4827         #print domain.quantities['elevation'].vertex_values
4828         answer = [[0., 8., 0.],
4829                   [0., 10., 8.]]
4830         assert allclose(domain.quantities['elevation'].vertex_values,
4831                        answer)
4832
4833         #print domain.quantities['stage'].vertex_values
4834         answer = [[0., 12., 10.],
4835                   [0., 10., 12.]]
4836         assert allclose(domain.quantities['stage'].vertex_values,
4837                        answer)
4838
4839         #print domain.quantities['friction'].vertex_values
4840         answer = [[0.01, 0.04, 0.03],
4841                   [0.01, 0.02, 0.04]]
4842         assert allclose(domain.quantities['friction'].vertex_values,
4843                        answer)
4844
4845         #print domain.quantities['friction'].vertex_values
4846         assert allclose(domain.tagged_elements['dsg'][0],0)
4847         assert allclose(domain.tagged_elements['ole nielsen'][0],1)
4848
4849         self.failUnless( domain.boundary[(1, 0)]  == '1',
4850                          "test_tags_to_boundaries  failed. Single boundary wasn't added.")
4851         self.failUnless( domain.boundary[(1, 2)]  == '2',
4852                          "test_tags_to_boundaries  failed. Single boundary wasn't added.")
4853         self.failUnless( domain.boundary[(0, 1)]  == '3',
4854                          "test_tags_to_boundaries  failed. Single boundary wasn't added.")
4855         self.failUnless( domain.boundary[(0, 0)]  == 'exterior',
4856                          "test_tags_to_boundaries  failed. Single boundary wasn't added.")
4857         #print "domain.boundary",domain.boundary
4858         self.failUnless( len(domain.boundary)  == 4,
4859                          "test_pmesh2Domain Too many boundaries")
4860         #FIXME change to use get_xllcorner
4861         #print "d.geo_reference.xllcorner",domain.geo_reference.xllcorner
4862         self.failUnless(domain.geo_reference.xllcorner  == 140.0,
4863                          "bad geo_referece")
4864         #************
4865         os.remove(fileName)
4866
4867        #-------------------------------------------------------------
4868       
4869if __name__ == "__main__":
4870
4871    suite = unittest.makeSuite(Test_Shallow_Water,'test')   
4872    #suite = unittest.makeSuite(Test_Shallow_Water,'test_tight_slope_limiters')
4873    #suite = unittest.makeSuite(Test_Shallow_Water,'test_get_maximum_inundation_from_sww')
4874    #suite = unittest.makeSuite(Test_Shallow_Water,'test_temp')   
4875   
4876
4877   
4878    runner = unittest.TextTestRunner(verbosity=1)   
4879    runner.run(suite)
Note: See TracBrowser for help on using the repository browser.