source: inundation/ga/storm_surge/pyvolution/test_shallow_water.py @ 1697

Last change on this file since 1697 was 1671, checked in by ole, 20 years ago

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