source: inundation/pyvolution/test_shallow_water.py @ 2126

Last change on this file since 2126 was 1835, checked in by ole, 19 years ago

Implemented tms file format (like sww without x, y)
Fixed broken test in sww2ers

File size: 112.4 KB
Line 
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*#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*#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
199        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
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
996        filename = 'test_windstress_from_file'
997        start = time.mktime(time.strptime('2000', '%Y'))
998        fid = open(filename + '.txt', 'w')
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
1011
1012        #Convert ASCII file to NetCDF (Which is what we really like!)
1013        from data_manager import timefile2netcdf       
1014        timefile2netcdf(filename)
1015
1016
1017       
1018        #Setup wind stress
1019        F = file_function(filename + '.tms', 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       
1036        W = Wind_stress(F)
1037
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
1069        #os.remove(filename + '.txt')
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],
1241                            location='centroids')
1242        L = domain.quantities['stage'].vertex_values
1243
1244        #First order
1245        domain.order = 1
1246        domain.distribute_to_vertices_and_edges()
1247        assert allclose(L[1], val1)
1248
1249        #Second order
1250        domain.order = 2
1251        domain.distribute_to_vertices_and_edges()
1252        assert allclose(L[1], [2.2, 4.9, 4.9])
1253
1254
1255
1256    def test_distribute_away_from_bed(self):
1257        #Using test data generated by pyvolution-2
1258        #Assuming no friction and flat bed (0.0)
1259
1260        a = [0.0, 0.0]
1261        b = [0.0, 2.0]
1262        c = [2.0, 0.0]
1263        d = [0.0, 4.0]
1264        e = [2.0, 2.0]
1265        f = [4.0, 0.0]
1266
1267        points = [a, b, c, d, e, f]
1268        #bac, bce, ecf, dbe
1269        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
1270
1271        domain = Domain(points, vertices)
1272        L = domain.quantities['stage'].vertex_values
1273
1274        def stage(x,y):
1275            return x**2
1276
1277        domain.set_quantity('stage', stage, location='centroids')
1278
1279        a, b = domain.quantities['stage'].compute_gradients()
1280        assert allclose(a[1], 3.33333334)
1281        assert allclose(b[1], 0.0)
1282
1283        domain.order = 1
1284        domain.distribute_to_vertices_and_edges()
1285        assert allclose(L[1], 1.77777778)
1286
1287        domain.order = 2
1288        domain.distribute_to_vertices_and_edges()
1289        assert allclose(L[1], [0.57777777, 2.37777778, 2.37777778])
1290
1291
1292
1293    def test_distribute_away_from_bed1(self):
1294        #Using test data generated by pyvolution-2
1295        #Assuming no friction and flat bed (0.0)
1296
1297        a = [0.0, 0.0]
1298        b = [0.0, 2.0]
1299        c = [2.0, 0.0]
1300        d = [0.0, 4.0]
1301        e = [2.0, 2.0]
1302        f = [4.0, 0.0]
1303
1304        points = [a, b, c, d, e, f]
1305        #bac, bce, ecf, dbe
1306        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
1307
1308        domain = Domain(points, vertices)
1309        L = domain.quantities['stage'].vertex_values
1310
1311        def stage(x,y):
1312            return x**4+y**2
1313
1314        domain.set_quantity('stage', stage, location='centroids')
1315        #print domain.quantities['stage'].centroid_values
1316
1317        a, b = domain.quantities['stage'].compute_gradients()
1318        assert allclose(a[1], 25.18518519)
1319        assert allclose(b[1], 3.33333333)
1320
1321        domain.order = 1
1322        domain.distribute_to_vertices_and_edges()
1323        assert allclose(L[1], 4.9382716)
1324
1325        domain.order = 2
1326        domain.distribute_to_vertices_and_edges()
1327        assert allclose(L[1], [1.07160494, 6.46058131, 7.28262855])
1328
1329
1330
1331    def test_distribute_near_bed(self):
1332        #Using test data generated by pyvolution-2
1333        #Assuming no friction and flat bed (0.0)
1334
1335        a = [0.0, 0.0]
1336        b = [0.0, 2.0]
1337        c = [2.0, 0.0]
1338        d = [0.0, 4.0]
1339        e = [2.0, 2.0]
1340        f = [4.0, 0.0]
1341
1342        points = [a, b, c, d, e, f]
1343        #bac, bce, ecf, dbe
1344        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
1345
1346        domain = Domain(points, vertices)
1347
1348
1349        #Set up for a gradient of (3,0) at mid triangle
1350        def slope(x, y):
1351            return 10*x
1352
1353        h = 0.1
1354        def stage(x,y):
1355            return slope(x,y)+h
1356
1357        domain.set_quantity('elevation', slope)
1358        domain.set_quantity('stage', stage, location='centroids')
1359
1360        #print domain.quantities['elevation'].centroid_values
1361        #print domain.quantities['stage'].centroid_values
1362
1363        E = domain.quantities['elevation'].vertex_values
1364        L = domain.quantities['stage'].vertex_values
1365
1366        #print E
1367        domain.order = 1
1368        domain.distribute_to_vertices_and_edges()
1369        ##assert allclose(L[1], [0.19999999, 20.05, 20.05])
1370        assert allclose(L[1], [0.1, 20.1, 20.1])
1371
1372        domain.order = 2
1373        domain.distribute_to_vertices_and_edges()
1374        assert allclose(L[1], [0.1, 20.1, 20.1])
1375
1376    def test_distribute_near_bed1(self):
1377        #Using test data generated by pyvolution-2
1378        #Assuming no friction and flat bed (0.0)
1379
1380        a = [0.0, 0.0]
1381        b = [0.0, 2.0]
1382        c = [2.0, 0.0]
1383        d = [0.0, 4.0]
1384        e = [2.0, 2.0]
1385        f = [4.0, 0.0]
1386
1387        points = [a, b, c, d, e, f]
1388        #bac, bce, ecf, dbe
1389        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
1390
1391        domain = Domain(points, vertices)
1392
1393
1394        #Set up for a gradient of (3,0) at mid triangle
1395        def slope(x, y):
1396            return x**4+y**2
1397
1398        h = 0.1
1399        def stage(x,y):
1400            return slope(x,y)+h
1401
1402        domain.set_quantity('elevation', slope)
1403        domain.set_quantity('stage', stage)
1404
1405        #print domain.quantities['elevation'].centroid_values
1406        #print domain.quantities['stage'].centroid_values
1407
1408        E = domain.quantities['elevation'].vertex_values
1409        L = domain.quantities['stage'].vertex_values
1410
1411        #print E
1412        domain.order = 1
1413        domain.distribute_to_vertices_and_edges()
1414        ##assert allclose(L[1], [4.19999999, 16.07142857, 20.02857143])
1415        assert allclose(L[1], [4.1, 16.1, 20.1])
1416
1417        domain.order = 2
1418        domain.distribute_to_vertices_and_edges()
1419        assert allclose(L[1], [4.1, 16.1, 20.1])
1420
1421
1422
1423    def test_second_order_distribute_real_data(self):
1424        #Using test data generated by pyvolution-2
1425        #Assuming no friction and flat bed (0.0)
1426
1427        a = [0.0, 0.0]
1428        b = [0.0, 1.0/5]
1429        c = [0.0, 2.0/5]
1430        d = [1.0/5, 0.0]
1431        e = [1.0/5, 1.0/5]
1432        f = [1.0/5, 2.0/5]
1433        g = [2.0/5, 2.0/5]
1434
1435        points = [a, b, c, d, e, f, g]
1436        #bae, efb, cbf, feg
1437        vertices = [ [1,0,4], [4,5,1], [2,1,5], [5,4,6]]
1438
1439        domain = Domain(points, vertices)
1440
1441        def slope(x, y):
1442            return -x/3
1443
1444        domain.set_quantity('elevation', slope)
1445        domain.set_quantity('stage',
1446                            [0.01298164, 0.00365611,
1447                             0.01440365, -0.0381856437096],
1448                            location='centroids')
1449        domain.set_quantity('xmomentum',
1450                            [0.00670439, 0.01263789,
1451                             0.00647805, 0.0178180740668],
1452                            location='centroids')
1453        domain.set_quantity('ymomentum',
1454                            [-7.23510980e-004, -6.30413883e-005,
1455                             6.30413883e-005, 0.000200907255866],
1456                            location='centroids')
1457
1458        E = domain.quantities['elevation'].vertex_values
1459        L = domain.quantities['stage'].vertex_values
1460        X = domain.quantities['xmomentum'].vertex_values
1461        Y = domain.quantities['ymomentum'].vertex_values
1462
1463        #print E
1464        domain.order = 2
1465        domain.beta_h = 0.0 #Use first order in h-limiter
1466        domain.distribute_to_vertices_and_edges()
1467
1468        #print L[1,:]
1469        #print X[1,:]
1470        #print Y[1,:]
1471
1472        assert allclose(L[1,:], [-0.00825735775384,
1473                                 -0.00801881482869,
1474                                 0.0272445025825])
1475        assert allclose(X[1,:], [0.0143507718962,
1476                                 0.0142502147066,
1477                                 0.00931268339717])
1478        assert allclose(Y[1,:], [-0.000117062180693,
1479                                 7.94434448109e-005,
1480                                 -0.000151505429018])
1481
1482
1483
1484    def test_balance_deep_and_shallow(self):
1485        """Test that balanced limiters preserve conserved quantites.
1486        """
1487        import copy
1488
1489        a = [0.0, 0.0]
1490        b = [0.0, 2.0]
1491        c = [2.0, 0.0]
1492        d = [0.0, 4.0]
1493        e = [2.0, 2.0]
1494        f = [4.0, 0.0]
1495
1496        points = [a, b, c, d, e, f]
1497
1498        #bac, bce, ecf, dbe
1499        elements = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4] ]
1500
1501        mesh = Domain(points, elements)
1502        mesh.check_integrity()
1503
1504        #Create a deliberate overshoot
1505        mesh.set_quantity('stage', [[3,0,3], [2,2,6], [5,3,8], [8,3,5]])
1506        mesh.set_quantity('elevation', 0) #Flat bed
1507        stage = mesh.quantities['stage']
1508
1509        ref_centroid_values = copy.copy(stage.centroid_values[:]) #Copy
1510
1511        #Limit
1512        mesh.distribute_to_vertices_and_edges()
1513
1514        #Assert that quantities are conserved
1515        from Numeric import sum
1516        for k in range(mesh.number_of_elements):
1517            assert allclose (ref_centroid_values[k],
1518                             sum(stage.vertex_values[k,:])/3)
1519
1520
1521        #Now try with a non-flat bed - closely hugging initial stage in places
1522        #This will create alphas in the range [0, 0.478260, 1]
1523        mesh.set_quantity('stage', [[3,0,3], [2,2,6], [5,3,8], [8,3,5]])
1524        mesh.set_quantity('elevation', [[0,0,0],
1525                                        [1.8,1.9,5.9],
1526                                        [4.6,0,0],
1527                                        [0,2,4]])
1528        stage = mesh.quantities['stage']
1529
1530        ref_centroid_values = copy.copy(stage.centroid_values[:]) #Copy
1531        ref_vertex_values = copy.copy(stage.vertex_values[:]) #Copy
1532
1533        #Limit
1534        mesh.distribute_to_vertices_and_edges()
1535
1536
1537        #Assert that all vertex quantities have changed
1538        for k in range(mesh.number_of_elements):
1539            #print ref_vertex_values[k,:], stage.vertex_values[k,:]
1540            assert not allclose (ref_vertex_values[k,:], stage.vertex_values[k,:])
1541        #and assert that quantities are still conserved
1542        from Numeric import sum
1543        for k in range(mesh.number_of_elements):
1544            assert allclose (ref_centroid_values[k],
1545                             sum(stage.vertex_values[k,:])/3)
1546
1547
1548        #Also check that Python and C version produce the same
1549        assert allclose (stage.vertex_values,
1550                         [[2,2,2],
1551                          [1.93333333, 2.03333333, 6.03333333],
1552                          [6.93333333, 4.53333333, 4.53333333],
1553                          [5.33333333, 5.33333333, 5.33333333]])
1554
1555
1556
1557
1558    def test_conservation_1(self):
1559        """Test that stage is conserved globally
1560
1561        This one uses a flat bed, reflective bdries and a suitable
1562        initial condition
1563        """
1564        from mesh_factory import rectangular
1565        from shallow_water import Domain, Reflective_boundary,\
1566             Dirichlet_boundary, Constant_height
1567        from Numeric import array
1568
1569        #Create basic mesh
1570        points, vertices, boundary = rectangular(6, 6)
1571
1572        #Create shallow water domain
1573        domain = Domain(points, vertices, boundary)
1574        domain.smooth = False
1575        domain.default_order=2
1576
1577        #IC
1578        def x_slope(x, y):
1579            return x/3
1580
1581        domain.set_quantity('elevation', 0)
1582        domain.set_quantity('friction', 0)
1583        domain.set_quantity('stage', x_slope)
1584
1585        # Boundary conditions (reflective everywhere)
1586        Br = Reflective_boundary(domain)
1587        domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
1588
1589        domain.check_integrity()
1590
1591        #domain.visualise = True #If you want to take a sticky beak
1592
1593        initial_volume = domain.quantities['stage'].get_integral()
1594        initial_xmom = domain.quantities['xmomentum'].get_integral()
1595
1596        #print initial_xmom
1597
1598        #Evolution
1599        for t in domain.evolve(yieldstep = 0.05, finaltime = 5.0):
1600            volume =  domain.quantities['stage'].get_integral()
1601            assert allclose (volume, initial_volume)
1602
1603            #I don't believe that the total momentum should be the same
1604            #It starts with zero and ends with zero though
1605            #xmom = domain.quantities['xmomentum'].get_integral()
1606            #print xmom
1607            #assert allclose (xmom, initial_xmom)
1608
1609
1610
1611    def test_conservation_2(self):
1612        """Test that stage is conserved globally
1613
1614        This one uses a slopy bed, reflective bdries and a suitable
1615        initial condition
1616        """
1617        from mesh_factory import rectangular
1618        from shallow_water import Domain, Reflective_boundary,\
1619             Dirichlet_boundary, Constant_height
1620        from Numeric import array
1621
1622        #Create basic mesh
1623        points, vertices, boundary = rectangular(6, 6)
1624
1625        #Create shallow water domain
1626        domain = Domain(points, vertices, boundary)
1627        domain.smooth = False
1628        domain.default_order=2
1629
1630        #IC
1631        def x_slope(x, y):
1632            return x/3
1633
1634        domain.set_quantity('elevation', x_slope)
1635        domain.set_quantity('friction', 0)
1636        domain.set_quantity('stage', 0.4) #Steady
1637
1638        # Boundary conditions (reflective everywhere)
1639        Br = Reflective_boundary(domain)
1640        domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
1641
1642        domain.check_integrity()
1643
1644        #domain.visualise = True #If you want to take a sticky beak
1645
1646        initial_volume = domain.quantities['stage'].get_integral()
1647        initial_xmom = domain.quantities['xmomentum'].get_integral()
1648
1649        #print initial_xmom
1650
1651        #Evolution
1652        for t in domain.evolve(yieldstep = 0.05, finaltime = 5.0):
1653            volume =  domain.quantities['stage'].get_integral()
1654            assert allclose (volume, initial_volume)
1655
1656            #FIXME: What would we expect from momentum
1657            #xmom = domain.quantities['xmomentum'].get_integral()
1658            #print xmom
1659            #assert allclose (xmom, initial_xmom)
1660
1661
1662    def test_conservation_3(self):
1663        """Test that stage is conserved globally
1664
1665        This one uses a larger grid, convoluted bed, reflective bdries and a suitable
1666        initial condition
1667        """
1668        from mesh_factory import rectangular
1669        from shallow_water import Domain, Reflective_boundary,\
1670             Dirichlet_boundary, Constant_height
1671        from Numeric import array
1672
1673        #Create basic mesh
1674        points, vertices, boundary = rectangular(2, 1)
1675
1676        #Create shallow water domain
1677        domain = Domain(points, vertices, boundary)
1678        domain.smooth = False
1679        domain.default_order=2
1680        domain.beta_h = 0.2
1681
1682        #IC
1683        def x_slope(x, y):
1684            z = 0*x
1685            for i in range(len(x)):
1686                if x[i] < 0.3:
1687                    z[i] = x[i]/3
1688                if 0.3 <= x[i] < 0.5:
1689                    z[i] = -0.5
1690                if 0.5 <= x[i] < 0.7:
1691                    z[i] = 0.39
1692                if 0.7 <= x[i]:
1693                    z[i] = x[i]/3
1694            return z
1695
1696
1697
1698        domain.set_quantity('elevation', x_slope)
1699        domain.set_quantity('friction', 0)
1700        domain.set_quantity('stage', 0.4) #Steady
1701
1702        # Boundary conditions (reflective everywhere)
1703        Br = Reflective_boundary(domain)
1704        domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
1705
1706        domain.check_integrity()
1707
1708        #domain.visualise = True #If you want to take a sticky beak
1709
1710        initial_volume = domain.quantities['stage'].get_integral()
1711        initial_xmom = domain.quantities['xmomentum'].get_integral()
1712
1713        import copy
1714        ref_centroid_values =\
1715                 copy.copy(domain.quantities['stage'].centroid_values)
1716
1717        #print 'ORG', domain.quantities['stage'].centroid_values
1718        domain.distribute_to_vertices_and_edges()
1719
1720
1721        #print domain.quantities['stage'].centroid_values
1722        assert allclose(domain.quantities['stage'].centroid_values,
1723                        ref_centroid_values)
1724
1725
1726        #Check that initial limiter doesn't violate cons quan
1727        assert allclose (domain.quantities['stage'].get_integral(),
1728                         initial_volume)
1729
1730        #Evolution
1731        for t in domain.evolve(yieldstep = 0.05, finaltime = 10):
1732            volume =  domain.quantities['stage'].get_integral()
1733            #print t, volume, initial_volume
1734            assert allclose (volume, initial_volume)
1735
1736
1737    def test_conservation_4(self):
1738        """Test that stage is conserved globally
1739
1740        This one uses a larger grid, convoluted bed, reflective bdries and a suitable
1741        initial condition
1742        """
1743        from mesh_factory import rectangular
1744        from shallow_water import Domain, Reflective_boundary,\
1745             Dirichlet_boundary, Constant_height
1746        from Numeric import array
1747
1748        #Create basic mesh
1749        points, vertices, boundary = rectangular(6, 6)
1750
1751        #Create shallow water domain
1752        domain = Domain(points, vertices, boundary)
1753        domain.smooth = False
1754        domain.default_order=2
1755        domain.beta_h = 0.0
1756
1757        #IC
1758        def x_slope(x, y):
1759            z = 0*x
1760            for i in range(len(x)):
1761                if x[i] < 0.3:
1762                    z[i] = x[i]/3
1763                if 0.3 <= x[i] < 0.5:
1764                    z[i] = -0.5
1765                if 0.5 <= x[i] < 0.7:
1766                    #z[i] = 0.3 #OK with beta == 0.2
1767                    z[i] = 0.34 #OK with beta == 0.0
1768                    #z[i] = 0.35#Fails after 80 timesteps with an error
1769                                #of the order 1.0e-5
1770                if 0.7 <= x[i]:
1771                    z[i] = x[i]/3
1772            return z
1773
1774
1775
1776        domain.set_quantity('elevation', x_slope)
1777        domain.set_quantity('friction', 0)
1778        domain.set_quantity('stage', 0.4) #Steady
1779
1780        # Boundary conditions (reflective everywhere)
1781        Br = Reflective_boundary(domain)
1782        domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
1783
1784        domain.check_integrity()
1785
1786        #domain.visualise = True #If you want to take a sticky beak
1787
1788        initial_volume = domain.quantities['stage'].get_integral()
1789        initial_xmom = domain.quantities['xmomentum'].get_integral()
1790
1791        import copy
1792        ref_centroid_values =\
1793                 copy.copy(domain.quantities['stage'].centroid_values)
1794
1795        #Test limiter by itself
1796        domain.distribute_to_vertices_and_edges()
1797
1798        #Check that initial limiter doesn't violate cons quan
1799        assert allclose (domain.quantities['stage'].get_integral(),
1800                         initial_volume)
1801        #NOTE: This would fail if any initial stage was less than the
1802        #corresponding bed elevation - but that is reasonable.
1803
1804
1805        #Evolution
1806        for t in domain.evolve(yieldstep = 0.05, finaltime = 10.0):
1807            volume =  domain.quantities['stage'].get_integral()
1808
1809            #print t, volume, initial_volume
1810
1811
1812            #if not allclose (volume, initial_volume):
1813            #    print 't==4.05'
1814            #    for k in range(domain.number_of_elements):
1815            #       pass
1816            #       print domain.quantities['stage'].centroid_values[k] -\
1817            #             ref_centroid_values[k]
1818
1819            assert allclose (volume, initial_volume)
1820
1821
1822
1823
1824
1825    def test_conservation_5(self):
1826        """Test that momentum is conserved globally in
1827        steady state scenario
1828
1829        This one uses a slopy bed, dirichlet and reflective bdries
1830        """
1831        from mesh_factory import rectangular
1832        from shallow_water import Domain, Reflective_boundary,\
1833             Dirichlet_boundary, Constant_height
1834        from Numeric import array
1835
1836        #Create basic mesh
1837        points, vertices, boundary = rectangular(6, 6)
1838
1839        #Create shallow water domain
1840        domain = Domain(points, vertices, boundary)
1841        domain.smooth = False
1842        domain.default_order=2
1843
1844        #IC
1845        def x_slope(x, y):
1846            return x/3
1847
1848        domain.set_quantity('elevation', x_slope)
1849        domain.set_quantity('friction', 0)
1850        domain.set_quantity('stage', 0.4) #Steady
1851
1852        # Boundary conditions (reflective everywhere)
1853        Br = Reflective_boundary(domain)
1854        Bleft = Dirichlet_boundary([0.5,0,0])
1855        Bright = Dirichlet_boundary([0.1,0,0])
1856        domain.set_boundary({'left': Bleft, 'right': Bright,
1857                             'top': Br, 'bottom': Br})
1858
1859        domain.check_integrity()
1860
1861        #domain.visualise = True #If you want to take a sticky beak
1862
1863        initial_volume = domain.quantities['stage'].get_integral()
1864        initial_xmom = domain.quantities['xmomentum'].get_integral()
1865
1866
1867        #Evolution
1868        for t in domain.evolve(yieldstep = 0.05, finaltime = 15.0):
1869            stage =  domain.quantities['stage'].get_integral()
1870            xmom = domain.quantities['xmomentum'].get_integral()
1871            ymom = domain.quantities['ymomentum'].get_integral()
1872
1873            if allclose(t, 6):  #Steady state reached
1874                steady_xmom = domain.quantities['xmomentum'].get_integral()
1875                steady_ymom = domain.quantities['ymomentum'].get_integral()
1876                steady_stage = domain.quantities['stage'].get_integral()
1877
1878            if t > 6:
1879                #print '%.2f %14.8f %14.8f' %(t, ymom, steady_ymom)
1880                assert allclose(xmom, steady_xmom)
1881                assert allclose(ymom, steady_ymom)
1882                assert allclose(stage, steady_stage)
1883
1884
1885
1886
1887
1888
1889
1890
1891    def test_conservation_real(self):
1892        """Test that momentum is conserved globally
1893       
1894        Stephen finally made a test that revealed the problem.
1895        This test failed with code prior to 25 July 2005
1896        """
1897
1898        yieldstep = 0.01
1899        finaltime = 0.05
1900        min_depth = 1.0e-2
1901
1902
1903        import sys
1904        from os import sep; sys.path.append('..'+sep+'pyvolution')
1905        from mesh_factory import rectangular
1906        from shallow_water import Domain, Reflective_boundary
1907
1908
1909        #Create shallow water domain
1910        points, vertices, boundary = rectangular(10, 10, len1=500, len2=500)
1911        domain = Domain(points, vertices, boundary)
1912        domain.smooth = False
1913        domain.visualise = False
1914        domain.default_order = 1
1915        domain.minimum_allowed_height = min_depth
1916
1917        # Set initial condition
1918        class Set_IC:
1919            """Set an initial condition with a constant value, for x0<x<x1
1920            """
1921
1922            def __init__(self, x0=0.25, x1=0.5, h=1.0):
1923                self.x0 = x0
1924                self.x1 = x1
1925                self.= h
1926
1927            def __call__(self, x, y):
1928                return self.h*((x>self.x0)&(x<self.x1))
1929
1930
1931        domain.set_quantity('stage', Set_IC(200.0,300.0,5.0))
1932
1933
1934        #Boundaries
1935        R = Reflective_boundary(domain)
1936        domain.set_boundary( {'left': R, 'right': R, 'top':R, 'bottom': R})
1937
1938        ref = domain.quantities['stage'].get_integral()
1939
1940        # Evolution
1941        for t in domain.evolve(yieldstep = yieldstep, finaltime = finaltime):
1942            pass
1943            #print 'Integral stage = ',\
1944            #      domain.quantities['stage'].get_integral(),\
1945            #      ' Time = ',domain.time
1946
1947
1948        now = domain.quantities['stage'].get_integral()
1949
1950        msg = 'Stage not conserved: was %f, now %f' %(ref, now)
1951        assert allclose(ref, now), msg
1952
1953
1954
1955    def test_second_order_flat_bed_onestep(self):
1956
1957        from mesh_factory import rectangular
1958        from shallow_water import Domain, Reflective_boundary,\
1959             Dirichlet_boundary, Constant_height
1960        from Numeric import array
1961
1962        #Create basic mesh
1963        points, vertices, boundary = rectangular(6, 6)
1964
1965        #Create shallow water domain
1966        domain = Domain(points, vertices, boundary)
1967        domain.smooth = False
1968        domain.default_order=2
1969
1970        # Boundary conditions
1971        Br = Reflective_boundary(domain)
1972        Bd = Dirichlet_boundary([0.1, 0., 0.])
1973        domain.set_boundary({'left': Bd, 'right': Br, 'top': Br, 'bottom': Br})
1974
1975        domain.check_integrity()
1976
1977        #Evolution
1978        for t in domain.evolve(yieldstep = 0.05, finaltime = 0.05):
1979            pass# domain.write_time()
1980
1981        #Data from earlier version of pyvolution
1982        assert allclose(domain.min_timestep, 0.0396825396825)
1983        assert allclose(domain.max_timestep, 0.0396825396825)
1984
1985        assert allclose(domain.quantities['stage'].centroid_values[:12],
1986                        [0.00171396, 0.02656103, 0.00241523, 0.02656103,
1987                        0.00241523, 0.02656103, 0.00241523, 0.02656103,
1988                        0.00241523, 0.02656103, 0.00241523, 0.0272623])
1989
1990        domain.distribute_to_vertices_and_edges()
1991        assert allclose(domain.quantities['stage'].vertex_values[:12,0],
1992                        [0.0001714, 0.02656103, 0.00024152,
1993                        0.02656103, 0.00024152, 0.02656103,
1994                        0.00024152, 0.02656103, 0.00024152,
1995                        0.02656103, 0.00024152, 0.0272623])
1996
1997        assert allclose(domain.quantities['stage'].vertex_values[:12,1],
1998                        [0.00315012, 0.02656103, 0.00024152, 0.02656103,
1999                         0.00024152, 0.02656103, 0.00024152, 0.02656103,
2000                         0.00024152, 0.02656103, 0.00040506, 0.0272623])
2001
2002        assert allclose(domain.quantities['stage'].vertex_values[:12,2],
2003                        [0.00182037, 0.02656103, 0.00676264,
2004                         0.02656103, 0.00676264, 0.02656103,
2005                         0.00676264, 0.02656103, 0.00676264,
2006                         0.02656103, 0.0065991, 0.0272623])
2007
2008        assert allclose(domain.quantities['xmomentum'].centroid_values[:12],
2009                        [0.00113961, 0.01302432, 0.00148672,
2010                         0.01302432, 0.00148672, 0.01302432,
2011                         0.00148672, 0.01302432, 0.00148672 ,
2012                         0.01302432, 0.00148672, 0.01337143])
2013
2014        assert allclose(domain.quantities['ymomentum'].centroid_values[:12],
2015                        [-2.91240050e-004 , 1.22721531e-004,
2016                         -1.22721531e-004, 1.22721531e-004 ,
2017                         -1.22721531e-004, 1.22721531e-004,
2018                         -1.22721531e-004 , 1.22721531e-004,
2019                         -1.22721531e-004, 1.22721531e-004,
2020                         -1.22721531e-004, -4.57969873e-005])
2021
2022
2023
2024    def test_second_order_flat_bed_moresteps(self):
2025
2026        from mesh_factory import rectangular
2027        from shallow_water import Domain, Reflective_boundary,\
2028             Dirichlet_boundary, Constant_height
2029        from Numeric import array
2030
2031        #Create basic mesh
2032        points, vertices, boundary = rectangular(6, 6)
2033
2034        #Create shallow water domain
2035        domain = Domain(points, vertices, boundary)
2036        domain.smooth = False
2037        domain.default_order=2
2038
2039        # Boundary conditions
2040        Br = Reflective_boundary(domain)
2041        Bd = Dirichlet_boundary([0.1, 0., 0.])
2042        domain.set_boundary({'left': Bd, 'right': Br, 'top': Br, 'bottom': Br})
2043
2044        domain.check_integrity()
2045
2046        #Evolution
2047        for t in domain.evolve(yieldstep = 0.05, finaltime = 0.1):
2048            pass
2049
2050        #Data from earlier version of pyvolution
2051        #assert allclose(domain.min_timestep, 0.0396825396825)
2052        #assert allclose(domain.max_timestep, 0.0396825396825)
2053        #print domain.quantities['stage'].centroid_values
2054
2055
2056    def test_flatbed_first_order(self):
2057        from mesh_factory import rectangular
2058        from shallow_water import Domain,\
2059             Reflective_boundary, Dirichlet_boundary
2060
2061        from Numeric import array
2062
2063        #Create basic mesh
2064        N = 8
2065        points, vertices, boundary = rectangular(N, N)
2066
2067        #Create shallow water domain
2068        domain = Domain(points, vertices, boundary)
2069        domain.smooth = False
2070        domain.visualise = False
2071        domain.default_order=1
2072
2073        # Boundary conditions
2074        Br = Reflective_boundary(domain)
2075        Bd = Dirichlet_boundary([0.2,0.,0.])
2076
2077        domain.set_boundary({'left': Bd, 'right': Br, 'top': Br, 'bottom': Br})
2078        domain.check_integrity()
2079
2080
2081        #Evolution
2082        for t in domain.evolve(yieldstep = 0.02, finaltime = 0.5):
2083            pass
2084            #domain.write_time()
2085
2086        #FIXME: These numbers were from version before 25/10
2087        #assert allclose(domain.min_timestep, 0.0140413643926)
2088        #assert allclose(domain.max_timestep, 0.0140947355753)
2089
2090        for i in range(3):
2091            #assert allclose(domain.quantities['stage'].edge_values[:4,i],
2092            #                [0.10730244,0.12337617,0.11200126,0.12605666])
2093
2094            assert allclose(domain.quantities['xmomentum'].edge_values[:4,i],
2095                            [0.07610894,0.06901572,0.07284461,0.06819712])
2096
2097            #assert allclose(domain.quantities['ymomentum'].edge_values[:4,i],
2098            #                [-0.0060238, -0.00157404, -0.00309633, -0.0001637])
2099
2100    def test_flatbed_second_order(self):
2101        from mesh_factory import rectangular
2102        from shallow_water import Domain,\
2103             Reflective_boundary, Dirichlet_boundary
2104
2105        from Numeric import array
2106
2107        #Create basic mesh
2108        N = 8
2109        points, vertices, boundary = rectangular(N, N)
2110
2111        #Create shallow water domain
2112        domain = Domain(points, vertices, boundary)
2113        domain.smooth = False
2114        domain.visualise = False
2115        domain.default_order=2
2116        #domain.minimum_allowed_height = 0.0 #Makes it like the 'oldstyle' balance
2117
2118        # Boundary conditions
2119        Br = Reflective_boundary(domain)
2120        Bd = Dirichlet_boundary([0.2,0.,0.])
2121
2122        domain.set_boundary({'left': Bd, 'right': Br, 'top': Br, 'bottom': Br})
2123        domain.check_integrity()
2124
2125        #Evolution
2126        for t in domain.evolve(yieldstep = 0.01, finaltime = 0.03):
2127            pass
2128
2129
2130        assert allclose(domain.min_timestep, 0.0210448446782)
2131        assert allclose(domain.max_timestep, 0.0210448446782)
2132
2133        #print domain.quantities['stage'].vertex_values[:4,0]
2134        #print domain.quantities['xmomentum'].vertex_values[:4,0]
2135        #print domain.quantities['ymomentum'].vertex_values[:4,0]
2136
2137        #FIXME: These numbers were from version before 25/10
2138        #assert allclose(domain.quantities['stage'].vertex_values[:4,0],
2139        #                [0.00101913,0.05352143,0.00104852,0.05354394])
2140
2141        assert allclose(domain.quantities['xmomentum'].vertex_values[:4,0],
2142                        [ 0.00064835, 0.03685719, 0.00085073, 0.03687313])
2143
2144        #assert allclose(domain.quantities['xmomentum'].vertex_values[:4,0],
2145        #                [0.00090581,0.03685719,0.00088303,0.03687313])
2146
2147        assert allclose(domain.quantities['ymomentum'].vertex_values[:4,0],
2148                        [-0.00139463,0.0006156,-0.00060364,0.00061827])
2149
2150
2151
2152    def test_flatbed_second_order_distribute(self):
2153        #Use real data from pyvolution 2
2154        #painfully setup and extracted.
2155        from mesh_factory import rectangular
2156        from shallow_water import Domain,\
2157             Reflective_boundary, Dirichlet_boundary
2158
2159        from Numeric import array
2160
2161        #Create basic mesh
2162        N = 8
2163        points, vertices, boundary = rectangular(N, N)
2164
2165        #Create shallow water domain
2166        domain = Domain(points, vertices, boundary)
2167        domain.smooth = False
2168        domain.visualise = False
2169        domain.default_order=domain.order=2
2170
2171        # Boundary conditions
2172        Br = Reflective_boundary(domain)
2173        Bd = Dirichlet_boundary([0.2,0.,0.])
2174
2175        domain.set_boundary({'left': Bd, 'right': Br, 'top': Br, 'bottom': Br})
2176        domain.check_integrity()
2177
2178
2179
2180        for V in [False, True]:
2181            if V:
2182                #Set centroids as if system had been evolved
2183                L = zeros(2*N*N, Float)
2184                L[:32] = [7.21205592e-003, 5.35214298e-002, 1.00910824e-002,
2185                          5.35439433e-002, 1.00910824e-002, 5.35439433e-002,
2186                          1.00910824e-002, 5.35439433e-002, 1.00910824e-002,
2187                          5.35439433e-002, 1.00910824e-002, 5.35439433e-002,
2188                          1.00910824e-002, 5.35393928e-002, 1.02344264e-002,
2189                          5.59605058e-002, 0.00000000e+000, 3.31027800e-004,
2190                          0.00000000e+000, 4.37962142e-005, 0.00000000e+000,
2191                          4.37962142e-005, 0.00000000e+000, 4.37962142e-005,
2192                          0.00000000e+000, 4.37962142e-005, 0.00000000e+000,
2193                          4.37962142e-005, 0.00000000e+000, 4.37962142e-005,
2194                          0.00000000e+000, 5.57305948e-005]
2195
2196                X = zeros(2*N*N, Float)
2197                X[:32] = [6.48351607e-003, 3.68571894e-002, 8.50733285e-003,
2198                          3.68731327e-002, 8.50733285e-003, 3.68731327e-002,
2199                          8.50733285e-003, 3.68731327e-002, 8.50733285e-003,
2200                          3.68731327e-002, 8.50733285e-003, 3.68731327e-002,
2201                          8.50733285e-003, 3.68693861e-002, 8.65220973e-003,
2202                          3.85055387e-002, 0.00000000e+000, 2.86060840e-004,
2203                          0.00000000e+000, 3.58905503e-005, 0.00000000e+000,
2204                          3.58905503e-005, 0.00000000e+000, 3.58905503e-005,
2205                          0.00000000e+000, 3.58905503e-005, 0.00000000e+000,
2206                          3.58905503e-005, 0.00000000e+000, 3.58905503e-005,
2207                          0.00000000e+000, 4.57662812e-005]
2208
2209                Y = zeros(2*N*N, Float)
2210                Y[:32]=[-1.39463104e-003, 6.15600298e-004, -6.03637382e-004,
2211                        6.18272251e-004, -6.03637382e-004, 6.18272251e-004,
2212                        -6.03637382e-004, 6.18272251e-004, -6.03637382e-004,
2213                        6.18272251e-004, -6.03637382e-004, 6.18272251e-004,
2214                        -6.03637382e-004, 6.18599320e-004, -6.74622797e-004,
2215                        -1.48934756e-004, 0.00000000e+000, -5.35079969e-005,
2216                        0.00000000e+000, -2.57264987e-005, 0.00000000e+000,
2217                        -2.57264987e-005, 0.00000000e+000, -2.57264987e-005,
2218                        0.00000000e+000, -2.57264987e-005, 0.00000000e+000,
2219                        -2.57264987e-005, 0.00000000e+000, -2.57264987e-005,
2220                        0.00000000e+000, -2.57635178e-005]
2221
2222
2223                domain.set_quantity('stage', L, location='centroids')
2224                domain.set_quantity('xmomentum', X, location='centroids')
2225                domain.set_quantity('ymomentum', Y, location='centroids')
2226
2227                domain.check_integrity()
2228            else:
2229                #Evolution
2230                for t in domain.evolve(yieldstep = 0.01, finaltime = 0.03):
2231                    pass
2232                assert allclose(domain.min_timestep, 0.0210448446782)
2233                assert allclose(domain.max_timestep, 0.0210448446782)
2234
2235
2236            #Centroids were correct but not vertices.
2237            #Hence the check of distribute below.
2238            assert allclose(domain.quantities['stage'].centroid_values[:4],
2239                            [0.00721206,0.05352143,0.01009108,0.05354394])
2240
2241            assert allclose(domain.quantities['xmomentum'].centroid_values[:4],
2242                            [0.00648352,0.03685719,0.00850733,0.03687313])
2243
2244            assert allclose(domain.quantities['ymomentum'].centroid_values[:4],
2245                            [-0.00139463,0.0006156,-0.00060364,0.00061827])
2246
2247            #print 'C17=', domain.quantities['xmomentum'].centroid_values[17]
2248            #print 'C19=', domain.quantities['xmomentum'].centroid_values[19]
2249
2250            #assert allclose(domain.quantities['xmomentum'].centroid_values[17],0.00028606084)
2251            ##print domain.quantities['xmomentum'].centroid_values[17], V
2252            ##print
2253            if not V:
2254                assert allclose(domain.quantities['xmomentum'].centroid_values[17], 0.0)
2255            else:
2256                assert allclose(domain.quantities['xmomentum'].centroid_values[17], 0.00028606084)
2257
2258            import copy
2259            XX = copy.copy(domain.quantities['xmomentum'].centroid_values)
2260            assert allclose(domain.quantities['xmomentum'].centroid_values, XX)
2261
2262            domain.distribute_to_vertices_and_edges()
2263
2264            #assert allclose(domain.quantities['xmomentum'].centroid_values, XX)
2265
2266            assert allclose(domain.quantities['xmomentum'].centroid_values[17],
2267                            0.0)
2268
2269
2270            #FIXME: These numbers were from version before 25/10
2271            #assert allclose(domain.quantities['stage'].vertex_values[:4,0],
2272            #                [0.00101913,0.05352143,0.00104852,0.05354394])
2273
2274            assert allclose(domain.quantities['ymomentum'].vertex_values[:4,0],
2275                            [-0.00139463,0.0006156,-0.00060364,0.00061827])
2276
2277
2278            assert allclose(domain.quantities['xmomentum'].vertex_values[:4,0],
2279                            [0.00064835,0.03685719,0.00085073,0.03687313])
2280
2281
2282            #NB NO longer relvant:
2283
2284            #This was the culprit. First triangles vertex 0 had an
2285            #x-momentum of 0.0064835 instead of 0.00090581 and
2286            #third triangle had 0.00850733 instead of 0.00088303
2287            #print domain.quantities['xmomentum'].vertex_values[:4,0]
2288
2289            #print domain.quantities['xmomentum'].vertex_values[:4,0]
2290            #assert allclose(domain.quantities['xmomentum'].vertex_values[:4,0],
2291            #                [0.00090581,0.03685719,0.00088303,0.03687313])
2292
2293
2294
2295
2296
2297    def test_bedslope_problem_first_order(self):
2298
2299        from mesh_factory import rectangular
2300        from shallow_water import Domain, Reflective_boundary, Constant_height
2301        from Numeric import array
2302
2303        #Create basic mesh
2304        points, vertices, boundary = rectangular(6, 6)
2305
2306        #Create shallow water domain
2307        domain = Domain(points, vertices, boundary)
2308        domain.smooth = False
2309        domain.default_order=1
2310
2311        #Bed-slope and friction
2312        def x_slope(x, y):
2313            return -x/3
2314
2315        domain.set_quantity('elevation', x_slope)
2316
2317        # Boundary conditions
2318        Br = Reflective_boundary(domain)
2319        domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
2320
2321        #Initial condition
2322        domain.set_quantity('stage', Constant_height(x_slope, 0.05))
2323        domain.check_integrity()
2324
2325        #Evolution
2326        for t in domain.evolve(yieldstep = 0.05, finaltime = 0.05):
2327            pass# domain.write_time()
2328
2329        assert allclose(domain.min_timestep, 0.050010003001)
2330        assert allclose(domain.max_timestep, 0.050010003001)
2331
2332
2333    def test_bedslope_problem_first_order_moresteps(self):
2334
2335        from mesh_factory import rectangular
2336        from shallow_water import Domain, Reflective_boundary, Constant_height
2337        from Numeric import array
2338
2339        #Create basic mesh
2340        points, vertices, boundary = rectangular(6, 6)
2341
2342        #Create shallow water domain
2343        domain = Domain(points, vertices, boundary)
2344        domain.smooth = False
2345        domain.default_order=1
2346        domain.beta_h = 0.0 #Use first order in h-limiter
2347
2348        #Bed-slope and friction
2349        def x_slope(x, y):
2350            return -x/3
2351
2352        domain.set_quantity('elevation', x_slope)
2353
2354        # Boundary conditions
2355        Br = Reflective_boundary(domain)
2356        domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
2357
2358        #Initial condition
2359        domain.set_quantity('stage', Constant_height(x_slope, 0.05))
2360        domain.check_integrity()
2361
2362        #Evolution
2363        for t in domain.evolve(yieldstep = 0.05, finaltime = 0.5):
2364            pass# domain.write_time()
2365
2366        #Data from earlier version of pyvolution
2367        #print domain.quantities['stage'].centroid_values
2368
2369        assert allclose(domain.quantities['stage'].centroid_values,
2370                        [-0.02998628, -0.01520652, -0.03043492,
2371                        -0.0149132, -0.03004706, -0.01476251,
2372                        -0.0298215, -0.01467976, -0.02988158,
2373                        -0.01474662, -0.03036161, -0.01442995,
2374                        -0.07624583, -0.06297061, -0.07733792,
2375                        -0.06342237, -0.07695439, -0.06289595,
2376                        -0.07635559, -0.0626065, -0.07633628,
2377                        -0.06280072, -0.07739632, -0.06386738,
2378                        -0.12161738, -0.11028239, -0.1223796,
2379                        -0.11095953, -0.12189744, -0.11048616,
2380                        -0.12074535, -0.10987605, -0.12014311,
2381                        -0.10976691, -0.12096859, -0.11087692,
2382                        -0.16868259, -0.15868061, -0.16801135,
2383                        -0.1588003, -0.16674343, -0.15813323,
2384                        -0.16457595, -0.15693826, -0.16281096,
2385                        -0.15585154, -0.16283873, -0.15540068,
2386                        -0.17450362, -0.19919913, -0.18062882,
2387                        -0.19764131, -0.17783111, -0.19407213,
2388                        -0.1736915, -0.19053624, -0.17228678,
2389                        -0.19105634, -0.17920133, -0.1968828,
2390                        -0.14244395, -0.14604641, -0.14473537,
2391                        -0.1506107, -0.14510055, -0.14919522,
2392                        -0.14175896, -0.14560798, -0.13911658,
2393                        -0.14439383, -0.13924047, -0.14829043])
2394
2395
2396    def test_bedslope_problem_second_order_one_step(self):
2397
2398        from mesh_factory import rectangular
2399        from shallow_water import Domain, Reflective_boundary, Constant_height
2400        from Numeric import array
2401
2402        #Create basic mesh
2403        points, vertices, boundary = rectangular(6, 6)
2404
2405        #Create shallow water domain
2406        domain = Domain(points, vertices, boundary)
2407        domain.smooth = False
2408        domain.default_order=2
2409
2410        #Bed-slope and friction at vertices (and interpolated elsewhere)
2411        def x_slope(x, y):
2412            return -x/3
2413
2414        domain.set_quantity('elevation', x_slope)
2415
2416        # Boundary conditions
2417        Br = Reflective_boundary(domain)
2418        domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
2419
2420        #Initial condition
2421        domain.set_quantity('stage', Constant_height(x_slope, 0.05))
2422        domain.check_integrity()
2423
2424        assert allclose(domain.quantities['stage'].centroid_values,
2425                        [0.01296296, 0.03148148, 0.01296296,
2426                        0.03148148, 0.01296296, 0.03148148,
2427                        0.01296296, 0.03148148, 0.01296296,
2428                        0.03148148, 0.01296296, 0.03148148,
2429                        -0.04259259, -0.02407407, -0.04259259,
2430                        -0.02407407, -0.04259259, -0.02407407,
2431                        -0.04259259, -0.02407407, -0.04259259,
2432                        -0.02407407, -0.04259259, -0.02407407,
2433                        -0.09814815, -0.07962963, -0.09814815,
2434                        -0.07962963, -0.09814815, -0.07962963,
2435                        -0.09814815, -0.07962963, -0.09814815,
2436                        -0.07962963, -0.09814815, -0.07962963,
2437                        -0.1537037 , -0.13518519, -0.1537037,
2438                        -0.13518519, -0.1537037, -0.13518519,
2439                        -0.1537037 , -0.13518519, -0.1537037,
2440                        -0.13518519, -0.1537037, -0.13518519,
2441                        -0.20925926, -0.19074074, -0.20925926,
2442                        -0.19074074, -0.20925926, -0.19074074,
2443                        -0.20925926, -0.19074074, -0.20925926,
2444                        -0.19074074, -0.20925926, -0.19074074,
2445                        -0.26481481, -0.2462963, -0.26481481,
2446                        -0.2462963, -0.26481481, -0.2462963,
2447                        -0.26481481, -0.2462963, -0.26481481,
2448                        -0.2462963, -0.26481481, -0.2462963])
2449
2450
2451        #print domain.quantities['stage'].extrapolate_second_order()
2452        #domain.distribute_to_vertices_and_edges()
2453        #print domain.quantities['stage'].vertex_values[:,0]
2454
2455        #Evolution
2456        for t in domain.evolve(yieldstep = 0.05, finaltime = 0.05):
2457            #domain.write_time()
2458            pass
2459
2460
2461        #print domain.quantities['stage'].centroid_values
2462        assert allclose(domain.quantities['stage'].centroid_values,
2463                 [0.01290985, 0.02356019, 0.01619096, 0.02356019, 0.01619096,
2464                  0.02356019, 0.01619096, 0.02356019, 0.01619096, 0.02356019,
2465                  0.01619096, 0.0268413, -0.04411074, -0.0248011, -0.04186556,
2466                  -0.0248011, -0.04186556, -0.0248011, -0.04186556, -0.0248011,
2467                  -0.04186556, -0.0248011, -0.04186556, -0.02255593,
2468                  -0.09966629, -0.08035666, -0.09742112, -0.08035666,
2469                  -0.09742112, -0.08035666, -0.09742112, -0.08035666,
2470                  -0.09742112, -0.08035666, -0.09742112, -0.07811149,
2471                  -0.15522185, -0.13591222, -0.15297667, -0.13591222,
2472                  -0.15297667, -0.13591222, -0.15297667, -0.13591222,
2473                  -0.15297667, -0.13591222, -0.15297667, -0.13366704,
2474                  -0.2107774, -0.19146777, -0.20853223, -0.19146777,
2475                  -0.20853223, -0.19146777, -0.20853223, -0.19146777,
2476                  -0.20853223, -0.19146777, -0.20853223, -0.1892226,
2477                  -0.26120669, -0.24776246, -0.25865535, -0.24776246,
2478                  -0.25865535, -0.24776246, -0.25865535, -0.24776246,
2479                  -0.25865535, -0.24776246, -0.25865535, -0.24521113])
2480
2481
2482
2483    def test_bedslope_problem_second_order_two_steps(self):
2484
2485        from mesh_factory import rectangular
2486        from shallow_water import Domain, Reflective_boundary, Constant_height
2487        from Numeric import array
2488
2489        #Create basic mesh
2490        points, vertices, boundary = rectangular(6, 6)
2491
2492        #Create shallow water domain
2493        domain = Domain(points, vertices, boundary)
2494        domain.smooth = False
2495        domain.default_order=2
2496        domain.beta_h = 0.0 #Use first order in h-limiter
2497
2498        #Bed-slope and friction at vertices (and interpolated elsewhere)
2499        def x_slope(x, y):
2500            return -x/3
2501
2502        domain.set_quantity('elevation', x_slope)
2503
2504        # Boundary conditions
2505        Br = Reflective_boundary(domain)
2506        domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
2507
2508        #Initial condition
2509        domain.set_quantity('stage', Constant_height(x_slope, 0.05))
2510        domain.check_integrity()
2511
2512        assert allclose(domain.quantities['stage'].centroid_values,
2513                        [0.01296296, 0.03148148, 0.01296296,
2514                        0.03148148, 0.01296296, 0.03148148,
2515                        0.01296296, 0.03148148, 0.01296296,
2516                        0.03148148, 0.01296296, 0.03148148,
2517                        -0.04259259, -0.02407407, -0.04259259,
2518                        -0.02407407, -0.04259259, -0.02407407,
2519                        -0.04259259, -0.02407407, -0.04259259,
2520                        -0.02407407, -0.04259259, -0.02407407,
2521                        -0.09814815, -0.07962963, -0.09814815,
2522                        -0.07962963, -0.09814815, -0.07962963,
2523                        -0.09814815, -0.07962963, -0.09814815,
2524                        -0.07962963, -0.09814815, -0.07962963,
2525                        -0.1537037 , -0.13518519, -0.1537037,
2526                        -0.13518519, -0.1537037, -0.13518519,
2527                        -0.1537037 , -0.13518519, -0.1537037,
2528                        -0.13518519, -0.1537037, -0.13518519,
2529                        -0.20925926, -0.19074074, -0.20925926,
2530                        -0.19074074, -0.20925926, -0.19074074,
2531                        -0.20925926, -0.19074074, -0.20925926,
2532                        -0.19074074, -0.20925926, -0.19074074,
2533                        -0.26481481, -0.2462963, -0.26481481,
2534                        -0.2462963, -0.26481481, -0.2462963,
2535                        -0.26481481, -0.2462963, -0.26481481,
2536                        -0.2462963, -0.26481481, -0.2462963])
2537
2538
2539        #print domain.quantities['stage'].extrapolate_second_order()
2540        #domain.distribute_to_vertices_and_edges()
2541        #print domain.quantities['stage'].vertex_values[:,0]
2542
2543        #Evolution
2544        for t in domain.evolve(yieldstep = 0.05, finaltime = 0.1):
2545            pass
2546
2547
2548        #Data from earlier version of pyvolution ft=0.1
2549        assert allclose(domain.min_timestep, 0.0376895634803)
2550        assert allclose(domain.max_timestep, 0.0415635655309)
2551
2552
2553        assert allclose(domain.quantities['stage'].centroid_values,
2554                        [0.00855788, 0.01575204, 0.00994606, 0.01717072,
2555                         0.01005985, 0.01716362, 0.01005985, 0.01716299,
2556                         0.01007098, 0.01736248, 0.01216452, 0.02026776,
2557                         -0.04462374, -0.02479045, -0.04199789, -0.0229465,
2558                         -0.04184033, -0.02295693, -0.04184013, -0.02295675,
2559                         -0.04184486, -0.0228168, -0.04028876, -0.02036486,
2560                         -0.10029444, -0.08170809, -0.09772846, -0.08021704,
2561                         -0.09760006, -0.08022143, -0.09759984, -0.08022124,
2562                         -0.09760261, -0.08008893, -0.09603914, -0.07758209,
2563                         -0.15584152, -0.13723138, -0.15327266, -0.13572906,
2564                         -0.15314427, -0.13573349, -0.15314405, -0.13573331,
2565                         -0.15314679, -0.13560104, -0.15158523, -0.13310701,
2566                         -0.21208605, -0.19283913, -0.20955631, -0.19134189,
2567                         -0.20942821, -0.19134598, -0.20942799, -0.1913458,
2568                         -0.20943005, -0.19120952, -0.20781177, -0.18869401,
2569                         -0.25384082, -0.2463294, -0.25047649, -0.24464654,
2570                         -0.25031159, -0.24464253, -0.25031112, -0.24464253,
2571                         -0.25031463, -0.24454764, -0.24885323, -0.24286438])
2572
2573
2574
2575
2576    def test_bedslope_problem_second_order_two_yieldsteps(self):
2577
2578        from mesh_factory import rectangular
2579        from shallow_water import Domain, Reflective_boundary, Constant_height
2580        from Numeric import array
2581
2582        #Create basic mesh
2583        points, vertices, boundary = rectangular(6, 6)
2584
2585        #Create shallow water domain
2586        domain = Domain(points, vertices, boundary)
2587        domain.smooth = False
2588        domain.default_order=2
2589        domain.beta_h = 0.0 #Use first order in h-limiter
2590
2591        #Bed-slope and friction at vertices (and interpolated elsewhere)
2592        def x_slope(x, y):
2593            return -x/3
2594
2595        domain.set_quantity('elevation', x_slope)
2596
2597        # Boundary conditions
2598        Br = Reflective_boundary(domain)
2599        domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
2600
2601        #Initial condition
2602        domain.set_quantity('stage', Constant_height(x_slope, 0.05))
2603        domain.check_integrity()
2604
2605        assert allclose(domain.quantities['stage'].centroid_values,
2606                        [0.01296296, 0.03148148, 0.01296296,
2607                        0.03148148, 0.01296296, 0.03148148,
2608                        0.01296296, 0.03148148, 0.01296296,
2609                        0.03148148, 0.01296296, 0.03148148,
2610                        -0.04259259, -0.02407407, -0.04259259,
2611                        -0.02407407, -0.04259259, -0.02407407,
2612                        -0.04259259, -0.02407407, -0.04259259,
2613                        -0.02407407, -0.04259259, -0.02407407,
2614                        -0.09814815, -0.07962963, -0.09814815,
2615                        -0.07962963, -0.09814815, -0.07962963,
2616                        -0.09814815, -0.07962963, -0.09814815,
2617                        -0.07962963, -0.09814815, -0.07962963,
2618                        -0.1537037 , -0.13518519, -0.1537037,
2619                        -0.13518519, -0.1537037, -0.13518519,
2620                        -0.1537037 , -0.13518519, -0.1537037,
2621                        -0.13518519, -0.1537037, -0.13518519,
2622                        -0.20925926, -0.19074074, -0.20925926,
2623                        -0.19074074, -0.20925926, -0.19074074,
2624                        -0.20925926, -0.19074074, -0.20925926,
2625                        -0.19074074, -0.20925926, -0.19074074,
2626                        -0.26481481, -0.2462963, -0.26481481,
2627                        -0.2462963, -0.26481481, -0.2462963,
2628                        -0.26481481, -0.2462963, -0.26481481,
2629                        -0.2462963, -0.26481481, -0.2462963])
2630
2631
2632        #print domain.quantities['stage'].extrapolate_second_order()
2633        #domain.distribute_to_vertices_and_edges()
2634        #print domain.quantities['stage'].vertex_values[:,0]
2635
2636        #Evolution
2637        for t in domain.evolve(yieldstep = 0.05, finaltime = 0.1):   #0.05??
2638            #domain.write_time()
2639            pass
2640
2641
2642
2643        assert allclose(domain.quantities['stage'].centroid_values,
2644                 [0.00855788, 0.01575204, 0.00994606, 0.01717072, 0.01005985,
2645                  0.01716362, 0.01005985, 0.01716299, 0.01007098, 0.01736248,
2646                  0.01216452, 0.02026776, -0.04462374, -0.02479045, -0.04199789,
2647                  -0.0229465, -0.04184033, -0.02295693, -0.04184013,
2648                  -0.02295675, -0.04184486, -0.0228168, -0.04028876,
2649                  -0.02036486, -0.10029444, -0.08170809, -0.09772846,
2650                  -0.08021704, -0.09760006, -0.08022143, -0.09759984,
2651                  -0.08022124, -0.09760261, -0.08008893, -0.09603914,
2652                  -0.07758209, -0.15584152, -0.13723138, -0.15327266,
2653                  -0.13572906, -0.15314427, -0.13573349, -0.15314405,
2654                  -0.13573331, -0.15314679, -0.13560104, -0.15158523,
2655                  -0.13310701, -0.21208605, -0.19283913, -0.20955631,
2656                  -0.19134189, -0.20942821, -0.19134598, -0.20942799,
2657                  -0.1913458, -0.20943005, -0.19120952, -0.20781177,
2658                  -0.18869401, -0.25384082, -0.2463294, -0.25047649,
2659                  -0.24464654, -0.25031159, -0.24464253, -0.25031112,
2660                  -0.24464253, -0.25031463, -0.24454764, -0.24885323,
2661                  -0.24286438])
2662
2663
2664
2665    def test_bedslope_problem_second_order_more_steps(self):
2666
2667        from mesh_factory import rectangular
2668        from shallow_water import Domain, Reflective_boundary, Constant_height
2669        from Numeric import array
2670
2671        #Create basic mesh
2672        points, vertices, boundary = rectangular(6, 6)
2673
2674        #Create shallow water domain
2675        domain = Domain(points, vertices, boundary)
2676        domain.smooth = False
2677        domain.default_order=2
2678        domain.beta_h = 0.0 #Use first order in h-limiter
2679
2680        #Bed-slope and friction at vertices (and interpolated elsewhere)
2681        def x_slope(x, y):
2682            return -x/3
2683
2684        domain.set_quantity('elevation', x_slope)
2685
2686        # Boundary conditions
2687        Br = Reflective_boundary(domain)
2688        domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
2689
2690        #Initial condition
2691        domain.set_quantity('stage', Constant_height(x_slope, 0.05))
2692        domain.check_integrity()
2693
2694        assert allclose(domain.quantities['stage'].centroid_values,
2695                        [0.01296296, 0.03148148, 0.01296296,
2696                        0.03148148, 0.01296296, 0.03148148,
2697                        0.01296296, 0.03148148, 0.01296296,
2698                        0.03148148, 0.01296296, 0.03148148,
2699                        -0.04259259, -0.02407407, -0.04259259,
2700                        -0.02407407, -0.04259259, -0.02407407,
2701                        -0.04259259, -0.02407407, -0.04259259,
2702                        -0.02407407, -0.04259259, -0.02407407,
2703                        -0.09814815, -0.07962963, -0.09814815,
2704                        -0.07962963, -0.09814815, -0.07962963,
2705                        -0.09814815, -0.07962963, -0.09814815,
2706                        -0.07962963, -0.09814815, -0.07962963,
2707                        -0.1537037 , -0.13518519, -0.1537037,
2708                        -0.13518519, -0.1537037, -0.13518519,
2709                        -0.1537037 , -0.13518519, -0.1537037,
2710                        -0.13518519, -0.1537037, -0.13518519,
2711                        -0.20925926, -0.19074074, -0.20925926,
2712                        -0.19074074, -0.20925926, -0.19074074,
2713                        -0.20925926, -0.19074074, -0.20925926,
2714                        -0.19074074, -0.20925926, -0.19074074,
2715                        -0.26481481, -0.2462963, -0.26481481,
2716                        -0.2462963, -0.26481481, -0.2462963,
2717                        -0.26481481, -0.2462963, -0.26481481,
2718                        -0.2462963, -0.26481481, -0.2462963])
2719
2720
2721        #print domain.quantities['stage'].extrapolate_second_order()
2722        #domain.distribute_to_vertices_and_edges()
2723        #print domain.quantities['stage'].vertex_values[:,0]
2724
2725        #Evolution
2726        for t in domain.evolve(yieldstep = 0.05, finaltime = 0.5):
2727            pass
2728
2729
2730        assert allclose(domain.quantities['stage'].centroid_values,
2731     [-0.02907028, -0.01475478, -0.02973417, -0.01447186, -0.02932665, -0.01428285,
2732      -0.02901975, -0.0141361,  -0.02898816, -0.01418135, -0.02961409, -0.01403487,
2733      -0.07597998, -0.06252591, -0.07664854, -0.06312532, -0.07638287, -0.06265139,
2734      -0.07571145, -0.06235231, -0.0756817,  -0.06245309, -0.07652292, -0.06289946,
2735      -0.12367464, -0.11088981, -0.12237277, -0.11115338, -0.1218934,  -0.1107174,
2736      -0.12081485, -0.11000491, -0.12038451, -0.11010335, -0.12102113, -0.11012105,
2737      -0.16909116, -0.15831543, -0.16730214, -0.15786249, -0.1665493,  -0.15697919,
2738      -0.16496618, -0.15559852, -0.16338679, -0.15509088, -0.16364092, -0.15424423,
2739      -0.18771107, -0.19903904, -0.18903759, -0.19858437, -0.18701552, -0.19697797,
2740      -0.1833593,  -0.19505871, -0.1818806,  -0.19418042, -0.18586159, -0.19576946,
2741      -0.13986873, -0.14170053, -0.14132188, -0.14560674, -0.14095617, -0.14373292,
2742      -0.13785933, -0.14033364, -0.13592955, -0.13936356, -0.13596008, -0.14216296])
2743
2744        assert allclose(domain.quantities['xmomentum'].centroid_values,
2745     [ 0.00831121,  0.00317948,  0.00731797,  0.00334939,  0.00764717,  0.00348053,
2746       0.00788729,  0.00356522,  0.00780649,  0.00341919,  0.00693525,  0.00310375,
2747       0.02166196,  0.01421475,  0.02017737,  0.01316839,  0.02037015,  0.01368659,
2748       0.02106,     0.01399161,  0.02074514,  0.01354935,  0.01887407,  0.0123113,
2749       0.03775083,  0.02855197,  0.03689337,  0.02759782,  0.03732848,  0.02812072,
2750       0.03872545,  0.02913348,  0.03880939,  0.02803804,  0.03546499,  0.0260039,
2751       0.0632131,   0.04730634,  0.0576324,   0.04592336,  0.05790921,  0.04690514,
2752       0.05986467,  0.04871165,  0.06170068,  0.04811572,  0.05657041,  0.04416292,
2753       0.08489642,  0.07188097,  0.07835261,  0.06843406,  0.07986412,  0.0698247,
2754       0.08201071,  0.07216756,  0.08378418,  0.07273624,  0.080399,    0.06645841,
2755       0.01631548,  0.04691608,  0.0206632,   0.044409,    0.02115518,  0.04560305,
2756       0.02160608,  0.04663725,  0.02174734,  0.04795559,  0.02281427,  0.05667111])
2757
2758
2759        assert allclose(domain.quantities['ymomentum'].centroid_values,
2760     [ 1.45876601e-004, -3.24627393e-004, -1.57572719e-004, -2.92790187e-004,
2761      -9.90988382e-005, -3.06677335e-004, -1.62493106e-004, -3.71310004e-004,
2762      -1.99445058e-004, -3.28493467e-004,  6.68217349e-005, -8.42042805e-006,
2763       5.05093371e-004, -1.42842214e-004, -6.81454718e-005, -5.02084057e-004,
2764      -8.50583861e-005, -4.65443981e-004, -1.96406564e-004, -5.88889562e-004,
2765      -2.70160173e-004, -5.35485454e-004,  2.60780997e-004,  3.12145471e-005,
2766       5.16189608e-004,  1.07069062e-004,  9.29989252e-005, -3.71211119e-004,
2767       1.16350246e-004, -3.82407830e-004, -1.62077969e-004, -6.30906636e-004,
2768      -4.74025708e-004, -6.94463009e-004,  6.15092843e-005,  2.22106820e-004,
2769      -6.29589294e-004,  2.43611937e-004, -5.88125094e-004, -6.94293192e-005,
2770      -4.17914641e-004,  6.64609019e-005, -7.68334577e-004, -3.40232101e-004,
2771      -1.67424308e-003, -7.39485066e-004, -1.59966988e-003,  5.68262838e-005,
2772      -1.48470633e-003, -1.84554882e-003, -2.27200099e-003, -1.67506848e-003,
2773      -1.95610258e-003, -1.47638801e-003, -1.73779477e-003, -1.85498791e-003,
2774      -2.01357843e-003, -2.17675471e-003, -1.65783870e-003, -1.15818681e-003,
2775      -1.18663036e-003, -2.94229849e-003, -3.59309018e-003, -5.13496584e-003,
2776      -6.17359400e-003, -5.98761937e-003, -6.00540116e-003, -5.01121966e-003,
2777      -4.50964850e-003, -3.06319963e-003,  6.08950810e-004, -4.79537921e-004])
2778
2779
2780
2781
2782    def test_temp_play(self):
2783
2784        from mesh_factory import rectangular
2785        from shallow_water import Domain, Reflective_boundary, Constant_height
2786        from Numeric import array
2787
2788        #Create basic mesh
2789        points, vertices, boundary = rectangular(5, 5)
2790
2791        #Create shallow water domain
2792        domain = Domain(points, vertices, boundary)
2793        domain.smooth = False
2794        domain.default_order=2
2795        domain.beta_h = 0.0 #Use first order in h-limiter
2796
2797        #Bed-slope and friction at vertices (and interpolated elsewhere)
2798        def x_slope(x, y):
2799            return -x/3
2800
2801        domain.set_quantity('elevation', x_slope)
2802
2803        # Boundary conditions
2804        Br = Reflective_boundary(domain)
2805        domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
2806
2807        #Initial condition
2808        domain.set_quantity('stage', Constant_height(x_slope, 0.05))
2809        domain.check_integrity()
2810
2811        #Evolution
2812        for t in domain.evolve(yieldstep = 0.05, finaltime = 0.1):
2813            pass
2814
2815        assert allclose(domain.quantities['stage'].centroid_values[:4],
2816                        [0.00206836, 0.01296714, 0.00363415, 0.01438924])
2817        assert allclose(domain.quantities['xmomentum'].centroid_values[:4],
2818                        [0.01360154, 0.00671133, 0.01264578, 0.00648503])
2819        assert allclose(domain.quantities['ymomentum'].centroid_values[:4],
2820                        [-1.19201077e-003, -7.23647546e-004,
2821                        -6.39083123e-005, 6.29815168e-005])
2822
2823
2824    def test_complex_bed(self):
2825        #No friction is tested here
2826
2827        from shallow_water import Domain, Reflective_boundary, Dirichlet_boundary,\
2828             Transmissive_boundary, Time_boundary,\
2829             Weir_simple as Weir, Constant_height
2830
2831        from mesh_factory import rectangular
2832        from Numeric import array
2833
2834        N = 12
2835        points, vertices, boundary = rectangular(N, N/2, len1=1.2,len2=0.6,
2836                                                 origin=(-0.07, 0))
2837
2838
2839        domain = Domain(points, vertices, boundary)
2840        domain.smooth = False
2841        domain.visualise = False
2842        domain.default_order=2
2843
2844
2845        inflow_stage = 0.1
2846        Z = Weir(inflow_stage)
2847        domain.set_quantity('elevation', Z)
2848
2849        Br = Reflective_boundary(domain)
2850        Bd = Dirichlet_boundary([inflow_stage, 0.0, 0.0])
2851        domain.set_boundary({'left': Bd, 'right': Br, 'bottom': Br, 'top': Br})
2852
2853        domain.set_quantity('stage', Constant_height(Z, 0.))
2854
2855        for t in domain.evolve(yieldstep = 0.02, finaltime = 0.2):
2856            pass
2857
2858
2859        #print domain.quantities['stage'].centroid_values
2860
2861        #FIXME: These numbers were from version before 25/10
2862        #assert allclose(domain.quantities['stage'].centroid_values,
2863# [3.95822638e-002,  5.61022588e-002,  4.66437868e-002,  5.73081011e-002,
2864#  4.72394613e-002,  5.74684939e-002,  4.74309483e-002,  5.77458084e-002,
2865#  4.80628177e-002,  5.85656225e-002,  4.90498542e-002,  6.02609831e-002,
2866#  1.18470315e-002,  1.75136443e-002,  1.18035266e-002,  2.15565695e-002,
2867#  1.31620268e-002,  2.14351640e-002,  1.32351076e-002,  2.15450687e-002,
2868#  1.36414028e-002,  2.24274619e-002,  1.51689511e-002,  2.21789655e-002,
2869# -7.54337535e-003, -6.86362021e-004, -7.74146760e-003, -1.83756530e-003,
2870# -8.16773628e-003, -4.49916813e-004, -8.08202599e-003, -3.91118720e-004,
2871# -8.10292716e-003, -3.88584984e-004, -7.35226124e-003,  2.73985295e-004,
2872#  1.86166683e-001,  8.74070369e-002,  1.86166712e-001,  8.74035875e-002,
2873#  6.11666935e-002, -3.76173225e-002, -6.38333276e-002, -3.76147365e-002,
2874#  6.11666725e-002,  8.73846774e-002,  1.86166697e-001,  8.74171550e-002,
2875# -4.83333333e-002,  1.18333333e-001, -4.83333333e-002,  1.18333333e-001,
2876# -4.83333333e-002, -6.66666667e-003, -1.73333333e-001, -1.31666667e-001,
2877# -1.73333333e-001, -6.66666667e-003, -4.83333333e-002,  1.18333333e-001,
2878# -2.48333333e-001, -2.31666667e-001, -2.48333333e-001, -2.31666667e-001,
2879# -2.48333333e-001, -2.31666667e-001, -2.48333333e-001, -2.31666667e-001,
2880# -2.48333333e-001, -2.31666667e-001, -2.48333333e-001, -2.31666667e-001,
2881# -4.65000000e-001, -3.65000000e-001, -4.65000000e-001, -3.65000000e-001,
2882# -4.65000000e-001, -3.65000000e-001, -4.65000000e-001, -3.65000000e-001,
2883# -4.65000000e-001, -3.65000000e-001, -4.65000000e-001, -3.65000000e-001,
2884# -5.98333333e-001, -5.81666667e-001, -5.98333333e-001, -5.81666667e-001,
2885# -5.98333333e-001, -5.81666667e-001, -5.98333333e-001, -5.81666667e-001,
2886# -5.98333333e-001, -5.81666667e-001, -5.98333333e-001, -5.81666667e-001,
2887# -6.48333333e-001, -6.31666667e-001, -6.48333333e-001, -6.31666667e-001,
2888# -6.48333333e-001, -6.31666667e-001, -6.48333333e-001, -6.31666667e-001,
2889# -6.48333333e-001, -6.31666667e-001, -6.48333333e-001, -6.31666667e-001,
2890# -5.31666667e-001, -5.98333333e-001, -5.31666667e-001, -5.98333333e-001,
2891# -5.31666667e-001, -5.98333333e-001, -5.31666667e-001, -5.98333333e-001,
2892# -5.31666667e-001, -5.98333333e-001, -5.31666667e-001, -5.98333333e-001,
2893# -4.98333333e-001, -4.81666667e-001, -4.98333333e-001, -4.81666667e-001,
2894# -4.98333333e-001, -4.81666667e-001, -4.98333333e-001, -4.81666667e-001,
2895# -4.98333333e-001, -4.81666667e-001, -4.98333333e-001, -4.81666667e-001,
2896# -5.48333333e-001, -5.31666667e-001, -5.48333333e-001, -5.31666667e-001,
2897# -5.48333333e-001, -5.31666667e-001, -5.48333333e-001, -5.31666667e-001,
2898# -5.48333333e-001, -5.31666667e-001, -5.48333333e-001, -5.31666667e-001])
2899
2900
2901
2902    def test_spatio_temporal_boundary(self):
2903        """Test that boundary values can be read from file and interpolated
2904        in both time and space.
2905
2906        Verify that the same steady state solution is arrived at and that
2907        time interpolation works.
2908
2909        The full solution history is not exactly the same as
2910        file boundary must read and interpolate from *smoothed* version
2911        as stored in sww.
2912        """
2913        import time
2914
2915        #Create sww file of simple propagation from left to right
2916        #through rectangular domain
2917
2918        from mesh_factory import rectangular
2919
2920        #Create basic mesh
2921        points, vertices, boundary = rectangular(3, 3)
2922
2923        #Create shallow water domain
2924        domain1 = Domain(points, vertices, boundary)
2925
2926        from util import mean
2927        domain1.reduction = mean
2928        domain1.smooth = False  #Exact result
2929
2930        domain1.default_order = 2
2931        domain1.store = True
2932        domain1.set_datadir('.')
2933        domain1.set_name('spatio_temporal_boundary_source' + str(time.time()))
2934
2935        #FIXME: This is extremely important!
2936        #How can we test if they weren't stored?
2937        domain1.quantities_to_be_stored = ['stage', 'xmomentum', 'ymomentum']
2938
2939
2940        #Bed-slope and friction at vertices (and interpolated elsewhere)
2941        domain1.set_quantity('elevation', 0)
2942        domain1.set_quantity('friction', 0)
2943
2944        # Boundary conditions
2945        Br = Reflective_boundary(domain1)
2946        Bd = Dirichlet_boundary([0.3,0,0])
2947        domain1.set_boundary({'left': Bd, 'top': Bd, 'right': Br, 'bottom': Br})
2948        #Initial condition
2949        domain1.set_quantity('stage', 0)
2950        domain1.check_integrity()
2951
2952        finaltime = 5
2953        #Evolution  (full domain - large steps)
2954        for t in domain1.evolve(yieldstep = 0.671, finaltime = finaltime):
2955            pass
2956            #domain1.write_time()
2957
2958        cv1 = domain1.quantities['stage'].centroid_values
2959
2960
2961        #Create an triangle shaped domain (reusing coordinates from domain 1),
2962        #formed from the lower and right hand  boundaries and
2963        #the sw-ne diagonal
2964        #from domain 1. Call it domain2
2965
2966        points = [ [0,0], [1.0/3,0], [1.0/3,1.0/3],
2967                   [2.0/3,0], [2.0/3,1.0/3], [2.0/3,2.0/3],
2968                   [1,0], [1,1.0/3], [1,2.0/3], [1,1]]
2969
2970        vertices = [ [1,2,0], [3,4,1], [2,1,4], [4,5,2],
2971                     [6,7,3], [4,3,7], [7,8,4], [5,4,8], [8,9,5]]
2972
2973        boundary = { (0,1):'bottom', (1,1):'bottom', (4,1): 'bottom',
2974                     (4,2):'right', (6,2):'right', (8,2):'right',
2975                     (0,0):'diagonal', (3,0):'diagonal', (8,0):'diagonal'}
2976
2977        domain2 = Domain(points, vertices, boundary)
2978
2979        domain2.reduction = domain1.reduction
2980        domain2.smooth = False
2981        domain2.default_order = 2
2982
2983        #Bed-slope and friction at vertices (and interpolated elsewhere)
2984        domain2.set_quantity('elevation', 0)
2985        domain2.set_quantity('friction', 0)
2986        domain2.set_quantity('stage', 0)
2987
2988        # Boundary conditions
2989        Br = Reflective_boundary(domain2)
2990        Bf = Spatio_temporal_boundary(domain1.filename + '.' + domain1.format,
2991                                      domain2)
2992        domain2.set_boundary({'right':Br, 'bottom':Br, 'diagonal':Bf})
2993        domain2.check_integrity()
2994
2995
2996
2997        #Evolution (small steps)
2998        for t in domain2.evolve(yieldstep = 0.0711, finaltime = finaltime):
2999            pass
3000
3001
3002        #Use output from domain1 as spatio-temporal boundary for domain2
3003        #and verify that results at right hand side are close.
3004
3005        cv2 = domain2.quantities['stage'].centroid_values
3006
3007        #print take(cv1, (12,14,16))  #Right
3008        #print take(cv2, (4,6,8))
3009        #print take(cv1, (0,6,12))  #Bottom
3010        #print take(cv2, (0,1,4))
3011        #print take(cv1, (0,8,16))  #Diag
3012        #print take(cv2, (0,3,8))
3013
3014        assert allclose( take(cv1, (0,8,16)), take(cv2, (0,3,8))) #Diag
3015        assert allclose( take(cv1, (0,6,12)), take(cv2, (0,1,4))) #Bottom
3016        assert allclose( take(cv1, (12,14,16)), take(cv2, (4,6,8))) #RHS
3017
3018        #Cleanup
3019        os.remove(domain1.filename + '.' + domain1.format)
3020
3021
3022
3023    def test_spatio_temporal_boundary_2(self):
3024        """Test that boundary values can be read from file and interpolated
3025        in both time and space.
3026        This is a more basic test, verifying that boundary object
3027        produces the expected results
3028
3029
3030        """
3031        import time
3032
3033        #Create sww file of simple propagation from left to right
3034        #through rectangular domain
3035
3036        from mesh_factory import rectangular
3037
3038        #Create basic mesh
3039        points, vertices, boundary = rectangular(3, 3)
3040
3041        #Create shallow water domain
3042        domain1 = Domain(points, vertices, boundary)
3043
3044        from util import mean
3045        domain1.reduction = mean
3046        domain1.smooth = True #To mimic MOST output
3047
3048        domain1.default_order = 2
3049        domain1.store = True
3050        domain1.set_datadir('.')
3051        domain1.set_name('spatio_temporal_boundary_source' + str(time.time()))
3052
3053        #FIXME: This is extremely important!
3054        #How can we test if they weren't stored?
3055        domain1.quantities_to_be_stored = ['stage', 'xmomentum', 'ymomentum']
3056
3057
3058        #Bed-slope and friction at vertices (and interpolated elsewhere)
3059        domain1.set_quantity('elevation', 0)
3060        domain1.set_quantity('friction', 0)
3061
3062        # Boundary conditions
3063        Br = Reflective_boundary(domain1)
3064        Bd = Dirichlet_boundary([0.3,0,0])
3065        domain1.set_boundary({'left': Bd, 'top': Bd, 'right': Br, 'bottom': Br})
3066        #Initial condition
3067        domain1.set_quantity('stage', 0)
3068        domain1.check_integrity()
3069
3070        finaltime = 5
3071        #Evolution  (full domain - large steps)
3072        for t in domain1.evolve(yieldstep = 1, finaltime = finaltime):
3073            pass
3074            #domain1.write_time()
3075
3076
3077        #Create an triangle shaped domain (coinciding with some
3078        #coordinates from domain 1),
3079        #formed from the lower and right hand  boundaries and
3080        #the sw-ne diagonal
3081        #from domain 1. Call it domain2
3082
3083        points = [ [0,0], [1.0/3,0], [1.0/3,1.0/3],
3084                   [2.0/3,0], [2.0/3,1.0/3], [2.0/3,2.0/3],
3085                   [1,0], [1,1.0/3], [1,2.0/3], [1,1]]
3086
3087        vertices = [ [1,2,0],
3088                     [3,4,1], [2,1,4], [4,5,2],
3089                     [6,7,3], [4,3,7], [7,8,4], [5,4,8], [8,9,5]]
3090
3091        boundary = { (0,1):'bottom', (1,1):'bottom', (4,1): 'bottom',
3092                     (4,2):'right', (6,2):'right', (8,2):'right',
3093                     (0,0):'diagonal', (3,0):'diagonal', (8,0):'diagonal'}
3094
3095        domain2 = Domain(points, vertices, boundary)
3096
3097        domain2.reduction = domain1.reduction
3098        domain2.smooth = False
3099        domain2.default_order = 2
3100
3101        #Bed-slope and friction at vertices (and interpolated elsewhere)
3102        domain2.set_quantity('elevation', 0)
3103        domain2.set_quantity('friction', 0)
3104        domain2.set_quantity('stage', 0)
3105
3106
3107        #Read results for specific timesteps t=1 and t=2
3108        from Scientific.IO.NetCDF import NetCDFFile
3109        fid = NetCDFFile(domain1.filename + '.' + domain1.format)
3110
3111        x = fid.variables['x'][:]
3112        y = fid.variables['y'][:]
3113        s1 = fid.variables['stage'][1,:]
3114        s2 = fid.variables['stage'][2,:]
3115        fid.close()
3116
3117        from Numeric import take, reshape, concatenate
3118        shp = (len(x), 1)
3119        points = concatenate( (reshape(x, shp), reshape(y, shp)), axis=1)
3120        #The diagonal points of domain 1 are 0, 5, 10, 15
3121
3122        #print points[0], points[5], points[10], points[15]
3123        assert allclose( take(points, [0,5,10,15]),
3124                         [[0,0], [1.0/3, 1.0/3], [2.0/3, 2.0/3], [1,1]])
3125
3126
3127        # Boundary conditions
3128        Br = Reflective_boundary(domain2)
3129        Bf = Spatio_temporal_boundary(domain1.filename + '.' + domain1.format,
3130                                      domain2)
3131        domain2.set_boundary({'right':Br, 'bottom':Br, 'diagonal':Bf})
3132        domain2.check_integrity()
3133
3134        #Test that interpolation points are the mid points of the all boundary
3135        #segments
3136
3137        boundary_midpoints = [[1.0/6, 0], [1.0/2, 0], [5.0/6,0],
3138                              [1.0, 1.0/6], [1.0, 1.0/2], [1.0, 5.0/6],
3139                              [1.0/6, 1.0/6], [0.5, 0.5], [5.0/6, 5.0/6]]
3140
3141        boundary_midpoints.sort()
3142        R = Bf.F.interpolation_points.tolist()
3143        R.sort()
3144        assert allclose(boundary_midpoints, R)
3145
3146        #Check spatially interpolated output at time == 1
3147        domain2.time = 1
3148
3149        #First diagonal midpoint
3150        R0 = Bf.evaluate(0,0)
3151        assert allclose(R0[0], (s1[0] + s1[5])/2)
3152
3153        #Second diagonal midpoint
3154        R0 = Bf.evaluate(3,0)
3155        assert allclose(R0[0], (s1[5] + s1[10])/2)
3156
3157        #First diagonal midpoint
3158        R0 = Bf.evaluate(8,0)
3159        assert allclose(R0[0], (s1[10] + s1[15])/2)
3160
3161        #Check spatially interpolated output at time == 2
3162        domain2.time = 2
3163
3164        #First diagonal midpoint
3165        R0 = Bf.evaluate(0,0)
3166        assert allclose(R0[0], (s2[0] + s2[5])/2)
3167
3168        #Second diagonal midpoint
3169        R0 = Bf.evaluate(3,0)
3170        assert allclose(R0[0], (s2[5] + s2[10])/2)
3171
3172        #First diagonal midpoint
3173        R0 = Bf.evaluate(8,0)
3174        assert allclose(R0[0], (s2[10] + s2[15])/2)
3175
3176
3177        #Now check temporal interpolation
3178
3179        domain2.time = 1 + 2.0/3
3180
3181        #First diagonal midpoint
3182        R0 = Bf.evaluate(0,0)
3183        assert allclose(R0[0], ((s1[0] + s1[5])/2 + 2.0*(s2[0] + s2[5])/2)/3)
3184
3185        #Second diagonal midpoint
3186        R0 = Bf.evaluate(3,0)
3187        assert allclose(R0[0], ((s1[5] + s1[10])/2 + 2.0*(s2[5] + s2[10])/2)/3)
3188
3189        #First diagonal midpoint
3190        R0 = Bf.evaluate(8,0)
3191        assert allclose(R0[0], ((s1[10] + s1[15])/2 + 2.0*(s2[10] + s2[15])/2)/3)
3192
3193
3194
3195        #Cleanup
3196        os.remove(domain1.filename + '.' + domain1.format)
3197
3198
3199        #-------------------------------------------------------------
3200if __name__ == "__main__":
3201    suite = unittest.makeSuite(Test_Shallow_Water,'test')
3202    runner = unittest.TextTestRunner()
3203    runner.run(suite)
Note: See TracBrowser for help on using the repository browser.