source: inundation/pyvolution/test_shallow_water.py @ 2366

Last change on this file since 2366 was 2305, checked in by ole, 19 years ago

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