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

Last change on this file since 5442 was 5442, checked in by ole, 15 years ago

Retired h-limiter and beta_h as per ticket:194.
All unit tests and validation tests pass.

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