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

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

Prepared tests for making new slope limiters default. Still a few unprepared tests remaining and
one example in test_data_manager where small timesteps are generated.

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