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

Last change on this file since 1027 was 1027, checked in by steve, 19 years ago

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