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

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

Ensured all beta_h == 0.0, test and validations passed.

File size: 209.5 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.beta_h = 0.0 # Backwards compatibility (14/2/7)
1439        domain.use_centroid_velocities = 0 # Backwards compatibility (7/5/8)
1440       
1441
1442        #-----------------------------------------------------------------
1443        # Setup initial conditions
1444        #-----------------------------------------------------------------
1445
1446        def topography(x,y): 
1447            return -x/2                              # linear bed slope
1448
1449        domain.set_quantity('elevation', topography) 
1450        domain.set_quantity('friction', 0.0)         
1451        domain.set_quantity('stage', expression='elevation')           
1452
1453
1454        #----------------------------------------------------------------
1455        # Setup boundary conditions
1456        #----------------------------------------------------------------
1457
1458        Br = Reflective_boundary(domain)      # Solid reflective wall
1459        Bd = Dirichlet_boundary([-0.2,0.,0.]) # Constant boundary values
1460        domain.set_boundary({'left': Br, 'right': Bd, 'top': Br, 'bottom': Br})
1461
1462
1463        #----------------------------------------------------------------
1464        # Evolve system through time
1465        #----------------------------------------------------------------
1466
1467        interpolation_points = [[0.4,0.5], [0.6,0.5], [0.8,0.5], [0.9,0.5]]
1468        gauge_values = []
1469        for _ in interpolation_points:
1470            gauge_values.append([])
1471
1472        time = []
1473        for t in domain.evolve(yieldstep = 0.1, finaltime = 5.0):
1474            # Record time series at known points
1475            time.append(domain.get_time())
1476           
1477            stage = domain.get_quantity('stage')
1478            w = stage.get_values(interpolation_points=interpolation_points)
1479           
1480            for i, _ in enumerate(interpolation_points):
1481                gauge_values[i].append(w[i])
1482
1483
1484        #print
1485        #print time
1486        #print
1487        #for i, (x,y) in enumerate(interpolation_points):
1488        #    print i, gauge_values[i]
1489        #    print
1490
1491        #Reference (nautilus 26/6/2008) (beta_h == 0)
1492       
1493        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]
1494
1495        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]
1496       
1497        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]
1498       
1499        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]
1500       
1501        #FIXME (DSG):This is a hack so the anuga install, not precompiled
1502        # works on DSG's win2000, python 2.3
1503        #The problem is the gauge_values[X] are 52 long, not 51.
1504        #
1505        # This was probably fixed by Stephen in changeset:3804
1506        #if len(gauge_values[0]) == 52: gauge_values[0].pop()
1507        #if len(gauge_values[1]) == 52: gauge_values[1].pop()
1508        #if len(gauge_values[2]) == 52: gauge_values[2].pop()
1509        #if len(gauge_values[3]) == 52: gauge_values[3].pop()
1510
1511##         print len(G0), len(gauge_values[0])
1512##         print len(G1), len(gauge_values[1])
1513       
1514        #print gauge_values[3]
1515        #print G0[:4]
1516        #print array(gauge_values[0])-array(G0)
1517       
1518       
1519        assert allclose(gauge_values[0], G0)
1520        assert allclose(gauge_values[1], G1)
1521        assert allclose(gauge_values[2], G2)
1522        assert allclose(gauge_values[3], G3)       
1523
1524
1525
1526
1527
1528
1529
1530    #####################################################
1531
1532    def test_flux_optimisation(self):
1533        """test_flux_optimisation
1534        Test that fluxes are correctly computed using
1535        dry and still cell exclusions
1536        """
1537
1538        from anuga.config import g
1539        import copy
1540
1541        a = [0.0, 0.0]
1542        b = [0.0, 2.0]
1543        c = [2.0, 0.0]
1544        d = [0.0, 4.0]
1545        e = [2.0, 2.0]
1546        f = [4.0, 0.0]
1547
1548        points = [a, b, c, d, e, f]
1549        #bac, bce, ecf, dbe
1550        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
1551
1552        domain = Domain(points, vertices)
1553
1554        #Set up for a gradient of (3,0) at mid triangle (bce)
1555        def slope(x, y):
1556            return 3*x
1557
1558        h = 0.1
1559        def stage(x,y):
1560            return slope(x,y)+h
1561
1562        domain.set_quantity('elevation', slope)
1563        domain.set_quantity('stage', stage)
1564
1565        # Allow slope limiters to work (FIXME (Ole): Shouldn't this be automatic in ANUGA?)     
1566        domain.distribute_to_vertices_and_edges()       
1567
1568        initial_stage = copy.copy(domain.quantities['stage'].vertex_values)
1569
1570        domain.set_boundary({'exterior': Reflective_boundary(domain)})
1571
1572
1573        #  Check that update arrays are initialised to zero
1574        assert allclose(domain.get_quantity('stage').explicit_update, 0)
1575        assert allclose(domain.get_quantity('xmomentum').explicit_update, 0)
1576        assert allclose(domain.get_quantity('ymomentum').explicit_update, 0)               
1577
1578
1579        # Get true values
1580        domain.optimise_dry_cells = False
1581        domain.compute_fluxes()
1582        stage_ref = copy.copy(domain.get_quantity('stage').explicit_update)
1583        xmom_ref = copy.copy(domain.get_quantity('xmomentum').explicit_update)
1584        ymom_ref = copy.copy(domain.get_quantity('ymomentum').explicit_update)       
1585
1586        # Try with flux optimisation
1587        domain.optimise_dry_cells = True
1588        domain.compute_fluxes()
1589
1590        assert allclose(stage_ref, domain.get_quantity('stage').explicit_update)
1591        assert allclose(xmom_ref, domain.get_quantity('xmomentum').explicit_update)
1592        assert allclose(ymom_ref, domain.get_quantity('ymomentum').explicit_update)
1593       
1594   
1595       
1596    def test_initial_condition(self):
1597        """test_initial_condition
1598        Test that initial condition is output at time == 0 and that
1599        computed values change as system evolves
1600        """
1601
1602        from anuga.config import g
1603        import copy
1604
1605        a = [0.0, 0.0]
1606        b = [0.0, 2.0]
1607        c = [2.0, 0.0]
1608        d = [0.0, 4.0]
1609        e = [2.0, 2.0]
1610        f = [4.0, 0.0]
1611
1612        points = [a, b, c, d, e, f]
1613        #bac, bce, ecf, dbe
1614        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
1615
1616        domain = Domain(points, vertices)
1617
1618        #Set up for a gradient of (3,0) at mid triangle (bce)
1619        def slope(x, y):
1620            return 3*x
1621
1622        h = 0.1
1623        def stage(x,y):
1624            return slope(x,y)+h
1625
1626        domain.set_quantity('elevation', slope)
1627        domain.set_quantity('stage', stage)
1628
1629        # Allow slope limiters to work (FIXME (Ole): Shouldn't this be automatic in ANUGA?)     
1630        domain.distribute_to_vertices_and_edges()       
1631
1632        initial_stage = copy.copy(domain.quantities['stage'].vertex_values)
1633
1634        domain.set_boundary({'exterior': Reflective_boundary(domain)})
1635
1636        domain.optimise_dry_cells = True
1637        #Evolution
1638        for t in domain.evolve(yieldstep = 0.5, finaltime = 2.0):
1639            stage = domain.quantities['stage'].vertex_values
1640
1641            if t == 0.0:
1642                assert allclose(stage, initial_stage)
1643            else:
1644                assert not allclose(stage, initial_stage)
1645
1646
1647        os.remove(domain.get_name() + '.sww')
1648
1649
1650
1651    #####################################################
1652    def test_gravity(self):
1653        #Assuming no friction
1654
1655        from anuga.config import g
1656
1657        a = [0.0, 0.0]
1658        b = [0.0, 2.0]
1659        c = [2.0, 0.0]
1660        d = [0.0, 4.0]
1661        e = [2.0, 2.0]
1662        f = [4.0, 0.0]
1663
1664        points = [a, b, c, d, e, f]
1665        #bac, bce, ecf, dbe
1666        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
1667
1668        domain = Domain(points, vertices)
1669
1670        #Set up for a gradient of (3,0) at mid triangle (bce)
1671        def slope(x, y):
1672            return 3*x
1673
1674        h = 0.1
1675        def stage(x,y):
1676            return slope(x,y)+h
1677
1678        domain.set_quantity('elevation', slope)
1679        domain.set_quantity('stage', stage)
1680
1681        for name in domain.conserved_quantities:
1682            assert allclose(domain.quantities[name].explicit_update, 0)
1683            assert allclose(domain.quantities[name].semi_implicit_update, 0)
1684
1685        domain.compute_forcing_terms()
1686
1687        assert allclose(domain.quantities['stage'].explicit_update, 0)
1688        assert allclose(domain.quantities['xmomentum'].explicit_update, -g*h*3)
1689        assert allclose(domain.quantities['ymomentum'].explicit_update, 0)
1690
1691
1692    def test_manning_friction(self):
1693        from anuga.config import g
1694
1695        a = [0.0, 0.0]
1696        b = [0.0, 2.0]
1697        c = [2.0, 0.0]
1698        d = [0.0, 4.0]
1699        e = [2.0, 2.0]
1700        f = [4.0, 0.0]
1701
1702        points = [a, b, c, d, e, f]
1703        #bac, bce, ecf, dbe
1704        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
1705
1706        domain = Domain(points, vertices)
1707
1708        #Set up for a gradient of (3,0) at mid triangle (bce)
1709        def slope(x, y):
1710            return 3*x
1711
1712        h = 0.1
1713        def stage(x,y):
1714            return slope(x,y)+h
1715
1716        eta = 0.07
1717        domain.set_quantity('elevation', slope)
1718        domain.set_quantity('stage', stage)
1719        domain.set_quantity('friction', eta)
1720
1721        for name in domain.conserved_quantities:
1722            assert allclose(domain.quantities[name].explicit_update, 0)
1723            assert allclose(domain.quantities[name].semi_implicit_update, 0)
1724
1725        domain.compute_forcing_terms()
1726
1727        assert allclose(domain.quantities['stage'].explicit_update, 0)
1728        assert allclose(domain.quantities['xmomentum'].explicit_update, -g*h*3)
1729        assert allclose(domain.quantities['ymomentum'].explicit_update, 0)
1730
1731        assert allclose(domain.quantities['stage'].semi_implicit_update, 0)
1732        assert allclose(domain.quantities['xmomentum'].semi_implicit_update, 0)
1733        assert allclose(domain.quantities['ymomentum'].semi_implicit_update, 0)
1734
1735        #Create some momentum for friction to work with
1736        domain.set_quantity('xmomentum', 1)
1737        S = -g * eta**2 / h**(7.0/3)
1738
1739        domain.compute_forcing_terms()
1740        assert allclose(domain.quantities['stage'].semi_implicit_update, 0)
1741        assert allclose(domain.quantities['xmomentum'].semi_implicit_update, S)
1742        assert allclose(domain.quantities['ymomentum'].semi_implicit_update, 0)
1743
1744        #A more complex example
1745        domain.quantities['stage'].semi_implicit_update[:] = 0.0
1746        domain.quantities['xmomentum'].semi_implicit_update[:] = 0.0
1747        domain.quantities['ymomentum'].semi_implicit_update[:] = 0.0
1748
1749        domain.set_quantity('xmomentum', 3)
1750        domain.set_quantity('ymomentum', 4)
1751
1752        S = -g * eta**2 * 5 / h**(7.0/3)
1753
1754
1755        domain.compute_forcing_terms()
1756
1757        assert allclose(domain.quantities['stage'].semi_implicit_update, 0)
1758        assert allclose(domain.quantities['xmomentum'].semi_implicit_update, 3*S)
1759        assert allclose(domain.quantities['ymomentum'].semi_implicit_update, 4*S)
1760
1761    def test_constant_wind_stress(self):
1762        from anuga.config import rho_a, rho_w, eta_w
1763        from math import pi, cos, sin
1764
1765        a = [0.0, 0.0]
1766        b = [0.0, 2.0]
1767        c = [2.0, 0.0]
1768        d = [0.0, 4.0]
1769        e = [2.0, 2.0]
1770        f = [4.0, 0.0]
1771
1772        points = [a, b, c, d, e, f]
1773        #bac, bce, ecf, dbe
1774        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
1775
1776
1777        domain = Domain(points, vertices)
1778
1779        #Flat surface with 1m of water
1780        domain.set_quantity('elevation', 0)
1781        domain.set_quantity('stage', 1.0)
1782        domain.set_quantity('friction', 0)
1783
1784        Br = Reflective_boundary(domain)
1785        domain.set_boundary({'exterior': Br})
1786
1787        #Setup only one forcing term, constant wind stress
1788        s = 100
1789        phi = 135
1790        domain.forcing_terms = []
1791        domain.forcing_terms.append( Wind_stress(s, phi) )
1792
1793        domain.compute_forcing_terms()
1794
1795
1796        const = eta_w*rho_a/rho_w
1797
1798        #Convert to radians
1799        phi = phi*pi/180
1800
1801        #Compute velocity vector (u, v)
1802        u = s*cos(phi)
1803        v = s*sin(phi)
1804
1805        #Compute wind stress
1806        S = const * sqrt(u**2 + v**2)
1807
1808        assert allclose(domain.quantities['stage'].explicit_update, 0)
1809        assert allclose(domain.quantities['xmomentum'].explicit_update, S*u)
1810        assert allclose(domain.quantities['ymomentum'].explicit_update, S*v)
1811
1812
1813    def test_variable_wind_stress(self):
1814        from anuga.config import rho_a, rho_w, eta_w
1815        from math import pi, cos, sin
1816
1817        a = [0.0, 0.0]
1818        b = [0.0, 2.0]
1819        c = [2.0, 0.0]
1820        d = [0.0, 4.0]
1821        e = [2.0, 2.0]
1822        f = [4.0, 0.0]
1823
1824        points = [a, b, c, d, e, f]
1825        #bac, bce, ecf, dbe
1826        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
1827
1828        domain = Domain(points, vertices)
1829
1830        #Flat surface with 1m of water
1831        domain.set_quantity('elevation', 0)
1832        domain.set_quantity('stage', 1.0)
1833        domain.set_quantity('friction', 0)
1834
1835        Br = Reflective_boundary(domain)
1836        domain.set_boundary({'exterior': Br})
1837
1838
1839        domain.time = 5.54 #Take a random time (not zero)
1840
1841        #Setup only one forcing term, constant wind stress
1842        s = 100
1843        phi = 135
1844        domain.forcing_terms = []
1845        domain.forcing_terms.append( Wind_stress(s = speed, phi = angle) )
1846
1847        domain.compute_forcing_terms()
1848
1849        #Compute reference solution
1850        const = eta_w*rho_a/rho_w
1851
1852        N = len(domain) # number_of_triangles
1853
1854        xc = domain.get_centroid_coordinates()
1855        t = domain.time
1856
1857        x = xc[:,0]
1858        y = xc[:,1]
1859        s_vec = speed(t,x,y)
1860        phi_vec = angle(t,x,y)
1861
1862
1863        for k in range(N):
1864            #Convert to radians
1865            phi = phi_vec[k]*pi/180
1866            s = s_vec[k]
1867
1868            #Compute velocity vector (u, v)
1869            u = s*cos(phi)
1870            v = s*sin(phi)
1871
1872            #Compute wind stress
1873            S = const * sqrt(u**2 + v**2)
1874
1875            assert allclose(domain.quantities['stage'].explicit_update[k], 0)
1876            assert allclose(domain.quantities['xmomentum'].explicit_update[k], S*u)
1877            assert allclose(domain.quantities['ymomentum'].explicit_update[k], S*v)
1878
1879
1880
1881
1882
1883
1884    def test_windfield_from_file(self):
1885        from anuga.config import rho_a, rho_w, eta_w
1886        from math import pi, cos, sin
1887        from anuga.config import time_format
1888        from anuga.abstract_2d_finite_volumes.util import file_function
1889        import time
1890
1891
1892        a = [0.0, 0.0]
1893        b = [0.0, 2.0]
1894        c = [2.0, 0.0]
1895        d = [0.0, 4.0]
1896        e = [2.0, 2.0]
1897        f = [4.0, 0.0]
1898
1899        points = [a, b, c, d, e, f]
1900        #bac, bce, ecf, dbe
1901        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
1902
1903        domain = Domain(points, vertices)
1904
1905        #Flat surface with 1m of water
1906        domain.set_quantity('elevation', 0)
1907        domain.set_quantity('stage', 1.0)
1908        domain.set_quantity('friction', 0)
1909
1910        Br = Reflective_boundary(domain)
1911        domain.set_boundary({'exterior': Br})
1912
1913
1914        domain.time = 7 #Take a time that is represented in file (not zero)
1915
1916        #Write wind stress file (ensure that domain.time is covered)
1917        #Take x=1 and y=0
1918        filename = 'test_windstress_from_file'
1919        start = time.mktime(time.strptime('2000', '%Y'))
1920        fid = open(filename + '.txt', 'w')
1921        dt = 1  #One second interval
1922        t = 0.0
1923        while t <= 10.0:
1924            t_string = time.strftime(time_format, time.gmtime(t+start))
1925
1926            fid.write('%s, %f %f\n' %(t_string,
1927                                      speed(t,[1],[0])[0],
1928                                      angle(t,[1],[0])[0]))
1929            t += dt
1930
1931        fid.close()
1932
1933
1934        #Convert ASCII file to NetCDF (Which is what we really like!)
1935        from data_manager import timefile2netcdf       
1936        timefile2netcdf(filename)
1937        os.remove(filename + '.txt')
1938
1939       
1940        #Setup wind stress
1941        F = file_function(filename + '.tms', quantities = ['Attribute0',
1942                                                           'Attribute1'])
1943        os.remove(filename + '.tms')
1944       
1945
1946        #print 'F(5)', F(5)
1947
1948        #print 'F(5,x,y)', F(5,x=zeros(3),y=zeros(3))
1949       
1950        #print dir(F)
1951        #print F.T
1952        #print F.precomputed_values
1953        #
1954        #F = file_function(filename + '.txt')       
1955        #
1956        #print dir(F)
1957        #print F.T       
1958        #print F.Q
1959       
1960        W = Wind_stress(F)
1961
1962        domain.forcing_terms = []
1963        domain.forcing_terms.append(W)
1964
1965        domain.compute_forcing_terms()
1966
1967        #Compute reference solution
1968        const = eta_w*rho_a/rho_w
1969
1970        N = len(domain) # number_of_triangles
1971
1972        t = domain.time
1973
1974        s = speed(t,[1],[0])[0]
1975        phi = angle(t,[1],[0])[0]
1976
1977        #Convert to radians
1978        phi = phi*pi/180
1979
1980
1981        #Compute velocity vector (u, v)
1982        u = s*cos(phi)
1983        v = s*sin(phi)
1984
1985        #Compute wind stress
1986        S = const * sqrt(u**2 + v**2)
1987
1988        for k in range(N):
1989            assert allclose(domain.quantities['stage'].explicit_update[k], 0)
1990            assert allclose(domain.quantities['xmomentum'].explicit_update[k], S*u)
1991            assert allclose(domain.quantities['ymomentum'].explicit_update[k], S*v)
1992
1993
1994    def test_windfield_from_file_seconds(self):
1995        from anuga.config import rho_a, rho_w, eta_w
1996        from math import pi, cos, sin
1997        from anuga.config import time_format
1998        from anuga.abstract_2d_finite_volumes.util import file_function
1999        import time
2000
2001
2002        a = [0.0, 0.0]
2003        b = [0.0, 2.0]
2004        c = [2.0, 0.0]
2005        d = [0.0, 4.0]
2006        e = [2.0, 2.0]
2007        f = [4.0, 0.0]
2008
2009        points = [a, b, c, d, e, f]
2010        #bac, bce, ecf, dbe
2011        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
2012
2013        domain = Domain(points, vertices)
2014
2015        #Flat surface with 1m of water
2016        domain.set_quantity('elevation', 0)
2017        domain.set_quantity('stage', 1.0)
2018        domain.set_quantity('friction', 0)
2019
2020        Br = Reflective_boundary(domain)
2021        domain.set_boundary({'exterior': Br})
2022
2023
2024        domain.time = 7 #Take a time that is represented in file (not zero)
2025
2026        #Write wind stress file (ensure that domain.time is covered)
2027        #Take x=1 and y=0
2028        filename = 'test_windstress_from_file'
2029        start = time.mktime(time.strptime('2000', '%Y'))
2030        fid = open(filename + '.txt', 'w')
2031        dt = 0.5 #1  #One second interval
2032        t = 0.0
2033        while t <= 10.0:
2034            fid.write('%s, %f %f\n' %(str(t),
2035                                      speed(t,[1],[0])[0],
2036                                      angle(t,[1],[0])[0]))
2037            t += dt
2038
2039        fid.close()
2040
2041
2042        #Convert ASCII file to NetCDF (Which is what we really like!)
2043        from data_manager import timefile2netcdf       
2044        timefile2netcdf(filename, time_as_seconds=True)
2045        os.remove(filename + '.txt')
2046
2047       
2048        #Setup wind stress
2049        F = file_function(filename + '.tms', quantities = ['Attribute0',
2050                                                           'Attribute1'])
2051        os.remove(filename + '.tms')
2052       
2053
2054        #print 'F(5)', F(5)
2055
2056        #print 'F(5,x,y)', F(5,x=zeros(3),y=zeros(3))
2057       
2058        #print dir(F)
2059        #print F.T
2060        #print F.precomputed_values
2061        #
2062        #F = file_function(filename + '.txt')       
2063        #
2064        #print dir(F)
2065        #print F.T       
2066        #print F.Q
2067       
2068        W = Wind_stress(F)
2069
2070        domain.forcing_terms = []
2071        domain.forcing_terms.append(W)
2072
2073        domain.compute_forcing_terms()
2074
2075        #Compute reference solution
2076        const = eta_w*rho_a/rho_w
2077
2078        N = len(domain) # number_of_triangles
2079
2080        t = domain.time
2081
2082        s = speed(t,[1],[0])[0]
2083        phi = angle(t,[1],[0])[0]
2084
2085        #Convert to radians
2086        phi = phi*pi/180
2087
2088
2089        #Compute velocity vector (u, v)
2090        u = s*cos(phi)
2091        v = s*sin(phi)
2092
2093        #Compute wind stress
2094        S = const * sqrt(u**2 + v**2)
2095
2096        for k in range(N):
2097            assert allclose(domain.quantities['stage'].explicit_update[k], 0)
2098            assert allclose(domain.quantities['xmomentum'].explicit_update[k], S*u)
2099            assert allclose(domain.quantities['ymomentum'].explicit_update[k], S*v)
2100
2101
2102       
2103
2104    def test_wind_stress_error_condition(self):
2105        """Test that windstress reacts properly when forcing functions
2106        are wrong - e.g. returns a scalar
2107        """
2108
2109        from anuga.config import rho_a, rho_w, eta_w
2110        from math import pi, cos, sin
2111
2112        a = [0.0, 0.0]
2113        b = [0.0, 2.0]
2114        c = [2.0, 0.0]
2115        d = [0.0, 4.0]
2116        e = [2.0, 2.0]
2117        f = [4.0, 0.0]
2118
2119        points = [a, b, c, d, e, f]
2120        #bac, bce, ecf, dbe
2121        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
2122
2123        domain = Domain(points, vertices)
2124
2125        #Flat surface with 1m of water
2126        domain.set_quantity('elevation', 0)
2127        domain.set_quantity('stage', 1.0)
2128        domain.set_quantity('friction', 0)
2129
2130        Br = Reflective_boundary(domain)
2131        domain.set_boundary({'exterior': Br})
2132
2133
2134        domain.time = 5.54 #Take a random time (not zero)
2135
2136        #Setup only one forcing term, bad func
2137        domain.forcing_terms = []
2138
2139        try:
2140            domain.forcing_terms.append(Wind_stress(s = scalar_func,
2141                                                    phi = angle))
2142        except AssertionError:
2143            pass
2144        else:
2145            msg = 'Should have raised exception'
2146            raise msg
2147
2148
2149        try:
2150            domain.forcing_terms.append(Wind_stress(s = speed,
2151                                                    phi = scalar_func))
2152        except AssertionError:
2153            pass
2154        else:
2155            msg = 'Should have raised exception'
2156            raise msg
2157
2158        try:
2159            domain.forcing_terms.append(Wind_stress(s = speed,
2160                                                    phi = 'xx'))
2161        except:
2162            pass
2163        else:
2164            msg = 'Should have raised exception'
2165            raise msg
2166
2167
2168
2169    def test_rainfall(self):
2170        from math import pi, cos, sin
2171
2172        a = [0.0, 0.0]
2173        b = [0.0, 2.0]
2174        c = [2.0, 0.0]
2175        d = [0.0, 4.0]
2176        e = [2.0, 2.0]
2177        f = [4.0, 0.0]
2178
2179        points = [a, b, c, d, e, f]
2180        #bac, bce, ecf, dbe
2181        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
2182
2183
2184        domain = Domain(points, vertices)
2185
2186        #Flat surface with 1m of water
2187        domain.set_quantity('elevation', 0)
2188        domain.set_quantity('stage', 1.0)
2189        domain.set_quantity('friction', 0)
2190
2191        Br = Reflective_boundary(domain)
2192        domain.set_boundary({'exterior': Br})
2193
2194        # Setup only one forcing term, constant rainfall
2195        domain.forcing_terms = []
2196        domain.forcing_terms.append( Rainfall(domain, rate=2.0) )
2197
2198        domain.compute_forcing_terms()
2199        assert allclose(domain.quantities['stage'].explicit_update, 2.0/1000)
2200
2201
2202        # FIXME: Do Time dependent rainfall
2203
2204
2205
2206    def test_rainfall_restricted_by_polygon(self):
2207        from math import pi, cos, sin
2208
2209        a = [0.0, 0.0]
2210        b = [0.0, 2.0]
2211        c = [2.0, 0.0]
2212        d = [0.0, 4.0]
2213        e = [2.0, 2.0]
2214        f = [4.0, 0.0]
2215
2216        points = [a, b, c, d, e, f]
2217        #bac, bce, ecf, dbe
2218        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
2219
2220
2221        domain = Domain(points, vertices)
2222
2223        #Flat surface with 1m of water
2224        domain.set_quantity('elevation', 0)
2225        domain.set_quantity('stage', 1.0)
2226        domain.set_quantity('friction', 0)
2227
2228        Br = Reflective_boundary(domain)
2229        domain.set_boundary({'exterior': Br})
2230
2231        # Setup only one forcing term, constant rainfall restricted to a polygon enclosing triangle #1 (bce)
2232        domain.forcing_terms = []
2233        R = Rainfall(domain, rate=2.0, polygon = [[1,1], [2,1], [2,2], [1,2]])
2234
2235        assert allclose(R.exchange_area, 1)
2236       
2237        domain.forcing_terms.append(R)
2238
2239        domain.compute_forcing_terms()
2240        #print domain.quantities['stage'].explicit_update
2241       
2242        assert allclose(domain.quantities['stage'].explicit_update[1], 2.0/1000)
2243        assert allclose(domain.quantities['stage'].explicit_update[0], 0)
2244        assert allclose(domain.quantities['stage'].explicit_update[2:], 0)       
2245       
2246        # FIXME: Do Time dependent rainfall with poly
2247
2248
2249
2250
2251    def test_inflow_using_circle(self):
2252        from math import pi, cos, sin
2253
2254        a = [0.0, 0.0]
2255        b = [0.0, 2.0]
2256        c = [2.0, 0.0]
2257        d = [0.0, 4.0]
2258        e = [2.0, 2.0]
2259        f = [4.0, 0.0]
2260
2261        points = [a, b, c, d, e, f]
2262        #bac, bce, ecf, dbe
2263        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
2264
2265
2266        domain = Domain(points, vertices)
2267
2268        # Flat surface with 1m of water
2269        domain.set_quantity('elevation', 0)
2270        domain.set_quantity('stage', 1.0)
2271        domain.set_quantity('friction', 0)
2272
2273        Br = Reflective_boundary(domain)
2274        domain.set_boundary({'exterior': Br})
2275
2276        # Setup only one forcing term, constant inflow of 2 m^3/s on a circle affecting triangles #0 and #1 (bac and bce)
2277        domain.forcing_terms = []
2278        domain.forcing_terms.append( Inflow(domain, rate=2.0, center=(1,1), radius=1) )
2279
2280        domain.compute_forcing_terms()
2281        #print domain.quantities['stage'].explicit_update
2282       
2283        assert allclose(domain.quantities['stage'].explicit_update[1], 2.0/pi)
2284        assert allclose(domain.quantities['stage'].explicit_update[0], 2.0/pi)
2285        assert allclose(domain.quantities['stage'].explicit_update[2:], 0)       
2286
2287
2288    def test_inflow_using_circle_function(self):
2289        from math import pi, cos, sin
2290
2291        a = [0.0, 0.0]
2292        b = [0.0, 2.0]
2293        c = [2.0, 0.0]
2294        d = [0.0, 4.0]
2295        e = [2.0, 2.0]
2296        f = [4.0, 0.0]
2297
2298        points = [a, b, c, d, e, f]
2299        #bac, bce, ecf, dbe
2300        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
2301
2302
2303        domain = Domain(points, vertices)
2304
2305        # Flat surface with 1m of water
2306        domain.set_quantity('elevation', 0)
2307        domain.set_quantity('stage', 1.0)
2308        domain.set_quantity('friction', 0)
2309
2310        Br = Reflective_boundary(domain)
2311        domain.set_boundary({'exterior': Br})
2312
2313        # Setup only one forcing term, time dependent inflow of 2 m^3/s on a circle affecting triangles #0 and #1 (bac and bce)
2314        domain.forcing_terms = []
2315        domain.forcing_terms.append( Inflow(domain, rate=lambda t: 2., center=(1,1), radius=1) )
2316
2317        domain.compute_forcing_terms()
2318       
2319        assert allclose(domain.quantities['stage'].explicit_update[1], 2.0/pi)
2320        assert allclose(domain.quantities['stage'].explicit_update[0], 2.0/pi)
2321        assert allclose(domain.quantities['stage'].explicit_update[2:], 0)       
2322       
2323
2324
2325    #####################################################
2326    def test_first_order_extrapolator_const_z(self):
2327
2328        a = [0.0, 0.0]
2329        b = [0.0, 2.0]
2330        c = [2.0, 0.0]
2331        d = [0.0, 4.0]
2332        e = [2.0, 2.0]
2333        f = [4.0, 0.0]
2334
2335        points = [a, b, c, d, e, f]
2336        #bac, bce, ecf, dbe
2337        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
2338
2339        domain = Domain(points, vertices)
2340        val0 = 2.+2.0/3
2341        val1 = 4.+4.0/3
2342        val2 = 8.+2.0/3
2343        val3 = 2.+8.0/3
2344
2345        zl=zr=-3.75 #Assume constant bed (must be less than stage)
2346        domain.set_quantity('elevation', zl*ones( (4,3) ))
2347        domain.set_quantity('stage', [[val0, val0-1, val0-2],
2348                                      [val1, val1+1, val1],
2349                                      [val2, val2-2, val2],
2350                                      [val3-0.5, val3, val3]])
2351
2352
2353
2354        domain._order_ = 1
2355        domain.distribute_to_vertices_and_edges()
2356
2357        #Check that centroid values were distributed to vertices
2358        C = domain.quantities['stage'].centroid_values
2359        for i in range(3):
2360            assert allclose( domain.quantities['stage'].vertex_values[:,i], C)
2361
2362
2363    def test_first_order_limiter_variable_z(self):
2364        #Check that first order limiter follows bed_slope
2365        from Numeric import alltrue, greater_equal
2366        from anuga.config import epsilon
2367
2368        a = [0.0, 0.0]
2369        b = [0.0, 2.0]
2370        c = [2.0,0.0]
2371        d = [0.0, 4.0]
2372        e = [2.0, 2.0]
2373        f = [4.0,0.0]
2374
2375        points = [a, b, c, d, e, f]
2376        #bac, bce, ecf, dbe
2377        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
2378
2379        domain = Domain(points, vertices)
2380        val0 = 2.+2.0/3
2381        val1 = 4.+4.0/3
2382        val2 = 8.+2.0/3
2383        val3 = 2.+8.0/3
2384
2385        domain.set_quantity('elevation', [[0,0,0], [6,0,0],
2386                                          [6,6,6], [6,6,6]])
2387        domain.set_quantity('stage', [[val0, val0, val0],
2388                                      [val1, val1, val1],
2389                                      [val2, val2, val2],
2390                                      [val3, val3, val3]])
2391
2392        E = domain.quantities['elevation'].vertex_values
2393        L = domain.quantities['stage'].vertex_values
2394
2395
2396        #Check that some stages are not above elevation (within eps)
2397        #- so that the limiter has something to work with
2398        assert not alltrue(alltrue(greater_equal(L,E-epsilon)))
2399
2400        domain._order_ = 1
2401        domain.distribute_to_vertices_and_edges()
2402
2403        #Check that all stages are above elevation (within eps)
2404        assert alltrue(alltrue(greater_equal(L,E-epsilon)))
2405
2406
2407    #####################################################
2408    def test_distribute_basic(self):
2409        #Using test data generated by abstract_2d_finite_volumes-2
2410        #Assuming no friction and flat bed (0.0)
2411
2412        a = [0.0, 0.0]
2413        b = [0.0, 2.0]
2414        c = [2.0, 0.0]
2415        d = [0.0, 4.0]
2416        e = [2.0, 2.0]
2417        f = [4.0, 0.0]
2418
2419        points = [a, b, c, d, e, f]
2420        #bac, bce, ecf, dbe
2421        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
2422
2423        domain = Domain(points, vertices)
2424
2425        val0 = 2.
2426        val1 = 4.
2427        val2 = 8.
2428        val3 = 2.
2429
2430        domain.set_quantity('stage', [val0, val1, val2, val3],
2431                            location='centroids')
2432        L = domain.quantities['stage'].vertex_values
2433
2434        #First order
2435        domain._order_ = 1
2436        domain.distribute_to_vertices_and_edges()
2437        assert allclose(L[1], val1)
2438
2439        #Second order
2440        domain._order_ = 2
2441        domain.beta_w      = 0.9
2442        domain.beta_w_dry  = 0.9
2443        domain.beta_uh     = 0.9
2444        domain.beta_uh_dry = 0.9
2445        domain.beta_vh     = 0.9
2446        domain.beta_vh_dry = 0.9
2447        domain.distribute_to_vertices_and_edges()
2448        assert allclose(L[1], [2.2, 4.9, 4.9])
2449
2450
2451
2452    def test_distribute_away_from_bed(self):
2453        #Using test data generated by abstract_2d_finite_volumes-2
2454        #Assuming no friction and flat bed (0.0)
2455
2456        a = [0.0, 0.0]
2457        b = [0.0, 2.0]
2458        c = [2.0, 0.0]
2459        d = [0.0, 4.0]
2460        e = [2.0, 2.0]
2461        f = [4.0, 0.0]
2462
2463        points = [a, b, c, d, e, f]
2464        #bac, bce, ecf, dbe
2465        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
2466
2467        domain = Domain(points, vertices)
2468        L = domain.quantities['stage'].vertex_values
2469
2470        def stage(x,y):
2471            return x**2
2472
2473        domain.set_quantity('stage', stage, location='centroids')
2474
2475        domain.quantities['stage'].compute_gradients()
2476
2477        a, b = domain.quantities['stage'].get_gradients()
2478               
2479        assert allclose(a[1], 3.33333334)
2480        assert allclose(b[1], 0.0)
2481
2482        domain._order_ = 1
2483        domain.distribute_to_vertices_and_edges()
2484        assert allclose(L[1], 1.77777778)
2485
2486        domain._order_ = 2
2487        domain.beta_w      = 0.9
2488        domain.beta_w_dry  = 0.9
2489        domain.beta_uh     = 0.9
2490        domain.beta_uh_dry = 0.9
2491        domain.beta_vh     = 0.9
2492        domain.beta_vh_dry = 0.9
2493        domain.distribute_to_vertices_and_edges()
2494        assert allclose(L[1], [0.57777777, 2.37777778, 2.37777778])
2495
2496
2497
2498    def test_distribute_away_from_bed1(self):
2499        #Using test data generated by abstract_2d_finite_volumes-2
2500        #Assuming no friction and flat bed (0.0)
2501
2502        a = [0.0, 0.0]
2503        b = [0.0, 2.0]
2504        c = [2.0, 0.0]
2505        d = [0.0, 4.0]
2506        e = [2.0, 2.0]
2507        f = [4.0, 0.0]
2508
2509        points = [a, b, c, d, e, f]
2510        #bac, bce, ecf, dbe
2511        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
2512
2513        domain = Domain(points, vertices)
2514        L = domain.quantities['stage'].vertex_values
2515
2516        def stage(x,y):
2517            return x**4+y**2
2518
2519        domain.set_quantity('stage', stage, location='centroids')
2520        #print domain.quantities['stage'].centroid_values
2521
2522        domain.quantities['stage'].compute_gradients()
2523        a, b = domain.quantities['stage'].get_gradients()
2524        assert allclose(a[1], 25.18518519)
2525        assert allclose(b[1], 3.33333333)
2526
2527        domain._order_ = 1
2528        domain.distribute_to_vertices_and_edges()
2529        assert allclose(L[1], 4.9382716)
2530
2531        domain._order_ = 2
2532        domain.beta_w      = 0.9
2533        domain.beta_w_dry  = 0.9
2534        domain.beta_uh     = 0.9
2535        domain.beta_uh_dry = 0.9
2536        domain.beta_vh     = 0.9
2537        domain.beta_vh_dry = 0.9
2538        domain.distribute_to_vertices_and_edges()
2539        assert allclose(L[1], [1.07160494, 6.46058131, 7.28262855])
2540
2541
2542
2543    def test_distribute_near_bed(self):
2544
2545        a = [0.0, 0.0]
2546        b = [0.0, 2.0]
2547        c = [2.0, 0.0]
2548        d = [0.0, 4.0]
2549        e = [2.0, 2.0]
2550        f = [4.0, 0.0]
2551
2552        points = [a, b, c, d, e, f]
2553        #bac, bce, ecf, dbe
2554        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
2555
2556        domain = Domain(points, vertices)
2557
2558
2559        #Set up for a gradient of (10,0) at mid triangle (bce)
2560        def slope(x, y):
2561            return 10*x
2562
2563        h = 0.1
2564        def stage(x, y):
2565            return slope(x, y) + h
2566
2567        domain.set_quantity('elevation', slope)
2568        domain.set_quantity('stage', stage, location='centroids')
2569
2570        #print domain.quantities['elevation'].centroid_values
2571        #print domain.quantities['stage'].centroid_values
2572
2573        E = domain.quantities['elevation'].vertex_values
2574        L = domain.quantities['stage'].vertex_values
2575
2576        # Get reference values
2577        volumes = []
2578        for i in range(len(L)):
2579            volumes.append(sum(L[i])/3)
2580            assert allclose(volumes[i], domain.quantities['stage'].centroid_values[i]) 
2581       
2582       
2583        domain._order_ = 1
2584       
2585        domain.tight_slope_limiters = 0
2586        domain.distribute_to_vertices_and_edges()
2587        assert allclose(L[1], [0.1, 20.1, 20.1])
2588        for i in range(len(L)):
2589            assert allclose(volumes[i], sum(L[i])/3)                   
2590       
2591        domain.tight_slope_limiters = 1 # Allow triangle to be flatter (closer to bed)
2592        domain.distribute_to_vertices_and_edges()
2593        assert allclose(L[1], [0.298, 20.001, 20.001])
2594        for i in range(len(L)):
2595            assert allclose(volumes[i], sum(L[i])/3)   
2596
2597        domain._order_ = 2
2598       
2599        domain.tight_slope_limiters = 0
2600        domain.distribute_to_vertices_and_edges()
2601        assert allclose(L[1], [0.1, 20.1, 20.1])       
2602        for i in range(len(L)):
2603            assert allclose(volumes[i], sum(L[i])/3)           
2604       
2605        domain.tight_slope_limiters = 1 # Allow triangle to be flatter (closer to bed)
2606        domain.distribute_to_vertices_and_edges()
2607        assert allclose(L[1], [0.298, 20.001, 20.001])
2608        for i in range(len(L)):
2609            assert allclose(volumes[i], sum(L[i])/3)   
2610       
2611
2612
2613    def test_distribute_near_bed1(self):
2614
2615        a = [0.0, 0.0]
2616        b = [0.0, 2.0]
2617        c = [2.0, 0.0]
2618        d = [0.0, 4.0]
2619        e = [2.0, 2.0]
2620        f = [4.0, 0.0]
2621
2622        points = [a, b, c, d, e, f]
2623        #bac, bce, ecf, dbe
2624        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
2625
2626        domain = Domain(points, vertices)
2627
2628
2629        #Set up for a gradient of (8,2) at mid triangle (bce)
2630        def slope(x, y):
2631            return x**4+y**2
2632
2633        h = 0.1
2634        def stage(x,y):
2635            return slope(x,y)+h
2636
2637        domain.set_quantity('elevation', slope)
2638        domain.set_quantity('stage', stage)
2639
2640        #print domain.quantities['elevation'].centroid_values
2641        #print domain.quantities['stage'].centroid_values
2642
2643        E = domain.quantities['elevation'].vertex_values
2644        L = domain.quantities['stage'].vertex_values
2645
2646        # Get reference values
2647        volumes = []
2648        for i in range(len(L)):
2649            volumes.append(sum(L[i])/3)
2650            assert allclose(volumes[i], domain.quantities['stage'].centroid_values[i]) 
2651       
2652        #print E
2653        domain._order_ = 1
2654       
2655        domain.tight_slope_limiters = 0
2656        domain.distribute_to_vertices_and_edges()
2657        assert allclose(L[1], [4.1, 16.1, 20.1])       
2658        for i in range(len(L)):
2659            assert allclose(volumes[i], sum(L[i])/3)
2660       
2661               
2662        domain.tight_slope_limiters = 1 # Allow triangle to be flatter (closer to bed)
2663        domain.distribute_to_vertices_and_edges()
2664        assert allclose(L[1], [4.2386, 16.0604, 20.001])
2665        for i in range(len(L)):
2666            assert allclose(volumes[i], sum(L[i])/3)   
2667       
2668
2669        domain._order_ = 2
2670       
2671        domain.tight_slope_limiters = 0   
2672        domain.distribute_to_vertices_and_edges()
2673        assert allclose(L[1], [4.1, 16.1, 20.1])
2674        for i in range(len(L)):
2675            assert allclose(volumes[i], sum(L[i])/3)   
2676       
2677        domain.tight_slope_limiters = 1 # Allow triangle to be flatter (closer to bed)
2678        domain.distribute_to_vertices_and_edges()
2679        #print L[1]
2680        assert allclose(L[1], [4.23370103, 16.06529897, 20.001]) or\
2681               allclose(L[1], [4.18944138, 16.10955862, 20.001]) or\
2682               allclose(L[1], [4.19351461, 16.10548539, 20.001]) # old limiters
2683       
2684        for i in range(len(L)):
2685            assert allclose(volumes[i], sum(L[i])/3)
2686
2687
2688    def test_second_order_distribute_real_data(self):
2689        #Using test data generated by abstract_2d_finite_volumes-2
2690        #Assuming no friction and flat bed (0.0)
2691
2692        a = [0.0, 0.0]
2693        b = [0.0, 1.0/5]
2694        c = [0.0, 2.0/5]
2695        d = [1.0/5, 0.0]
2696        e = [1.0/5, 1.0/5]
2697        f = [1.0/5, 2.0/5]
2698        g = [2.0/5, 2.0/5]
2699
2700        points = [a, b, c, d, e, f, g]
2701        #bae, efb, cbf, feg
2702        vertices = [ [1,0,4], [4,5,1], [2,1,5], [5,4,6]]
2703
2704        domain = Domain(points, vertices)
2705
2706        def slope(x, y):
2707            return -x/3
2708
2709        domain.set_quantity('elevation', slope)
2710        domain.set_quantity('stage',
2711                            [0.01298164, 0.00365611,
2712                             0.01440365, -0.0381856437096],
2713                            location='centroids')
2714        domain.set_quantity('xmomentum',
2715                            [0.00670439, 0.01263789,
2716                             0.00647805, 0.0178180740668],
2717                            location='centroids')
2718        domain.set_quantity('ymomentum',
2719                            [-7.23510980e-004, -6.30413883e-005,
2720                             6.30413883e-005, 0.000200907255866],
2721                            location='centroids')
2722
2723        E = domain.quantities['elevation'].vertex_values
2724        L = domain.quantities['stage'].vertex_values
2725        X = domain.quantities['xmomentum'].vertex_values
2726        Y = domain.quantities['ymomentum'].vertex_values
2727
2728        #print E
2729        domain._order_ = 2
2730        domain.beta_w      = 0.9
2731        domain.beta_w_dry  = 0.9
2732        domain.beta_uh     = 0.9
2733        domain.beta_uh_dry = 0.9
2734        domain.beta_vh     = 0.9
2735        domain.beta_vh_dry = 0.9
2736        domain.beta_h = 0.0 #Use first order in h-limiter
2737       
2738        # FIXME (Ole): Need tests where this is commented out
2739        domain.tight_slope_limiters = 0 # Backwards compatibility (14/4/7)                 
2740        domain.use_centroid_velocities = 0 # Backwards compatibility (7/5/8)
2741       
2742               
2743        domain.distribute_to_vertices_and_edges()
2744
2745        #print L[1,:]
2746        #print X[1,:]
2747        #print Y[1,:]
2748
2749        assert allclose(L[1,:], [-0.00825735775384,
2750                                 -0.00801881482869,
2751                                 0.0272445025825])
2752        assert allclose(X[1,:], [0.0143507718962,
2753                                 0.0142502147066,
2754                                 0.00931268339717])
2755        assert allclose(Y[1,:], [-0.000117062180693,
2756                                 7.94434448109e-005,
2757                                 -0.000151505429018])
2758
2759
2760
2761    def test_balance_deep_and_shallow(self):
2762        """Test that balanced limiters preserve conserved quantites.
2763        This test is using old depth based balanced limiters
2764        """
2765        import copy
2766
2767        a = [0.0, 0.0]
2768        b = [0.0, 2.0]
2769        c = [2.0, 0.0]
2770        d = [0.0, 4.0]
2771        e = [2.0, 2.0]
2772        f = [4.0, 0.0]
2773
2774        points = [a, b, c, d, e, f]
2775
2776        #bac, bce, ecf, dbe
2777        elements = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4] ]
2778
2779        domain = Domain(points, elements)
2780        domain.check_integrity()
2781
2782        #Create a deliberate overshoot
2783        domain.set_quantity('stage', [[3,0,3], [2,2,6], [5,3,8], [8,3,5]])
2784        domain.set_quantity('elevation', 0) #Flat bed
2785        stage = domain.quantities['stage']
2786
2787        ref_centroid_values = copy.copy(stage.centroid_values[:]) #Copy
2788
2789        #Limit
2790        domain.tight_slope_limiters = 0               
2791        domain.distribute_to_vertices_and_edges()
2792
2793        #Assert that quantities are conserved
2794        from Numeric import sum
2795        for k in range(len(domain)):
2796            assert allclose (ref_centroid_values[k],
2797                             sum(stage.vertex_values[k,:])/3)
2798
2799
2800        #Now try with a non-flat bed - closely hugging initial stage in places
2801        #This will create alphas in the range [0, 0.478260, 1]
2802        domain.set_quantity('stage', [[3,0,3], [2,2,6], [5,3,8], [8,3,5]])
2803        domain.set_quantity('elevation', [[0,0,0],
2804                                        [1.8,1.9,5.9],
2805                                        [4.6,0,0],
2806                                        [0,2,4]])
2807        stage = domain.quantities['stage']
2808
2809        ref_centroid_values = copy.copy(stage.centroid_values[:]) #Copy
2810        ref_vertex_values = copy.copy(stage.vertex_values[:]) #Copy
2811
2812        #Limit
2813        domain.tight_slope_limiters = 0       
2814        domain.distribute_to_vertices_and_edges()
2815
2816
2817        #Assert that all vertex quantities have changed
2818        for k in range(len(domain)):
2819            #print ref_vertex_values[k,:], stage.vertex_values[k,:]
2820            assert not allclose (ref_vertex_values[k,:], stage.vertex_values[k,:])
2821        #and assert that quantities are still conserved
2822        from Numeric import sum
2823        for k in range(len(domain)):
2824            assert allclose (ref_centroid_values[k],
2825                             sum(stage.vertex_values[k,:])/3)
2826
2827
2828        # Check actual results
2829        assert allclose (stage.vertex_values,
2830                         [[2,2,2],
2831                          [1.93333333, 2.03333333, 6.03333333],
2832                          [6.93333333, 4.53333333, 4.53333333],
2833                          [5.33333333, 5.33333333, 5.33333333]])
2834
2835
2836    def test_balance_deep_and_shallow_tight_SL(self):
2837        """Test that balanced limiters preserve conserved quantites.
2838        This test is using Tight Slope Limiters
2839        """
2840        import copy
2841
2842        a = [0.0, 0.0]
2843        b = [0.0, 2.0]
2844        c = [2.0, 0.0]
2845        d = [0.0, 4.0]
2846        e = [2.0, 2.0]
2847        f = [4.0, 0.0]
2848
2849        points = [a, b, c, d, e, f]
2850
2851        #bac, bce, ecf, dbe
2852        elements = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4] ]
2853
2854        domain = Domain(points, elements)
2855        domain.check_integrity()
2856
2857        #Create a deliberate overshoot
2858        domain.set_quantity('stage', [[3,0,3], [2,2,6], [5,3,8], [8,3,5]])
2859        domain.set_quantity('elevation', 0) #Flat bed
2860        stage = domain.quantities['stage']
2861
2862        ref_centroid_values = copy.copy(stage.centroid_values[:]) #Copy
2863
2864        #Limit
2865        domain.tight_slope_limiters = 1               
2866        domain.distribute_to_vertices_and_edges()
2867
2868        #Assert that quantities are conserved
2869        from Numeric import sum
2870        for k in range(len(domain)):
2871            assert allclose (ref_centroid_values[k],
2872                             sum(stage.vertex_values[k,:])/3)
2873
2874
2875        #Now try with a non-flat bed - closely hugging initial stage in places
2876        #This will create alphas in the range [0, 0.478260, 1]
2877        domain.set_quantity('stage', [[3,0,3], [2,2,6], [5,3,8], [8,3,5]])
2878        domain.set_quantity('elevation', [[0,0,0],
2879                                        [1.8,1.9,5.9],
2880                                        [4.6,0,0],
2881                                        [0,2,4]])
2882        stage = domain.quantities['stage']
2883
2884        ref_centroid_values = copy.copy(stage.centroid_values[:]) #Copy
2885        ref_vertex_values = copy.copy(stage.vertex_values[:]) #Copy
2886
2887        #Limit
2888        domain.tight_slope_limiters = 1       
2889        domain.distribute_to_vertices_and_edges()
2890
2891
2892        #Assert that all vertex quantities have changed
2893        for k in range(len(domain)):
2894            #print ref_vertex_values[k,:], stage.vertex_values[k,:]
2895            assert not allclose (ref_vertex_values[k,:], stage.vertex_values[k,:])
2896        #and assert that quantities are still conserved
2897        from Numeric import sum
2898        for k in range(len(domain)):
2899            assert allclose (ref_centroid_values[k],
2900                             sum(stage.vertex_values[k,:])/3)
2901
2902
2903        #Also check that Python and C version produce the same
2904        # No longer applicable if tight_slope_limiters == 1
2905        #print stage.vertex_values
2906        #assert allclose (stage.vertex_values,
2907        #                 [[2,2,2],
2908        #                  [1.93333333, 2.03333333, 6.03333333],
2909        #                  [6.93333333, 4.53333333, 4.53333333],
2910        #                  [5.33333333, 5.33333333, 5.33333333]])
2911
2912
2913
2914    def test_balance_deep_and_shallow_Froude(self):
2915        """Test that balanced limiters preserve conserved quantites -
2916        and also that excessive Froude numbers are dealt with.
2917        This test is using tight slope limiters.
2918        """
2919        import copy
2920        from Numeric import sqrt, absolute
2921
2922        a = [0.0, 0.0]
2923        b = [0.0, 2.0]
2924        c = [2.0, 0.0]
2925        d = [0.0, 4.0]
2926        e = [2.0, 2.0]
2927        f = [4.0, 0.0]
2928
2929        points = [a, b, c, d, e, f]
2930
2931        # bac, bce, ecf, dbe
2932        elements = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4] ]
2933
2934        domain = Domain(points, elements)
2935        domain.check_integrity()
2936        domain.tight_slope_limiters = True
2937        domain.use_centroid_velocities = True               
2938
2939        # Create non-flat bed - closely hugging initial stage in places
2940        # This will create alphas in the range [0, 0.478260, 1]
2941        domain.set_quantity('stage', [[3,0,3], [2,2,6], [5,3,8], [8,3,5]])
2942        domain.set_quantity('elevation', [[0,0,0],
2943                                        [1.8,1.999,5.999],
2944                                        [4.6,0,0],
2945                                        [0,2,4]])
2946
2947        # Create small momenta, that nonetheless will generate large speeds
2948        # due to shallow depth at isolated vertices
2949        domain.set_quantity('xmomentum', -0.0058)
2950        domain.set_quantity('ymomentum', 0.0890)
2951
2952
2953
2954       
2955        stage = domain.quantities['stage']
2956        elevation = domain.quantities['elevation']
2957        xmomentum = domain.quantities['xmomentum']
2958        ymomentum = domain.quantities['ymomentum']       
2959
2960        # Setup triangle #1 to mimick real Froude explosion observed
2961        # in the Onslow example 13 Nov 2007.
2962
2963        stage.vertex_values[1,:] = [1.6385, 1.6361, 1.2953]
2964        elevation.vertex_values[1,:] = [1.6375, 1.6336, 0.4647]       
2965        xmomentum.vertex_values[1,:] = [-0.0058, -0.0050, -0.0066]
2966        ymomentum.vertex_values[1,:] = [0.0890, 0.0890, 0.0890]
2967
2968        xmomentum.interpolate()
2969        ymomentum.interpolate()       
2970        stage.interpolate()       
2971        elevation.interpolate()
2972
2973        # Verify interpolation
2974        assert allclose(stage.centroid_values[1], 1.5233)
2975        assert allclose(elevation.centroid_values[1], 1.2452667)
2976        assert allclose(xmomentum.centroid_values[1], -0.0058)
2977        assert allclose(ymomentum.centroid_values[1], 0.089)
2978
2979        # Derived quantities
2980        depth = stage-elevation
2981        u = xmomentum/depth
2982        v = ymomentum/depth
2983
2984        denom = (depth*g)**0.5 
2985        Fx = u/denom
2986        Fy = v/denom
2987       
2988   
2989        # Verify against Onslow example (14 Nov 2007)
2990        assert allclose(depth.centroid_values[1], 0.278033)
2991        assert allclose(u.centroid_values[1], -0.0208608)
2992        assert allclose(v.centroid_values[1], 0.3201055)
2993
2994        assert allclose(denom.centroid_values[1],
2995                        sqrt(depth.centroid_values[1]*g))
2996
2997        assert allclose(u.centroid_values[1]/denom.centroid_values[1],
2998                        -0.012637746977)
2999        assert allclose(Fx.centroid_values[1],
3000                        u.centroid_values[1]/denom.centroid_values[1])
3001
3002        # Check that Froude numbers are small at centroids.
3003        assert allclose(Fx.centroid_values[1], -0.012637746977)
3004        assert allclose(Fy.centroid_values[1], 0.193924048435)
3005
3006
3007        # But Froude numbers are huge at some vertices and edges
3008        assert allclose(Fx.vertex_values[1,:], [-5.85888475e+01,
3009                                                -1.27775313e+01,
3010                                                -2.78511420e-03])
3011
3012        assert allclose(Fx.edge_values[1,:], [-6.89150773e-03,
3013                                              -7.38672488e-03,
3014                                              -2.35626238e+01])
3015
3016        assert allclose(Fy.vertex_values[1,:], [8.99035764e+02,
3017                                                2.27440057e+02,
3018                                                3.75568430e-02])
3019
3020        assert allclose(Fy.edge_values[1,:], [1.05748998e-01,
3021                                              1.06035244e-01,
3022                                              3.88346947e+02])
3023       
3024       
3025        # The task is now to arrange the limiters such that Froude numbers
3026        # remain under control whil at the same time obeying the conservation
3027        # laws.
3028
3029       
3030        ref_centroid_values = copy.copy(stage.centroid_values[:]) #Copy
3031        ref_vertex_values = copy.copy(stage.vertex_values[:]) #Copy
3032
3033        # Limit (and invoke balance_deep_and_shallow)
3034        domain.tight_slope_limiters = 1
3035        domain.distribute_to_vertices_and_edges()
3036
3037        # Redo derived quantities
3038        depth = stage-elevation
3039        u = xmomentum/depth
3040        v = ymomentum/depth
3041
3042        # Assert that all vertex velocities stay within one
3043        # order of magnitude of centroid velocities.
3044        #print u.vertex_values[1,:]
3045        #print u.centroid_values[1]
3046       
3047        assert alltrue(absolute(u.vertex_values[1,:]) <= absolute(u.centroid_values[1])*10)
3048        assert alltrue(absolute(v.vertex_values[1,:]) <= absolute(v.centroid_values[1])*10) 
3049
3050        denom = (depth*g)**0.5 
3051        Fx = u/denom
3052        Fy = v/denom
3053
3054
3055        # Assert that Froude numbers are less than max value (TBA)
3056        # at vertices, edges and centroids.
3057        from anuga.config import maximum_froude_number
3058        assert alltrue(absolute(Fx.vertex_values[1,:]) < maximum_froude_number)
3059        assert alltrue(absolute(Fy.vertex_values[1,:]) < maximum_froude_number)
3060
3061
3062        # Assert that all vertex quantities have changed
3063        for k in range(len(domain)):
3064            #print ref_vertex_values[k,:], stage.vertex_values[k,:]
3065            assert not allclose (ref_vertex_values[k,:],
3066                                 stage.vertex_values[k,:])
3067           
3068        # Assert that quantities are still conserved
3069        from Numeric import sum
3070        for k in range(len(domain)):
3071            assert allclose (ref_centroid_values[k],
3072                             sum(stage.vertex_values[k,:])/3)
3073
3074
3075       
3076        return
3077   
3078        qwidth = 12
3079        for k in [1]: #range(len(domain)):
3080            print 'Triangle %d (C, V, E)' %k
3081           
3082            print 'stage'.ljust(qwidth), stage.centroid_values[k],\
3083                  stage.vertex_values[k,:], stage.edge_values[k,:]
3084            print 'elevation'.ljust(qwidth), elevation.centroid_values[k],\
3085                  elevation.vertex_values[k,:], elevation.edge_values[k,:]
3086            print 'depth'.ljust(qwidth), depth.centroid_values[k],\
3087                  depth.vertex_values[k,:], depth.edge_values[k,:]
3088            print 'xmomentum'.ljust(qwidth), xmomentum.centroid_values[k],\
3089                  xmomentum.vertex_values[k,:], xmomentum.edge_values[k,:]
3090            print 'ymomentum'.ljust(qwidth), ymomentum.centroid_values[k],\
3091                  ymomentum.vertex_values[k,:], ymomentum.edge_values[k,:]
3092            print 'u'.ljust(qwidth),u.centroid_values[k],\
3093                  u.vertex_values[k,:], u.edge_values[k,:]
3094            print 'v'.ljust(qwidth), v.centroid_values[k],\
3095                  v.vertex_values[k,:], v.edge_values[k,:]
3096            print 'Fx'.ljust(qwidth), Fx.centroid_values[k],\
3097                  Fx.vertex_values[k,:], Fx.edge_values[k,:]
3098            print 'Fy'.ljust(qwidth), Fy.centroid_values[k],\
3099                  Fy.vertex_values[k,:], Fy.edge_values[k,:]
3100           
3101           
3102       
3103
3104
3105
3106    def test_conservation_1(self):
3107        """Test that stage is conserved globally
3108
3109        This one uses a flat bed, reflective bdries and a suitable
3110        initial condition
3111        """
3112        from mesh_factory import rectangular
3113        from Numeric import array
3114
3115        #Create basic mesh
3116        points, vertices, boundary = rectangular(6, 6)
3117
3118        #Create shallow water domain
3119        domain = Domain(points, vertices, boundary)
3120        domain.smooth = False
3121        domain.default_order=2
3122
3123        #IC
3124        def x_slope(x, y):
3125            return x/3
3126
3127        domain.set_quantity('elevation', 0)
3128        domain.set_quantity('friction', 0)
3129        domain.set_quantity('stage', x_slope)
3130
3131        # Boundary conditions (reflective everywhere)
3132        Br = Reflective_boundary(domain)
3133        domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
3134
3135        domain.check_integrity()
3136
3137        initial_volume = domain.quantities['stage'].get_integral()
3138        initial_xmom = domain.quantities['xmomentum'].get_integral()
3139
3140        #print initial_xmom
3141
3142        #Evolution
3143        for t in domain.evolve(yieldstep = 0.05, finaltime = 5.0):
3144            volume =  domain.quantities['stage'].get_integral()
3145            assert allclose (volume, initial_volume)
3146
3147            #I don't believe that the total momentum should be the same
3148            #It starts with zero and ends with zero though
3149            #xmom = domain.quantities['xmomentum'].get_integral()
3150            #print xmom
3151            #assert allclose (xmom, initial_xmom)
3152
3153        os.remove(domain.get_name() + '.sww')
3154
3155
3156    def test_conservation_2(self):
3157        """Test that stage is conserved globally
3158
3159        This one uses a slopy bed, reflective bdries and a suitable
3160        initial condition
3161        """
3162        from mesh_factory import rectangular
3163        from Numeric import array
3164
3165        #Create basic mesh
3166        points, vertices, boundary = rectangular(6, 6)
3167
3168        #Create shallow water domain
3169        domain = Domain(points, vertices, boundary)
3170        domain.smooth = False
3171        domain.default_order=2
3172
3173        #IC
3174        def x_slope(x, y):
3175            return x/3
3176
3177        domain.set_quantity('elevation', x_slope)
3178        domain.set_quantity('friction', 0)
3179        domain.set_quantity('stage', 0.4) #Steady
3180
3181        # Boundary conditions (reflective everywhere)
3182        Br = Reflective_boundary(domain)
3183        domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
3184
3185        domain.check_integrity()
3186
3187        initial_volume = domain.quantities['stage'].get_integral()
3188        initial_xmom = domain.quantities['xmomentum'].get_integral()
3189
3190        #print initial_xmom
3191
3192        #Evolution
3193        for t in domain.evolve(yieldstep = 0.05, finaltime = 5.0):
3194            volume =  domain.quantities['stage'].get_integral()
3195            assert allclose (volume, initial_volume)
3196
3197            #FIXME: What would we expect from momentum
3198            #xmom = domain.quantities['xmomentum'].get_integral()
3199            #print xmom
3200            #assert allclose (xmom, initial_xmom)
3201
3202        os.remove(domain.get_name() + '.sww')
3203
3204    def test_conservation_3(self):
3205        """Test that stage is conserved globally
3206
3207        This one uses a larger grid, convoluted bed, reflective bdries and a suitable
3208        initial condition
3209        """
3210        from mesh_factory import rectangular
3211        from Numeric import array
3212
3213        #Create basic mesh
3214        points, vertices, boundary = rectangular(2, 1)
3215
3216        #Create shallow water domain
3217        domain = Domain(points, vertices, boundary)
3218        domain.smooth = False
3219        domain.default_order = 2
3220        domain.beta_h = 0.0
3221        domain.set_quantities_to_be_stored(['stage', 'xmomentum', 'ymomentum'])
3222
3223        #IC
3224        def x_slope(x, y):
3225            z = 0*x
3226            for i in range(len(x)):
3227                if x[i] < 0.3:
3228                    z[i] = x[i]/3
3229                if 0.3 <= x[i] < 0.5:
3230                    z[i] = -0.5
3231                if 0.5 <= x[i] < 0.7:
3232                    z[i] = 0.39
3233                if 0.7 <= x[i]:
3234                    z[i] = x[i]/3
3235            return z
3236
3237
3238
3239        domain.set_quantity('elevation', x_slope)
3240        domain.set_quantity('friction', 0)
3241        domain.set_quantity('stage', 0.4) #Steady
3242
3243        # Boundary conditions (reflective everywhere)
3244        Br = Reflective_boundary(domain)
3245        domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
3246
3247        domain.check_integrity()
3248
3249        initial_volume = domain.quantities['stage'].get_integral()
3250        initial_xmom = domain.quantities['xmomentum'].get_integral()
3251
3252        import copy
3253        ref_centroid_values =\
3254                 copy.copy(domain.quantities['stage'].centroid_values)
3255
3256        #print 'ORG', domain.quantities['stage'].centroid_values
3257        domain.distribute_to_vertices_and_edges()
3258
3259
3260        #print domain.quantities['stage'].centroid_values
3261        assert allclose(domain.quantities['stage'].centroid_values,
3262                        ref_centroid_values)
3263
3264
3265        #Check that initial limiter doesn't violate cons quan
3266        assert allclose(domain.quantities['stage'].get_integral(),
3267                        initial_volume)
3268
3269        #Evolution
3270        for t in domain.evolve(yieldstep = 0.05, finaltime = 10):
3271            volume =  domain.quantities['stage'].get_integral()
3272            #print t, volume, initial_volume
3273            assert allclose (volume, initial_volume)
3274
3275        os.remove(domain.get_name() + '.sww')
3276
3277    def test_conservation_4(self):
3278        """Test that stage is conserved globally
3279
3280        This one uses a larger grid, convoluted bed, reflective bdries and a suitable
3281        initial condition
3282        """
3283        from mesh_factory import rectangular
3284        from Numeric import array
3285
3286        #Create basic mesh
3287        points, vertices, boundary = rectangular(6, 6)
3288
3289        #Create shallow water domain
3290        domain = Domain(points, vertices, boundary)
3291        domain.smooth = False
3292        domain.default_order=2
3293        domain.beta_h = 0.0
3294
3295        #IC
3296        def x_slope(x, y):
3297            z = 0*x
3298            for i in range(len(x)):
3299                if x[i] < 0.3:
3300                    z[i] = x[i]/3
3301                if 0.3 <= x[i] < 0.5:
3302                    z[i] = -0.5
3303                if 0.5 <= x[i] < 0.7:
3304                    #z[i] = 0.3 #OK with beta == 0.2
3305                    z[i] = 0.34 #OK with beta == 0.0
3306                    #z[i] = 0.35#Fails after 80 timesteps with an error
3307                                #of the order 1.0e-5
3308                if 0.7 <= x[i]:
3309                    z[i] = x[i]/3
3310            return z
3311
3312
3313
3314        domain.set_quantity('elevation', x_slope)
3315        domain.set_quantity('friction', 0)
3316        domain.set_quantity('stage', 0.4) #Steady
3317
3318        # Boundary conditions (reflective everywhere)
3319        Br = Reflective_boundary(domain)
3320        domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
3321
3322        domain.check_integrity()
3323
3324        initial_volume = domain.quantities['stage'].get_integral()
3325        initial_xmom = domain.quantities['xmomentum'].get_integral()
3326
3327        import copy
3328        ref_centroid_values =\
3329                 copy.copy(domain.quantities['stage'].centroid_values)
3330
3331        #Test limiter by itself
3332        domain.distribute_to_vertices_and_edges()
3333
3334        #Check that initial limiter doesn't violate cons quan
3335        assert allclose (domain.quantities['stage'].get_integral(),
3336                         initial_volume)
3337        #NOTE: This would fail if any initial stage was less than the
3338        #corresponding bed elevation - but that is reasonable.
3339
3340
3341        #Evolution
3342        for t in domain.evolve(yieldstep = 0.05, finaltime = 10.0):
3343            volume =  domain.quantities['stage'].get_integral()
3344
3345            #print t, volume, initial_volume
3346
3347            assert allclose (volume, initial_volume)
3348
3349
3350        os.remove(domain.get_name() + '.sww')
3351
3352
3353    def test_conservation_5(self):
3354        """Test that momentum is conserved globally in
3355        steady state scenario
3356
3357        This one uses a slopy bed, dirichlet and reflective bdries
3358        """
3359        from mesh_factory import rectangular
3360        from Numeric import array
3361
3362        # Create basic mesh
3363        points, vertices, boundary = rectangular(6, 6)
3364
3365        # Create shallow water domain
3366        domain = Domain(points, vertices, boundary)
3367        domain.smooth = False
3368        domain.default_order = 2
3369
3370        # IC
3371        def x_slope(x, y):
3372            return x/3
3373
3374        domain.set_quantity('elevation', x_slope)
3375        domain.set_quantity('friction', 0)
3376        domain.set_quantity('stage', 0.4) # Steady
3377
3378        # Boundary conditions (reflective everywhere)
3379        Br = Reflective_boundary(domain)
3380        Bleft = Dirichlet_boundary([0.5,0,0])
3381        Bright = Dirichlet_boundary([0.1,0,0])
3382        domain.set_boundary({'left': Bleft, 'right': Bright,
3383                             'top': Br, 'bottom': Br})
3384
3385        domain.check_integrity()
3386
3387        initial_volume = domain.quantities['stage'].get_integral()
3388        initial_xmom = domain.quantities['xmomentum'].get_integral()
3389
3390
3391        # Evolution
3392        for t in domain.evolve(yieldstep = 0.05, finaltime = 15.0):
3393            stage =  domain.quantities['stage'].get_integral()
3394            xmom = domain.quantities['xmomentum'].get_integral()
3395            ymom = domain.quantities['ymomentum'].get_integral()
3396
3397            if allclose(t, 6):  # Steady state reached
3398                steady_xmom = domain.quantities['xmomentum'].get_integral()
3399                steady_ymom = domain.quantities['ymomentum'].get_integral()
3400                steady_stage = domain.quantities['stage'].get_integral()
3401
3402            if t > 6:
3403                #print '%.2f %14.8f %14.8f' %(t, ymom, steady_ymom)
3404                msg = 'xmom=%.2f, steady_xmom=%.2f' %(xmom, steady_xmom)
3405                assert allclose(xmom, steady_xmom), msg
3406                assert allclose(ymom, steady_ymom)
3407                assert allclose(stage, steady_stage)
3408
3409
3410        os.remove(domain.get_name() + '.sww')
3411
3412
3413
3414
3415
3416    def test_conservation_real(self):
3417        """Test that momentum is conserved globally
3418       
3419        Stephen finally made a test that revealed the problem.
3420        This test failed with code prior to 25 July 2005
3421        """
3422
3423        yieldstep = 0.01
3424        finaltime = 0.05
3425        min_depth = 1.0e-2
3426
3427
3428        import sys
3429        from os import sep; sys.path.append('..'+sep+'abstract_2d_finite_volumes')
3430        from mesh_factory import rectangular
3431
3432
3433        #Create shallow water domain
3434        points, vertices, boundary = rectangular(10, 10, len1=500, len2=500)
3435        domain = Domain(points, vertices, boundary)
3436        domain.smooth = False
3437        domain.default_order = 1
3438        domain.minimum_allowed_height = min_depth
3439
3440        # Set initial condition
3441        class Set_IC:
3442            """Set an initial condition with a constant value, for x0<x<x1
3443            """
3444
3445            def __init__(self, x0=0.25, x1=0.5, h=1.0):
3446                self.x0 = x0
3447                self.x1 = x1
3448                self.= h
3449
3450            def __call__(self, x, y):
3451                return self.h*((x>self.x0)&(x<self.x1))
3452
3453
3454        domain.set_quantity('stage', Set_IC(200.0,300.0,5.0))
3455
3456
3457        #Boundaries
3458        R = Reflective_boundary(domain)
3459        domain.set_boundary( {'left': R, 'right': R, 'top':R, 'bottom': R})
3460
3461        ref = domain.quantities['stage'].get_integral()
3462
3463        # Evolution
3464        for t in domain.evolve(yieldstep = yieldstep, finaltime = finaltime):
3465            pass
3466            #print 'Integral stage = ',\
3467            #      domain.quantities['stage'].get_integral(),\
3468            #      ' Time = ',domain.time
3469
3470
3471        now = domain.quantities['stage'].get_integral()
3472
3473        msg = 'Stage not conserved: was %f, now %f' %(ref, now)
3474        assert allclose(ref, now), msg
3475
3476        os.remove(domain.get_name() + '.sww')
3477
3478    def test_second_order_flat_bed_onestep(self):
3479
3480        from mesh_factory import rectangular
3481        from Numeric import array
3482
3483        #Create basic mesh
3484        points, vertices, boundary = rectangular(6, 6)
3485
3486        #Create shallow water domain
3487        domain = Domain(points, vertices, boundary)
3488        domain.smooth = False
3489        domain.default_order = 2
3490        domain.beta_w      = 0.9
3491        domain.beta_w_dry  = 0.9
3492        domain.beta_uh     = 0.9
3493        domain.beta_uh_dry = 0.9
3494        domain.beta_vh     = 0.9
3495        domain.beta_vh_dry = 0.9
3496        domain.H0 = 0 # Backwards compatibility (6/2/7)       
3497       
3498        # Boundary conditions
3499        Br = Reflective_boundary(domain)
3500        Bd = Dirichlet_boundary([0.1, 0., 0.])
3501        domain.set_boundary({'left': Bd, 'right': Br, 'top': Br, 'bottom': Br})
3502
3503        domain.check_integrity()
3504
3505        # Evolution
3506        for t in domain.evolve(yieldstep = 0.05, finaltime = 0.05):
3507            pass# domain.write_time()
3508
3509        # Data from earlier version of abstract_2d_finite_volumes
3510        assert allclose(domain.min_timestep, 0.0396825396825)
3511        assert allclose(domain.max_timestep, 0.0396825396825)
3512
3513        assert allclose(domain.quantities['stage'].centroid_values[:12],
3514                        [0.00171396, 0.02656103, 0.00241523, 0.02656103,
3515                        0.00241523, 0.02656103, 0.00241523, 0.02656103,
3516                        0.00241523, 0.02656103, 0.00241523, 0.0272623])
3517
3518        domain.distribute_to_vertices_and_edges()
3519
3520        assert allclose(domain.quantities['stage'].vertex_values[:12,0],
3521                        [0.0001714, 0.02656103, 0.00024152,
3522                        0.02656103, 0.00024152, 0.02656103,
3523                        0.00024152, 0.02656103, 0.00024152,
3524                        0.02656103, 0.00024152, 0.0272623])
3525
3526        assert allclose(domain.quantities['stage'].vertex_values[:12,1],
3527                        [0.00315012, 0.02656103, 0.00024152, 0.02656103,
3528                         0.00024152, 0.02656103, 0.00024152, 0.02656103,
3529                         0.00024152, 0.02656103, 0.00040506, 0.0272623])
3530
3531        assert allclose(domain.quantities['stage'].vertex_values[:12,2],
3532                        [0.00182037, 0.02656103, 0.00676264,
3533                         0.02656103, 0.00676264, 0.02656103,
3534                         0.00676264, 0.02656103, 0.00676264,
3535                         0.02656103, 0.0065991, 0.0272623])
3536
3537        assert allclose(domain.quantities['xmomentum'].centroid_values[:12],
3538                        [0.00113961, 0.01302432, 0.00148672,
3539                         0.01302432, 0.00148672, 0.01302432,
3540                         0.00148672, 0.01302432, 0.00148672 ,
3541                         0.01302432, 0.00148672, 0.01337143])
3542
3543        assert allclose(domain.quantities['ymomentum'].centroid_values[:12],
3544                        [-2.91240050e-004 , 1.22721531e-004,
3545                         -1.22721531e-004, 1.22721531e-004 ,
3546                         -1.22721531e-004, 1.22721531e-004,
3547                         -1.22721531e-004 , 1.22721531e-004,
3548                         -1.22721531e-004, 1.22721531e-004,
3549                         -1.22721531e-004, -4.57969873e-005])
3550
3551        os.remove(domain.get_name() + '.sww')
3552
3553
3554    def test_second_order_flat_bed_moresteps(self):
3555
3556        from mesh_factory import rectangular
3557        from Numeric import array
3558
3559        #Create basic mesh
3560        points, vertices, boundary = rectangular(6, 6)
3561
3562        #Create shallow water domain
3563        domain = Domain(points, vertices, boundary)
3564        domain.smooth = False
3565        domain.default_order=2
3566
3567        # Boundary conditions
3568        Br = Reflective_boundary(domain)
3569        Bd = Dirichlet_boundary([0.1, 0., 0.])
3570        domain.set_boundary({'left': Bd, 'right': Br, 'top': Br, 'bottom': Br})
3571
3572        domain.check_integrity()
3573
3574        #Evolution
3575        for t in domain.evolve(yieldstep = 0.05, finaltime = 0.1):
3576            pass
3577
3578        #Data from earlier version of abstract_2d_finite_volumes
3579        #assert allclose(domain.min_timestep, 0.0396825396825)
3580        #assert allclose(domain.max_timestep, 0.0396825396825)
3581        #print domain.quantities['stage'].centroid_values
3582
3583        os.remove(domain.get_name() + '.sww')       
3584
3585
3586    def test_flatbed_first_order(self):
3587        from mesh_factory import rectangular
3588        from Numeric import array
3589
3590        #Create basic mesh
3591        N = 8
3592        points, vertices, boundary = rectangular(N, N)
3593
3594        #Create shallow water domain
3595        domain = Domain(points, vertices, boundary)
3596        domain.smooth = False
3597        domain.default_order=1
3598        domain.H0 = 0 # Backwards compatibility (6/2/7)       
3599
3600        # Boundary conditions
3601        Br = Reflective_boundary(domain)
3602        Bd = Dirichlet_boundary([0.2,0.,0.])
3603
3604        domain.set_boundary({'left': Bd, 'right': Br, 'top': Br, 'bottom': Br})
3605        domain.check_integrity()
3606
3607
3608        #Evolution
3609        for t in domain.evolve(yieldstep = 0.02, finaltime = 0.5):
3610            pass
3611            #domain.write_time()
3612
3613        #FIXME: These numbers were from version before 25/10
3614        #assert allclose(domain.min_timestep, 0.0140413643926)
3615        #assert allclose(domain.max_timestep, 0.0140947355753)
3616
3617        for i in range(3):
3618            #assert allclose(domain.quantities['stage'].edge_values[:4,i],
3619            #                [0.10730244,0.12337617,0.11200126,0.12605666])
3620
3621            assert allclose(domain.quantities['xmomentum'].edge_values[:4,i],
3622                            [0.07610894,0.06901572,0.07284461,0.06819712])
3623
3624            #assert allclose(domain.quantities['ymomentum'].edge_values[:4,i],
3625            #                [-0.0060238, -0.00157404, -0.00309633, -0.0001637])
3626
3627
3628        os.remove(domain.get_name() + '.sww')           
3629
3630    def test_flatbed_second_order(self):
3631        from mesh_factory import rectangular
3632        from Numeric import array
3633
3634        #Create basic mesh
3635        N = 8
3636        points, vertices, boundary = rectangular(N, N)
3637
3638        #Create shallow water domain
3639        domain = Domain(points, vertices, boundary)
3640        domain.smooth = False
3641        domain.default_order=2
3642        domain.beta_w      = 0.9
3643        domain.beta_w_dry  = 0.9
3644        domain.beta_uh     = 0.9
3645        domain.beta_uh_dry = 0.9
3646        domain.beta_vh     = 0.9
3647        domain.beta_vh_dry = 0.9       
3648        #domain.minimum_allowed_height = 0.0 #Makes it like the 'oldstyle' balance
3649        domain.H0 = 0 # Backwards compatibility (6/2/7)
3650        domain.use_centroid_velocities = False # Backwards compatibility (8/5/8)
3651        domain.set_maximum_allowed_speed(1.0)       
3652
3653        # Boundary conditions
3654        Br = Reflective_boundary(domain)
3655        Bd = Dirichlet_boundary([0.2,0.,0.])
3656
3657        domain.set_boundary({'left': Bd, 'right': Br, 'top': Br, 'bottom': Br})
3658        domain.check_integrity()
3659
3660        # Evolution
3661        for t in domain.evolve(yieldstep = 0.01, finaltime = 0.03):
3662            pass
3663
3664        msg = 'min step was %f instead of %f' %(domain.min_timestep,
3665                                                0.0210448446782) 
3666
3667        assert allclose(domain.min_timestep, 0.0210448446782), msg
3668        assert allclose(domain.max_timestep, 0.0210448446782)
3669
3670        #print domain.quantities['stage'].vertex_values[:4,0]
3671        #print domain.quantities['xmomentum'].vertex_values[:4,0]
3672        #print domain.quantities['ymomentum'].vertex_values[:4,0]
3673
3674        #FIXME: These numbers were from version before 25/10
3675        #assert allclose(domain.quantities['stage'].vertex_values[:4,0],
3676        #                [0.00101913,0.05352143,0.00104852,0.05354394])
3677
3678        #FIXME: These numbers were from version before 21/3/6 -
3679        #could be recreated by setting maximum_allowed_speed to 0 maybe
3680        #assert allclose(domain.quantities['xmomentum'].vertex_values[:4,0],
3681        #                [ 0.00064835, 0.03685719, 0.00085073, 0.03687313])
3682       
3683        assert allclose(domain.quantities['xmomentum'].vertex_values[:4,0],
3684                        [ 0.00090581, 0.03685719, 0.00088303, 0.03687313])
3685                       
3686                       
3687
3688        #assert allclose(domain.quantities['xmomentum'].vertex_values[:4,0],
3689        #                [0.00090581,0.03685719,0.00088303,0.03687313])
3690
3691        assert allclose(domain.quantities['ymomentum'].vertex_values[:4,0],
3692                        [-0.00139463,0.0006156,-0.00060364,0.00061827])
3693
3694
3695        os.remove(domain.get_name() + '.sww')
3696
3697       
3698    def test_flatbed_second_order_vmax_0(self):
3699        from mesh_factory import rectangular
3700        from Numeric import array
3701
3702        #Create basic mesh
3703        N = 8
3704        points, vertices, boundary = rectangular(N, N)
3705
3706        #Create shallow water domain
3707        domain = Domain(points, vertices, boundary)
3708        domain.smooth = False
3709        domain.default_order=2
3710        domain.beta_w      = 0.9
3711        domain.beta_w_dry  = 0.9
3712        domain.beta_uh     = 0.9
3713        domain.beta_uh_dry = 0.9
3714        domain.beta_vh     = 0.9
3715        domain.beta_vh_dry = 0.9       
3716        domain.maximum_allowed_speed = 0.0 #Makes it like the 'oldstyle'
3717        domain.H0 = 0 # Backwards compatibility (6/2/7)
3718        domain.use_centroid_velocities = False # Backwards compatibility (8/5/8)       
3719
3720        # Boundary conditions
3721        Br = Reflective_boundary(domain)
3722        Bd = Dirichlet_boundary([0.2,0.,0.])
3723
3724        domain.set_boundary({'left': Bd, 'right': Br, 'top': Br, 'bottom': Br})
3725        domain.check_integrity()
3726
3727        #Evolution
3728        for t in domain.evolve(yieldstep = 0.01, finaltime = 0.03):
3729            pass
3730
3731
3732        assert allclose(domain.min_timestep, 0.0210448446782)
3733        assert allclose(domain.max_timestep, 0.0210448446782)
3734
3735        #FIXME: These numbers were from version before 21/3/6 -
3736        #could be recreated by setting maximum_allowed_speed to 0 maybe
3737        assert allclose(domain.quantities['xmomentum'].vertex_values[:4,0],
3738                        [ 0.00064835, 0.03685719, 0.00085073, 0.03687313]) 
3739       
3740
3741        assert allclose(domain.quantities['ymomentum'].vertex_values[:4,0],
3742                        [-0.00139463,0.0006156,-0.00060364,0.00061827])
3743
3744
3745        os.remove(domain.get_name() + '.sww')
3746
3747       
3748
3749    def test_flatbed_second_order_distribute(self):
3750        #Use real data from anuga.abstract_2d_finite_volumes 2
3751        #painfully setup and extracted.
3752        from mesh_factory import rectangular
3753        from Numeric import array
3754
3755        #Create basic mesh
3756        N = 8
3757        points, vertices, boundary = rectangular(N, N)
3758
3759        #Create shallow water domain
3760        domain = Domain(points, vertices, boundary)
3761        domain.smooth = False
3762        domain.default_order=domain._order_=2
3763        domain.beta_w      = 0.9
3764        domain.beta_w_dry  = 0.9
3765        domain.beta_uh     = 0.9
3766        domain.beta_uh_dry = 0.9
3767        domain.beta_vh     = 0.9
3768        domain.beta_vh_dry = 0.9
3769        domain.H0 = 0 # Backwards compatibility (6/2/7)
3770        domain.use_centroid_velocities = False # Backwards compatibility (8/5/8)       
3771        domain.set_maximum_allowed_speed(1.0)       
3772
3773        # Boundary conditions
3774        Br = Reflective_boundary(domain)
3775        Bd = Dirichlet_boundary([0.2,0.,0.])
3776
3777        domain.set_boundary({'left': Bd, 'right': Br, 'top': Br, 'bottom': Br})
3778        domain.check_integrity()
3779
3780
3781
3782        for V in [False, True]:
3783            if V:
3784                #Set centroids as if system had been evolved
3785                L = zeros(2*N*N, Float)
3786                L[:32] = [7.21205592e-003, 5.35214298e-002, 1.00910824e-002,
3787                          5.35439433e-002, 1.00910824e-002, 5.35439433e-002,
3788                          1.00910824e-002, 5.35439433e-002, 1.00910824e-002,
3789                          5.35439433e-002, 1.00910824e-002, 5.35439433e-002,
3790                          1.00910824e-002, 5.35393928e-002, 1.02344264e-002,
3791                          5.59605058e-002, 0.00000000e+000, 3.31027800e-004,
3792                          0.00000000e+000, 4.37962142e-005, 0.00000000e+000,
3793                          4.37962142e-005, 0.00000000e+000, 4.37962142e-005,
3794                          0.00000000e+000, 4.37962142e-005, 0.00000000e+000,
3795                          4.37962142e-005, 0.00000000e+000, 4.37962142e-005,
3796                          0.00000000e+000, 5.57305948e-005]
3797
3798                X = zeros(2*N*N, Float)
3799                X[:32] = [6.48351607e-003, 3.68571894e-002, 8.50733285e-003,
3800                          3.68731327e-002, 8.50733285e-003, 3.68731327e-002,
3801                          8.50733285e-003, 3.68731327e-002, 8.50733285e-003,
3802                          3.68731327e-002, 8.50733285e-003, 3.68731327e-002,
3803                          8.50733285e-003, 3.68693861e-002, 8.65220973e-003,
3804                          3.85055387e-002, 0.00000000e+000, 2.86060840e-004,
3805                          0.00000000e+000, 3.58905503e-005, 0.00000000e+000,
3806                          3.58905503e-005, 0.00000000e+000, 3.58905503e-005,
3807                          0.00000000e+000, 3.58905503e-005, 0.00000000e+000,
3808                          3.58905503e-005, 0.00000000e+000, 3.58905503e-005,
3809                          0.00000000e+000, 4.57662812e-005]
3810
3811                Y = zeros(2*N*N, Float)
3812                Y[:32]=[-1.39463104e-003, 6.15600298e-004, -6.03637382e-004,
3813                        6.18272251e-004, -6.03637382e-004, 6.18272251e-004,
3814                        -6.03637382e-004, 6.18272251e-004, -6.03637382e-004,
3815                        6.18272251e-004, -6.03637382e-004, 6.18272251e-004,
3816                        -6.03637382e-004, 6.18599320e-004, -6.74622797e-004,
3817                        -1.48934756e-004, 0.00000000e+000, -5.35079969e-005,
3818                        0.00000000e+000, -2.57264987e-005, 0.00000000e+000,
3819                        -2.57264987e-005, 0.00000000e+000, -2.57264987e-005,
3820                        0.00000000e+000, -2.57264987e-005, 0.00000000e+000,
3821                        -2.57264987e-005, 0.00000000e+000, -2.57264987e-005,
3822                        0.00000000e+000, -2.57635178e-005]
3823
3824
3825                domain.set_quantity('stage', L, location='centroids')
3826                domain.set_quantity('xmomentum', X, location='centroids')
3827                domain.set_quantity('ymomentum', Y, location='centroids')
3828
3829                domain.check_integrity()
3830            else:
3831                #Evolution
3832                for t in domain.evolve(yieldstep = 0.01, finaltime = 0.03):
3833                    pass
3834                assert allclose(domain.min_timestep, 0.0210448446782)
3835                assert allclose(domain.max_timestep, 0.0210448446782)
3836
3837
3838            #Centroids were correct but not vertices.
3839            #Hence the check of distribute below.
3840            assert allclose(domain.quantities['stage'].centroid_values[:4],
3841                            [0.00721206,0.05352143,0.01009108,0.05354394])
3842
3843            assert allclose(domain.quantities['xmomentum'].centroid_values[:4],
3844                            [0.00648352,0.03685719,0.00850733,0.03687313])
3845
3846            assert allclose(domain.quantities['ymomentum'].centroid_values[:4],
3847                            [-0.00139463,0.0006156,-0.00060364,0.00061827])
3848
3849            #print 'C17=', domain.quantities['xmomentum'].centroid_values[17]
3850            #print 'C19=', domain.quantities['xmomentum'].centroid_values[19]
3851
3852            #assert allclose(domain.quantities['xmomentum'].centroid_values[17],0.00028606084)
3853            ##print domain.quantities['xmomentum'].centroid_values[17], V
3854            ##print
3855            if not V:
3856                #FIXME: These numbers were from version before 21/3/6 -
3857                #could be recreated by setting maximum_allowed_speed to 0 maybe           
3858               
3859                #assert allclose(domain.quantities['xmomentum'].centroid_values[17], 0.0)               
3860                assert allclose(domain.quantities['xmomentum'].centroid_values[17], 0.000286060839592)                         
3861
3862            else:
3863                assert allclose(domain.quantities['xmomentum'].centroid_values[17], 0.00028606084)
3864
3865            import copy
3866            XX = copy.copy(domain.quantities['xmomentum'].centroid_values)
3867            assert allclose(domain.quantities['xmomentum'].centroid_values, XX)
3868
3869            domain.distribute_to_vertices_and_edges()
3870
3871            #assert allclose(domain.quantities['xmomentum'].centroid_values, XX)
3872
3873            #assert allclose(domain.quantities['xmomentum'].centroid_values[17],
3874            #                0.0)
3875            assert allclose(domain.quantities['xmomentum'].centroid_values[17], 0.000286060839592)                                                 
3876
3877
3878            #FIXME: These numbers were from version before 25/10
3879            #assert allclose(domain.quantities['stage'].vertex_values[:4,0],
3880            #                [0.00101913,0.05352143,0.00104852,0.05354394])
3881
3882            assert allclose(domain.quantities['ymomentum'].vertex_values[:4,0],
3883                            [-0.00139463,0.0006156,-0.00060364,0.00061827])
3884
3885
3886            assert allclose(domain.quantities['xmomentum'].vertex_values[:4,0],
3887                            [0.00090581,0.03685719,0.00088303,0.03687313])
3888
3889
3890            #NB NO longer relvant:
3891
3892            #This was the culprit. First triangles vertex 0 had an
3893            #x-momentum of 0.0064835 instead of 0.00090581 and
3894            #third triangle had 0.00850733 instead of 0.00088303
3895            #print domain.quantities['xmomentum'].vertex_values[:4,0]
3896
3897            #print domain.quantities['xmomentum'].vertex_values[:4,0]
3898            #assert allclose(domain.quantities['xmomentum'].vertex_values[:4,0],
3899            #                [0.00090581,0.03685719,0.00088303,0.03687313])
3900
3901        os.remove(domain.get_name() + '.sww')
3902
3903
3904
3905    def test_bedslope_problem_first_order(self):
3906
3907        from mesh_factory import rectangular
3908        from Numeric import array
3909
3910        #Create basic mesh
3911        points, vertices, boundary = rectangular(6, 6)
3912
3913        #Create shallow water domain
3914        domain = Domain(points, vertices, boundary)
3915        domain.smooth = False
3916        domain.default_order = 1
3917
3918        #Bed-slope and friction
3919        def x_slope(x, y):
3920            return -x/3
3921
3922        domain.set_quantity('elevation', x_slope)
3923
3924        # Boundary conditions
3925        Br = Reflective_boundary(domain)
3926        domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
3927
3928        #Initial condition
3929        #domain.set_quantity('stage', Constant_height(x_slope, 0.05))
3930        domain.set_quantity('stage', expression='elevation+0.05')
3931        domain.check_integrity()
3932
3933        #Evolution
3934        for t in domain.evolve(yieldstep = 0.05, finaltime = 0.05):
3935            pass# domain.write_time()
3936
3937        # FIXME (Ole): Need some other assertion here!
3938        #print domain.min_timestep, domain.max_timestep   
3939        #assert allclose(domain.min_timestep, 0.050010003001)
3940        #assert allclose(domain.max_timestep, 0.050010003001)
3941
3942
3943        os.remove(domain.get_name() + '.sww')
3944       
3945    def test_bedslope_problem_first_order_moresteps(self):
3946
3947        from mesh_factory import rectangular
3948        from Numeric import array
3949
3950        #Create basic mesh
3951        points, vertices, boundary = rectangular(6, 6)
3952
3953        #Create shallow water domain
3954        domain = Domain(points, vertices, boundary)
3955        domain.smooth = False
3956        domain.default_order = 1
3957        domain.beta_h = 0.0 # Use first order in h-limiter
3958       
3959        # FIXME (Ole): Need tests where these two are commented out
3960        domain.H0 = 0        # Backwards compatibility (6/2/7)       
3961        domain.tight_slope_limiters = 0 # Backwards compatibility (14/4/7)
3962        domain.use_centroid_velocities = 0 # Backwards compatibility (7/5/8)                 
3963
3964        #Bed-slope and friction
3965        def x_slope(x, y):
3966            return -x/3
3967
3968        domain.set_quantity('elevation', x_slope)
3969
3970        # Boundary conditions
3971        Br = Reflective_boundary(domain)
3972        domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
3973
3974        #Initial condition
3975        domain.set_quantity('stage', expression='elevation+0.05')
3976        domain.check_integrity()
3977
3978        #Evolution
3979        for t in domain.evolve(yieldstep = 0.05, finaltime = 0.5):
3980            pass# domain.write_time()
3981
3982        #Data from earlier version of abstract_2d_finite_volumes
3983        #print domain.quantities['stage'].centroid_values
3984
3985        assert allclose(domain.quantities['stage'].centroid_values,
3986                        [-0.02998628, -0.01520652, -0.03043492,
3987                        -0.0149132, -0.03004706, -0.01476251,
3988                        -0.0298215, -0.01467976, -0.02988158,
3989                        -0.01474662, -0.03036161, -0.01442995,
3990                        -0.07624583, -0.06297061, -0.07733792,
3991                        -0.06342237, -0.07695439, -0.06289595,
3992                        -0.07635559, -0.0626065, -0.07633628,
3993                        -0.06280072, -0.07739632, -0.06386738,
3994                        -0.12161738, -0.11028239, -0.1223796,
3995                        -0.11095953, -0.12189744, -0.11048616,
3996                        -0.12074535, -0.10987605, -0.12014311,
3997                        -0.10976691, -0.12096859, -0.11087692,
3998                        -0.16868259, -0.15868061, -0.16801135,
3999                        -0.1588003, -0.16674343, -0.15813323,
4000                        -0.16457595, -0.15693826, -0.16281096,
4001                        -0.15585154, -0.16283873, -0.15540068,
4002                        -0.17450362, -0.19919913, -0.18062882,
4003                        -0.19764131, -0.17783111, -0.19407213,
4004                        -0.1736915, -0.19053624, -0.17228678,
4005                        -0.19105634, -0.17920133, -0.1968828,
4006                        -0.14244395, -0.14604641, -0.14473537,
4007                        -0.1506107, -0.14510055, -0.14919522,
4008                        -0.14175896, -0.14560798, -0.13911658,
4009                        -0.14439383, -0.13924047, -0.14829043])
4010
4011        os.remove(domain.get_name() + '.sww')
4012       
4013    def test_bedslope_problem_second_order_one_step(self):
4014
4015        from mesh_factory import rectangular
4016        from Numeric import array
4017
4018        #Create basic mesh
4019        points, vertices, boundary = rectangular(6, 6)
4020
4021        #Create shallow water domain
4022        domain = Domain(points, vertices, boundary)
4023        domain.smooth = False
4024        domain.default_order=2
4025        domain.beta_w      = 0.9
4026        domain.beta_w_dry  = 0.9
4027        domain.beta_uh     = 0.9
4028        domain.beta_uh_dry = 0.9
4029        domain.beta_vh     = 0.9
4030        domain.beta_vh_dry = 0.9
4031
4032       
4033        # FIXME (Ole): Need tests where this is commented out
4034        domain.tight_slope_limiters = 0 # Backwards compatibility (14/4/7)
4035        domain.use_centroid_velocities = 0 # Backwards compatibility (7/5/8)     
4036       
4037        #Bed-slope and friction at vertices (and interpolated elsewhere)
4038        def x_slope(x, y):
4039            return -x/3
4040
4041        domain.set_quantity('elevation', x_slope)
4042
4043        # Boundary conditions
4044        Br = Reflective_boundary(domain)
4045        domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
4046
4047        #Initial condition
4048        domain.set_quantity('stage', expression='elevation+0.05')
4049        domain.check_integrity()
4050
4051        assert allclose(domain.quantities['stage'].centroid_values,
4052                        [0.01296296, 0.03148148, 0.01296296,
4053                        0.03148148, 0.01296296, 0.03148148,
4054                        0.01296296, 0.03148148, 0.01296296,
4055                        0.03148148, 0.01296296, 0.03148148,
4056                        -0.04259259, -0.02407407, -0.04259259,
4057                        -0.02407407, -0.04259259, -0.02407407,
4058                        -0.04259259, -0.02407407, -0.04259259,
4059                        -0.02407407, -0.04259259, -0.02407407,
4060                        -0.09814815, -0.07962963, -0.09814815,
4061                        -0.07962963, -0.09814815, -0.07962963,
4062                        -0.09814815, -0.07962963, -0.09814815,
4063                        -0.07962963, -0.09814815, -0.07962963,
4064                        -0.1537037 , -0.13518519, -0.1537037,
4065                        -0.13518519, -0.1537037, -0.13518519,
4066                        -0.1537037 , -0.13518519, -0.1537037,
4067                        -0.13518519, -0.1537037, -0.13518519,
4068                        -0.20925926, -0.19074074, -0.20925926,
4069                        -0.19074074, -0.20925926, -0.19074074,
4070                        -0.20925926, -0.19074074, -0.20925926,
4071                        -0.19074074, -0.20925926, -0.19074074,
4072                        -0.26481481, -0.2462963, -0.26481481,
4073                        -0.2462963, -0.26481481, -0.2462963,
4074                        -0.26481481, -0.2462963, -0.26481481,
4075                        -0.2462963, -0.26481481, -0.2462963])
4076
4077
4078        #print domain.quantities['stage'].extrapolate_second_order()
4079        #domain.distribute_to_vertices_and_edges()
4080        #print domain.quantities['stage'].vertex_values[:,0]
4081
4082        #Evolution
4083        for t in domain.evolve(yieldstep = 0.05, finaltime = 0.05):
4084            #domain.write_time()
4085            pass
4086
4087
4088        #print domain.quantities['stage'].centroid_values
4089        assert allclose(domain.quantities['stage'].centroid_values,
4090                 [0.01290985, 0.02356019, 0.01619096, 0.02356019, 0.01619096,
4091                  0.02356019, 0.01619096, 0.02356019, 0.01619096, 0.02356019,
4092                  0.01619096, 0.0268413, -0.04411074, -0.0248011, -0.04186556,
4093                  -0.0248011, -0.04186556, -0.0248011, -0.04186556, -0.0248011,
4094                  -0.04186556, -0.0248011, -0.04186556, -0.02255593,
4095                  -0.09966629, -0.08035666, -0.09742112, -0.08035666,
4096                  -0.09742112, -0.08035666, -0.09742112, -0.08035666,
4097                  -0.09742112, -0.08035666, -0.09742112, -0.07811149,
4098                  -0.15522185, -0.13591222, -0.15297667, -0.13591222,
4099                  -0.15297667, -0.13591222, -0.15297667, -0.13591222,
4100                  -0.15297667, -0.13591222, -0.15297667, -0.13366704,
4101                  -0.2107774, -0.19146777, -0.20853223, -0.19146777,
4102                  -0.20853223, -0.19146777, -0.20853223, -0.19146777,
4103                  -0.20853223, -0.19146777, -0.20853223, -0.1892226,
4104                  -0.26120669, -0.24776246, -0.25865535, -0.24776246,
4105                  -0.25865535, -0.24776246, -0.25865535, -0.24776246,
4106                  -0.25865535, -0.24776246, -0.25865535, -0.24521113])
4107
4108        os.remove(domain.get_name() + '.sww')
4109
4110    def test_bedslope_problem_second_order_two_steps(self):
4111
4112        from mesh_factory import rectangular
4113        from Numeric import array
4114
4115        #Create basic mesh
4116        points, vertices, boundary = rectangular(6, 6)
4117
4118        #Create shallow water domain
4119        domain = Domain(points, vertices, boundary)
4120        domain.smooth = False
4121        domain.default_order=2
4122        domain.beta_w      = 0.9
4123        domain.beta_w_dry  = 0.9
4124        domain.beta_uh     = 0.9
4125        domain.beta_uh_dry = 0.9
4126        domain.beta_vh     = 0.9
4127        domain.beta_vh_dry = 0.9
4128        domain.beta_h = 0.0 #Use first order in h-limiter
4129       
4130        # FIXME (Ole): Need tests where this is commented out
4131        domain.tight_slope_limiters = 0 # Backwards compatibility (14/4/7)                 
4132        domain.H0 = 0 # Backwards compatibility (6/2/7)
4133        domain.use_centroid_velocities = 0 # Backwards compatibility (7/5/8)
4134       
4135
4136        #Bed-slope and friction at vertices (and interpolated elsewhere)
4137        def x_slope(x, y):
4138            return -x/3
4139
4140        domain.set_quantity('elevation', x_slope)
4141
4142        # Boundary conditions
4143        Br = Reflective_boundary(domain)
4144        domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
4145
4146        #Initial condition
4147        domain.set_quantity('stage', expression='elevation+0.05')
4148        domain.check_integrity()
4149
4150        assert allclose(domain.quantities['stage'].centroid_values,
4151                        [0.01296296, 0.03148148, 0.01296296,
4152                        0.03148148, 0.01296296, 0.03148148,
4153                        0.01296296, 0.03148148, 0.01296296,
4154                        0.03148148, 0.01296296, 0.03148148,
4155                        -0.04259259, -0.02407407, -0.04259259,
4156                        -0.02407407, -0.04259259, -0.02407407,
4157                        -0.04259259, -0.02407407, -0.04259259,
4158                        -0.02407407, -0.04259259, -0.02407407,
4159                        -0.09814815, -0.07962963, -0.09814815,
4160                        -0.07962963, -0.09814815, -0.07962963,
4161                        -0.09814815, -0.07962963, -0.09814815,
4162                        -0.07962963, -0.09814815, -0.07962963,
4163                        -0.1537037 , -0.13518519, -0.1537037,
4164                        -0.13518519, -0.1537037, -0.13518519,
4165                        -0.1537037 , -0.13518519, -0.1537037,
4166                        -0.13518519, -0.1537037, -0.13518519,
4167                        -0.20925926, -0.19074074, -0.20925926,
4168                        -0.19074074, -0.20925926, -0.19074074,
4169                        -0.20925926, -0.19074074, -0.20925926,
4170                        -0.19074074, -0.20925926, -0.19074074,
4171                        -0.26481481, -0.2462963, -0.26481481,
4172                        -0.2462963, -0.26481481, -0.2462963,
4173                        -0.26481481, -0.2462963, -0.26481481,
4174                        -0.2462963, -0.26481481, -0.2462963])
4175
4176
4177        #print domain.quantities['stage'].extrapolate_second_order()
4178        #domain.distribute_to_vertices_and_edges()
4179        #print domain.quantities['stage'].vertex_values[:,0]
4180
4181        #Evolution
4182        for t in domain.evolve(yieldstep = 0.05, finaltime = 0.1):
4183            pass
4184
4185
4186        #Data from earlier version of abstract_2d_finite_volumes ft=0.1
4187        assert allclose(domain.min_timestep, 0.0376895634803)
4188        assert allclose(domain.max_timestep, 0.0415635655309)
4189
4190
4191        assert allclose(domain.quantities['stage'].centroid_values,
4192                        [0.00855788, 0.01575204, 0.00994606, 0.01717072,
4193                         0.01005985, 0.01716362, 0.01005985, 0.01716299,
4194                         0.01007098, 0.01736248, 0.01216452, 0.02026776,
4195                         -0.04462374, -0.02479045, -0.04199789, -0.0229465,
4196                         -0.04184033, -0.02295693, -0.04184013, -0.02295675,
4197                         -0.04184486, -0.0228168, -0.04028876, -0.02036486,
4198                         -0.10029444, -0.08170809, -0.09772846, -0.08021704,
4199                         -0.09760006, -0.08022143, -0.09759984, -0.08022124,
4200                         -0.09760261, -0.08008893, -0.09603914, -0.07758209,
4201                         -0.15584152, -0.13723138, -0.15327266, -0.13572906,
4202                         -0.15314427, -0.13573349, -0.15314405, -0.13573331,
4203                         -0.15314679, -0.13560104, -0.15158523, -0.13310701,
4204                         -0.21208605, -0.19283913, -0.20955631, -0.19134189,
4205                         -0.20942821, -0.19134598, -0.20942799, -0.1913458,
4206                         -0.20943005, -0.19120952, -0.20781177, -0.18869401,
4207                         -0.25384082, -0.2463294, -0.25047649, -0.24464654,
4208                         -0.25031159, -0.24464253, -0.25031112, -0.24464253,
4209                         -0.25031463, -0.24454764, -0.24885323, -0.24286438])
4210
4211
4212        os.remove(domain.get_name() + '.sww')
4213
4214    def test_bedslope_problem_second_order_two_yieldsteps(self):
4215
4216        from mesh_factory import rectangular
4217        from Numeric import array
4218
4219        #Create basic mesh
4220        points, vertices, boundary = rectangular(6, 6)
4221
4222        #Create shallow water domain
4223        domain = Domain(points, vertices, boundary)
4224        domain.smooth = False
4225        domain.default_order=2
4226        domain.beta_w      = 0.9
4227        domain.beta_w_dry  = 0.9
4228        domain.beta_uh     = 0.9
4229        domain.beta_uh_dry = 0.9
4230        domain.beta_vh     = 0.9
4231        domain.beta_vh_dry = 0.9
4232        domain.beta_h = 0.0 #Use first order in h-limiter
4233       
4234        # FIXME (Ole): Need tests where this is commented out
4235        domain.tight_slope_limiters = 0 # Backwards compatibility (14/4/7)                 
4236        domain.H0 = 0 # Backwards compatibility (6/2/7)
4237        domain.use_centroid_velocities = 0 # Backwards compatibility (7/5/8)
4238       
4239
4240        #Bed-slope and friction at vertices (and interpolated elsewhere)
4241        def x_slope(x, y):
4242            return -x/3
4243
4244        domain.set_quantity('elevation', x_slope)
4245
4246        # Boundary conditions
4247        Br = Reflective_boundary(domain)
4248        domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
4249
4250        #Initial condition
4251        domain.set_quantity('stage', expression='elevation+0.05')
4252        domain.check_integrity()
4253
4254        assert allclose(domain.quantities['stage'].centroid_values,
4255                        [0.01296296, 0.03148148, 0.01296296,
4256                        0.03148148, 0.01296296, 0.03148148,
4257                        0.01296296, 0.03148148, 0.01296296,
4258                        0.03148148, 0.01296296, 0.03148148,
4259                        -0.04259259, -0.02407407, -0.04259259,
4260                        -0.02407407, -0.04259259, -0.02407407,
4261                        -0.04259259, -0.02407407, -0.04259259,
4262                        -0.02407407, -0.04259259, -0.02407407,
4263                        -0.09814815, -0.07962963, -0.09814815,
4264                        -0.07962963, -0.09814815, -0.07962963,
4265                        -0.09814815, -0.07962963, -0.09814815,
4266                        -0.07962963, -0.09814815, -0.07962963,
4267                        -0.1537037 , -0.13518519, -0.1537037,
4268                        -0.13518519, -0.1537037, -0.13518519,
4269                        -0.1537037 , -0.13518519, -0.1537037,
4270                        -0.13518519, -0.1537037, -0.13518519,
4271                        -0.20925926, -0.19074074, -0.20925926,
4272                        -0.19074074, -0.20925926, -0.19074074,
4273                        -0.20925926, -0.19074074, -0.20925926,
4274                        -0.19074074, -0.20925926, -0.19074074,
4275                        -0.26481481, -0.2462963, -0.26481481,
4276                        -0.2462963, -0.26481481, -0.2462963,
4277                        -0.26481481, -0.2462963, -0.26481481,
4278                        -0.2462963, -0.26481481, -0.2462963])
4279
4280
4281        #print domain.quantities['stage'].extrapolate_second_order()
4282        #domain.distribute_to_vertices_and_edges()
4283        #print domain.quantities['stage'].vertex_values[:,0]
4284
4285        #Evolution
4286        for t in domain.evolve(yieldstep = 0.05, finaltime = 0.1):   #0.05??
4287            #domain.write_time()
4288            pass
4289
4290
4291
4292        assert allclose(domain.quantities['stage'].centroid_values,
4293                 [0.00855788, 0.01575204, 0.00994606, 0.01717072, 0.01005985,
4294                  0.01716362, 0.01005985, 0.01716299, 0.01007098, 0.01736248,
4295                  0.01216452, 0.02026776, -0.04462374, -0.02479045, -0.04199789,
4296                  -0.0229465, -0.04184033, -0.02295693, -0.04184013,
4297                  -0.02295675, -0.04184486, -0.0228168, -0.04028876,
4298                  -0.02036486, -0.10029444, -0.08170809, -0.09772846,
4299                  -0.08021704, -0.09760006, -0.08022143, -0.09759984,
4300                  -0.08022124, -0.09760261, -0.08008893, -0.09603914,
4301                  -0.07758209, -0.15584152, -0.13723138, -0.15327266,
4302                  -0.13572906, -0.15314427, -0.13573349, -0.15314405,
4303                  -0.13573331, -0.15314679, -0.13560104, -0.15158523,
4304                  -0.13310701, -0.21208605, -0.19283913, -0.20955631,
4305                  -0.19134189, -0.20942821, -0.19134598, -0.20942799,
4306                  -0.1913458, -0.20943005, -0.19120952, -0.20781177,
4307                  -0.18869401, -0.25384082, -0.2463294, -0.25047649,
4308                  -0.24464654, -0.25031159, -0.24464253, -0.25031112,
4309                  -0.24464253, -0.25031463, -0.24454764, -0.24885323,
4310                  -0.24286438])
4311
4312        os.remove(domain.get_name() + '.sww')
4313
4314    def test_bedslope_problem_second_order_more_steps(self):
4315
4316        from mesh_factory import rectangular
4317        from Numeric import array
4318
4319        #Create basic mesh
4320        points, vertices, boundary = rectangular(6, 6)
4321
4322        #Create shallow water domain
4323        domain = Domain(points, vertices, boundary)
4324        domain.smooth = False
4325        domain.default_order=2
4326        domain.beta_w      = 0.9
4327        domain.beta_w_dry  = 0.9
4328        domain.beta_uh     = 0.9
4329        domain.beta_uh_dry = 0.9
4330        domain.beta_vh     = 0.9
4331        domain.beta_vh_dry = 0.9
4332        domain.beta_h = 0.0 #Use first order in h-limiter
4333       
4334       
4335        # FIXME (Ole): Need tests where these two are commented out
4336        domain.H0 = 0        # Backwards compatibility (6/2/7)       
4337        domain.tight_slope_limiters = 0 # Backwards compatibility (14/4/7)                 
4338        domain.use_centroid_velocities = 0 # Backwards compatibility (7/5/8)
4339       
4340               
4341
4342        #Bed-slope and friction at vertices (and interpolated elsewhere)
4343        def x_slope(x, y):
4344            return -x/3
4345
4346        domain.set_quantity('elevation', x_slope)
4347
4348        # Boundary conditions
4349        Br = Reflective_boundary(domain)
4350        domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
4351
4352        #Initial condition
4353        domain.set_quantity('stage', expression = 'elevation + 0.05')
4354        domain.check_integrity()
4355
4356        assert allclose(domain.quantities['stage'].centroid_values,
4357                        [0.01296296, 0.03148148, 0.01296296,
4358                        0.03148148, 0.01296296, 0.03148148,
4359                        0.01296296, 0.03148148, 0.01296296,
4360                        0.03148148, 0.01296296, 0.03148148,
4361                        -0.04259259, -0.02407407, -0.04259259,
4362                        -0.02407407, -0.04259259, -0.02407407,
4363                        -0.04259259, -0.02407407, -0.04259259,
4364                        -0.02407407, -0.04259259, -0.02407407,
4365                        -0.09814815, -0.07962963, -0.09814815,
4366                        -0.07962963, -0.09814815, -0.07962963,
4367                        -0.09814815, -0.07962963, -0.09814815,
4368                        -0.07962963, -0.09814815, -0.07962963,
4369                        -0.1537037 , -0.13518519, -0.1537037,
4370                        -0.13518519, -0.1537037, -0.13518519,
4371                        -0.1537037 , -0.13518519, -0.1537037,
4372                        -0.13518519, -0.1537037, -0.13518519,
4373                        -0.20925926, -0.19074074, -0.20925926,
4374                        -0.19074074, -0.20925926, -0.19074074,
4375                        -0.20925926, -0.19074074, -0.20925926,
4376                        -0.19074074, -0.20925926, -0.19074074,
4377                        -0.26481481, -0.2462963, -0.26481481,
4378                        -0.2462963, -0.26481481, -0.2462963,
4379                        -0.26481481, -0.2462963, -0.26481481,
4380                        -0.2462963, -0.26481481, -0.2462963])
4381
4382
4383        #print domain.quantities['stage'].extrapolate_second_order()
4384        #domain.distribute_to_vertices_and_edges()
4385        #print domain.quantities['stage'].vertex_values[:,0]
4386
4387        #Evolution
4388        for t in domain.evolve(yieldstep = 0.05, finaltime = 0.5):
4389
4390            # Check that diagnostics works
4391            msg = domain.timestepping_statistics(track_speeds=True)
4392            #FIXME(Ole): One might check the contents of msg here.
4393
4394
4395
4396        assert allclose(domain.quantities['stage'].centroid_values,
4397     [-0.02907028, -0.01475478, -0.02973417, -0.01447186, -0.02932665, -0.01428285,
4398      -0.02901975, -0.0141361,  -0.02898816, -0.01418135, -0.02961409, -0.01403487,
4399      -0.07597998, -0.06252591, -0.07664854, -0.06312532, -0.07638287, -0.06265139,
4400      -0.07571145, -0.06235231, -0.0756817,  -0.06245309, -0.07652292, -0.06289946,
4401      -0.12367464, -0.11088981, -0.12237277, -0.11115338, -0.1218934,  -0.1107174,
4402      -0.12081485, -0.11000491, -0.12038451, -0.11010335, -0.12102113, -0.11012105,
4403      -0.16909116, -0.15831543, -0.16730214, -0.15786249, -0.1665493,  -0.15697919,
4404      -0.16496618, -0.15559852, -0.16338679, -0.15509088, -0.16364092, -0.15424423,
4405      -0.18771107, -0.19903904, -0.18903759, -0.19858437, -0.18701552, -0.19697797,
4406      -0.1833593,  -0.19505871, -0.1818806,  -0.19418042, -0.18586159, -0.19576946,
4407      -0.13986873, -0.14170053, -0.14132188, -0.14560674, -0.14095617, -0.14373292,
4408      -0.13785933, -0.14033364, -0.13592955, -0.13936356, -0.13596008, -0.14216296])
4409
4410        assert allclose(domain.quantities['xmomentum'].centroid_values,
4411     [ 0.00831121,  0.00317948,  0.00731797,  0.00334939,  0.00764717,  0.00348053,
4412       0.00788729,  0.00356522,  0.00780649,  0.00341919,  0.00693525,  0.00310375,
4413       0.02166196,  0.01421475,  0.02017737,  0.01316839,  0.02037015,  0.01368659,
4414       0.02106,     0.01399161,  0.02074514,  0.01354935,  0.01887407,  0.0123113,
4415       0.03775083,  0.02855197,  0.03689337,  0.02759782,  0.03732848,  0.02812072,
4416       0.03872545,  0.02913348,  0.03880939,  0.02803804,  0.03546499,  0.0260039,
4417       0.0632131,   0.04730634,  0.0576324,   0.04592336,  0.05790921,  0.04690514,
4418       0.05986467,  0.04871165,  0.06170068,  0.04811572,  0.05657041,  0.04416292,
4419       0.08489642,  0.07188097,  0.07835261,  0.06843406,  0.07986412,  0.0698247,
4420       0.08201071,  0.07216756,  0.08378418,  0.07273624,  0.080399,    0.06645841,
4421       0.01631548,  0.04691608,  0.0206632,   0.044409,    0.02115518,  0.04560305,
4422       0.02160608,  0.04663725,  0.02174734,  0.04795559,  0.02281427,  0.05667111])
4423
4424
4425        assert allclose(domain.quantities['ymomentum'].centroid_values,
4426     [ 1.45876601e-004, -3.24627393e-004, -1.57572719e-004, -2.92790187e-004,
4427      -9.90988382e-005, -3.06677335e-004, -1.62493106e-004, -3.71310004e-004,
4428      -1.99445058e-004, -3.28493467e-004,  6.68217349e-005, -8.42042805e-006,
4429       5.05093371e-004, -1.42842214e-004, -6.81454718e-005, -5.02084057e-004,
4430      -8.50583861e-005, -4.65443981e-004, -1.96406564e-004, -5.88889562e-004,
4431      -2.70160173e-004, -5.35485454e-004,  2.60780997e-004,  3.12145471e-005,
4432       5.16189608e-004,  1.07069062e-004,  9.29989252e-005, -3.71211119e-004,
4433       1.16350246e-004, -3.82407830e-004, -1.62077969e-004, -6.30906636e-004,
4434      -4.74025708e-004, -6.94463009e-004,  6.15092843e-005,  2.22106820e-004,
4435      -6.29589294e-004,  2.43611937e-004, -5.88125094e-004, -6.94293192e-005,
4436      -4.17914641e-004,  6.64609019e-005, -7.68334577e-004, -3.40232101e-004,
4437      -1.67424308e-003, -7.39485066e-004, -1.59966988e-003,  5.68262838e-005,
4438      -1.48470633e-003, -1.84554882e-003, -2.27200099e-003, -1.67506848e-003,
4439      -1.95610258e-003, -1.47638801e-003, -1.73779477e-003, -1.85498791e-003,
4440      -2.01357843e-003, -2.17675471e-003, -1.65783870e-003, -1.15818681e-003,
4441      -1.18663036e-003, -2.94229849e-003, -3.59309018e-003, -5.13496584e-003,
4442      -6.17359400e-003, -5.98761937e-003, -6.00540116e-003, -5.01121966e-003,
4443      -4.50964850e-003, -3.06319963e-003,  6.08950810e-004, -4.79537921e-004])
4444
4445        os.remove(domain.get_name() + '.sww')
4446
4447
4448
4449    def NOtest_bedslope_problem_second_order_more_steps_feb_2007(self):
4450        """test_bedslope_problem_second_order_more_steps_feb_2007
4451
4452        Test shallow water finite volumes, using parameters from
4453        feb 2007 rather than backward compatibility ad infinitum
4454       
4455        """
4456        from mesh_factory import rectangular
4457        from Numeric import array
4458
4459        #Create basic mesh
4460        points, vertices, boundary = rectangular(6, 6)
4461
4462        #Create shallow water domain
4463        domain = Domain(points, vertices, boundary)
4464        domain.smooth = False
4465        domain.default_order = 2
4466        domain.beta_w      = 0.9
4467        domain.beta_w_dry  = 0.9
4468        domain.beta_uh     = 0.9
4469        domain.beta_uh_dry = 0.9
4470        domain.beta_vh     = 0.9
4471        domain.beta_vh_dry = 0.9
4472        domain.beta_h = 0.0 #Use first order in h-limiter
4473        domain.H0 = 0.001
4474        domain.tight_slope_limiters = 1
4475
4476        #Bed-slope and friction at vertices (and interpolated elsewhere)
4477        def x_slope(x, y):
4478            return -x/3
4479
4480        domain.set_quantity('elevation', x_slope)
4481
4482        # Boundary conditions
4483        Br = Reflective_boundary(domain)
4484        domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
4485
4486        #Initial condition
4487        domain.set_quantity('stage', expression = 'elevation + 0.05')
4488        domain.check_integrity()
4489
4490        assert allclose(domain.quantities['stage'].centroid_values,
4491                        [0.01296296, 0.03148148, 0.01296296,
4492                        0.03148148, 0.01296296, 0.03148148,
4493                        0.01296296, 0.03148148, 0.01296296,
4494                        0.03148148, 0.01296296, 0.03148148,
4495                        -0.04259259, -0.02407407, -0.04259259,
4496                        -0.02407407, -0.04259259, -0.02407407,
4497                        -0.04259259, -0.02407407, -0.04259259,
4498                        -0.02407407, -0.04259259, -0.02407407,
4499                        -0.09814815, -0.07962963, -0.09814815,
4500                        -0.07962963, -0.09814815, -0.07962963,
4501                        -0.09814815, -0.07962963, -0.09814815,
4502                        -0.07962963, -0.09814815, -0.07962963,
4503                        -0.1537037 , -0.13518519, -0.1537037,
4504                        -0.13518519, -0.1537037, -0.13518519,
4505                        -0.1537037 , -0.13518519, -0.1537037,
4506                        -0.13518519, -0.1537037, -0.13518519,
4507                        -0.20925926, -0.19074074, -0.20925926,
4508                        -0.19074074, -0.20925926, -0.19074074,
4509                        -0.20925926, -0.19074074, -0.20925926,
4510                        -0.19074074, -0.20925926, -0.19074074,
4511                        -0.26481481, -0.2462963, -0.26481481,
4512                        -0.2462963, -0.26481481, -0.2462963,
4513                        -0.26481481, -0.2462963, -0.26481481,
4514                        -0.2462963, -0.26481481, -0.2462963])
4515
4516
4517        #print domain.quantities['stage'].extrapolate_second_order()
4518        #domain.distribute_to_vertices_and_edges()
4519        #print domain.quantities['stage'].vertex_values[:,0]
4520
4521        #Evolution
4522        for t in domain.evolve(yieldstep = 0.05, finaltime = 0.5):
4523            pass
4524
4525
4526        #print domain.quantities['stage'].centroid_values
4527           
4528        assert allclose(domain.quantities['stage'].centroid_values,
4529         [-0.03348416, -0.01749303, -0.03299091, -0.01739241, -0.03246447, -0.01732016,
4530          -0.03205390, -0.01717833, -0.03146383, -0.01699831, -0.03076577, -0.01671795,
4531          -0.07952656, -0.06684763, -0.07721455, -0.06668388, -0.07632976, -0.06600113,
4532          -0.07523678, -0.06546373, -0.07447040, -0.06508861, -0.07438723, -0.06359288,
4533          -0.12526729, -0.11205668, -0.12179433, -0.11068104, -0.12048395, -0.10968948,
4534          -0.11912023, -0.10862628, -0.11784090, -0.10803744, -0.11790629, -0.10742354,
4535          -0.16859613, -0.15427413, -0.16664444, -0.15464452, -0.16570816, -0.15327556,
4536          -0.16409162, -0.15204092, -0.16264608, -0.15102139, -0.16162736, -0.14969205,
4537          -0.18736511, -0.19874036, -0.18811230, -0.19758289, -0.18590182, -0.19580301,
4538          -0.18234588, -0.19423215, -0.18100376, -0.19380116, -0.18509710, -0.19501636,
4539          -0.13982382, -0.14166819, -0.14132775, -0.14528694, -0.14096905, -0.14351126,
4540          -0.13800356, -0.14027920, -0.13613538, -0.13936795, -0.13621902, -0.14204982])
4541
4542                     
4543        assert allclose(domain.quantities['xmomentum'].centroid_values,
4544      [0.00600290,  0.00175780,  0.00591905,  0.00190903,  0.00644462,  0.00203095,
4545       0.00684561,  0.00225089,  0.00708208,  0.00236235,  0.00649095,  0.00222343,
4546       0.02068693,  0.01164034,  0.01983343,  0.01159526,  0.02044611,  0.01233252,
4547       0.02135685,  0.01301289,  0.02161290,  0.01260280,  0.01867612,  0.01133078,
4548       0.04091313,  0.02668283,  0.03634781,  0.02733469,  0.03767692,  0.02836840,
4549       0.03906338,  0.02958073,  0.04025669,  0.02953292,  0.03665616,  0.02583565,
4550       0.06314558,  0.04830935,  0.05663609,  0.04564362,  0.05756200,  0.04739673,
4551       0.05967379,  0.04919083,  0.06124330,  0.04965808,  0.05879240,  0.04629319,
4552       0.08220739,  0.06924725,  0.07713556,  0.06782640,  0.07909499,  0.06992544,
4553       0.08116621,  0.07210181,  0.08281548,  0.07222669,  0.07941059,  0.06755612,
4554       0.01581588,  0.04533609,  0.02017939,  0.04342565,  0.02073232,  0.04476108,
4555       0.02117439,  0.04573358,  0.02129473,  0.04694267,  0.02220398,  0.05533458])
4556
4557       
4558        assert allclose(domain.quantities['ymomentum'].centroid_values,
4559     [-7.65882069e-005, -1.46087080e-004, -1.09630102e-004, -7.80950424e-005,
4560      -1.15922807e-005, -9.09134899e-005, -1.35994542e-004, -1.95673476e-004,
4561      -4.25779199e-004, -2.95890312e-004, -4.00060341e-004, -9.42021290e-005,
4562      -3.41372596e-004, -1.54560195e-004, -2.94810038e-004, -1.08844546e-004,
4563      -6.97240892e-005,  3.50299623e-005, -2.40159184e-004, -2.01805883e-004,
4564      -7.60732405e-004, -5.10897642e-004, -1.00940001e-003, -1.38037759e-004,
4565      -1.06169131e-003, -3.12307760e-004, -9.90602307e-004, -4.21634250e-005,
4566      -6.02424239e-004,  1.52230578e-004, -7.63833035e-004, -1.10273481e-004,
4567      -1.40187071e-003, -5.57831837e-004, -1.63988285e-003, -2.48018092e-004,
4568      -1.83309840e-003, -6.19360836e-004, -1.29955242e-003, -3.76237145e-004,
4569      -1.00613007e-003, -8.63641918e-005, -1.13604124e-003, -3.90589728e-004,
4570      -1.91457355e-003, -9.43783961e-004, -2.28090840e-003, -5.79107025e-004,
4571      -1.54091533e-003, -2.39785792e-003, -2.47947427e-003, -2.02694009e-003,
4572      -2.10441194e-003, -1.82082650e-003, -1.80229336e-003, -2.10418336e-003,
4573      -1.93104408e-003, -2.23200334e-003, -1.57239706e-003, -1.31486358e-003,
4574      -1.17564993e-003, -2.85846494e-003, -3.52956754e-003, -5.12658193e-003,
4575      -6.24238960e-003, -6.01820113e-003, -6.09602201e-003, -5.04787190e-003,
4576      -4.59373845e-003, -3.01393146e-003,  5.08550095e-004, -4.35896549e-004])
4577
4578        os.remove(domain.get_name() + '.sww')
4579
4580
4581    def test_temp_play(self):
4582
4583        from mesh_factory import rectangular
4584        from Numeric import array
4585
4586        #Create basic mesh
4587        points, vertices, boundary = rectangular(5, 5)
4588
4589        #Create shallow water domain
4590        domain = Domain(points, vertices, boundary)
4591        domain.smooth = False
4592        domain.default_order=2
4593        domain.beta_w      = 0.9
4594        domain.beta_w_dry  = 0.9
4595        domain.beta_uh     = 0.9
4596        domain.beta_uh_dry = 0.9
4597        domain.beta_vh     = 0.9
4598        domain.beta_vh_dry = 0.9
4599        domain.beta_h = 0.0 #Use first order in h-limiter
4600       
4601        # FIXME (Ole): Need tests where these two are commented out
4602        domain.H0 = 0        # Backwards compatibility (6/2/7)       
4603        domain.tight_slope_limiters = False # Backwards compatibility (14/4/7)
4604        domain.use_centroid_velocities = False # Backwards compatibility (7/5/8)
4605        domain.use_edge_limiter = False # Backwards compatibility (9/5/8)       
4606       
4607
4608        #Bed-slope and friction at vertices (and interpolated elsewhere)
4609        def x_slope(x, y):
4610            return -x/3
4611
4612        domain.set_quantity('elevation', x_slope)
4613
4614        # Boundary conditions
4615        Br = Reflective_boundary(domain)
4616        domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
4617
4618        #Initial condition
4619        domain.set_quantity('stage', expression='elevation+0.05')
4620        domain.check_integrity()
4621
4622        #Evolution
4623        for t in domain.evolve(yieldstep = 0.05, finaltime = 0.1):
4624            pass
4625
4626        assert allclose(domain.quantities['stage'].centroid_values[:4],
4627                        [0.00206836, 0.01296714, 0.00363415, 0.01438924])
4628        #print domain.quantities['xmomentum'].centroid_values[:4]
4629        assert allclose(domain.quantities['xmomentum'].centroid_values[:4],
4630                        [0.01360154, 0.00671133, 0.01264578, 0.00648503])
4631        assert allclose(domain.quantities['ymomentum'].centroid_values[:4],
4632                        [-1.19201077e-003, -7.23647546e-004,
4633                        -6.39083123e-005, 6.29815168e-005])
4634
4635        os.remove(domain.get_name() + '.sww')
4636
4637    def test_complex_bed(self):
4638        #No friction is tested here
4639
4640        from mesh_factory import rectangular
4641        from Numeric import array
4642
4643        N = 12
4644        points, vertices, boundary = rectangular(N, N/2, len1=1.2,len2=0.6,
4645                                                 origin=(-0.07, 0))
4646
4647
4648        domain = Domain(points, vertices, boundary)
4649        domain.smooth = False
4650        domain.default_order=2
4651
4652
4653        inflow_stage = 0.1
4654        Z = Weir(inflow_stage)
4655        domain.set_quantity('elevation', Z)
4656
4657        Br = Reflective_boundary(domain)
4658        Bd = Dirichlet_boundary([inflow_stage, 0.0, 0.0])
4659        domain.set_boundary({'left': Bd, 'right': Br, 'bottom': Br, 'top': Br})
4660
4661        domain.set_quantity('stage', expression='elevation')
4662
4663        for t in domain.evolve(yieldstep = 0.02, finaltime = 0.2):
4664            pass
4665
4666
4667        #print domain.quantities['stage'].centroid_values
4668
4669        #FIXME: These numbers were from version before 25/10
4670        #assert allclose(domain.quantities['stage'].centroid_values,
4671# [3.95822638e-002,  5.61022588e-002,  4.66437868e-002,  5.73081011e-002,
4672#  4.72394613e-002,  5.74684939e-002,  4.74309483e-002,  5.77458084e-002,
4673#  4.80628177e-002,  5.85656225e-002,  4.90498542e-002,  6.02609831e-002,
4674#  1.18470315e-002,  1.75136443e-002,  1.18035266e-002,  2.15565695e-002,
4675#  1.31620268e-002,  2.14351640e-002,  1.32351076e-002,  2.15450687e-002,
4676#  1.36414028e-002,  2.24274619e-002,  1.51689511e-002,  2.21789655e-002,
4677# -7.54337535e-003, -6.86362021e-004, -7.74146760e-003, -1.83756530e-003,
4678# -8.16773628e-003, -4.49916813e-004, -8.08202599e-003, -3.91118720e-004,
4679# -8.10292716e-003, -3.88584984e-004, -7.35226124e-003,  2.73985295e-004,
4680#  1.86166683e-001,  8.74070369e-002,  1.86166712e-001,  8.74035875e-002,
4681#  6.11666935e-002, -3.76173225e-002, -6.38333276e-002, -3.76147365e-002,
4682#  6.11666725e-002,  8.73846774e-002,  1.86166697e-001,  8.74171550e-002,
4683# -4.83333333e-002,  1.18333333e-001, -4.83333333e-002,  1.18333333e-001,
4684# -4.83333333e-002, -6.66666667e-003, -1.73333333e-001, -1.31666667e-001,
4685# -1.73333333e-001, -6.66666667e-003, -4.83333333e-002,  1.18333333e-001,
4686# -2.48333333e-001, -2.31666667e-001, -2.48333333e-001, -2.31666667e-001,
4687# -2.48333333e-001, -2.31666667e-001, -2.48333333e-001, -2.31666667e-001,
4688# -2.48333333e-001, -2.31666667e-001, -2.48333333e-001, -2.31666667e-001,
4689# -4.65000000e-001, -3.65000000e-001, -4.65000000e-001, -3.65000000e-001,
4690# -4.65000000e-001, -3.65000000e-001, -4.65000000e-001, -3.65000000e-001,
4691# -4.65000000e-001, -3.65000000e-001, -4.65000000e-001, -3.65000000e-001,
4692# -5.98333333e-001, -5.81666667e-001, -5.98333333e-001, -5.81666667e-001,
4693# -5.98333333e-001, -5.81666667e-001, -5.98333333e-001, -5.81666667e-001,
4694# -5.98333333e-001, -5.81666667e-001, -5.98333333e-001, -5.81666667e-001,
4695# -6.48333333e-001, -6.31666667e-001, -6.48333333e-001, -6.31666667e-001,
4696# -6.48333333e-001, -6.31666667e-001, -6.48333333e-001, -6.31666667e-001,
4697# -6.48333333e-001, -6.31666667e-001, -6.48333333e-001, -6.31666667e-001,
4698# -5.31666667e-001, -5.98333333e-001, -5.31666667e-001, -5.98333333e-001,
4699# -5.31666667e-001, -5.98333333e-001, -5.31666667e-001, -5.98333333e-001,
4700# -5.31666667e-001, -5.98333333e-001, -5.31666667e-001, -5.98333333e-001,
4701# -4.98333333e-001, -4.81666667e-001, -4.98333333e-001, -4.81666667e-001,
4702# -4.98333333e-001, -4.81666667e-001, -4.98333333e-001, -4.81666667e-001,
4703# -4.98333333e-001, -4.81666667e-001, -4.98333333e-001, -4.81666667e-001,
4704# -5.48333333e-001, -5.31666667e-001, -5.48333333e-001, -5.31666667e-001,
4705# -5.48333333e-001, -5.31666667e-001, -5.48333333e-001, -5.31666667e-001,
4706# -5.48333333e-001, -5.31666667e-001, -5.48333333e-001, -5.31666667e-001])
4707
4708        os.remove(domain.get_name() + '.sww')
4709
4710    def test_spatio_temporal_boundary_1(self):
4711        """Test that boundary values can be read from file and interpolated
4712        in both time and space.
4713
4714        Verify that the same steady state solution is arrived at and that
4715        time interpolation works.
4716
4717        The full solution history is not exactly the same as
4718        file boundary must read and interpolate from *smoothed* version
4719        as stored in sww.
4720        """
4721        import time
4722
4723        #Create sww file of simple propagation from left to right
4724        #through rectangular domain
4725
4726        from mesh_factory import rectangular
4727
4728        #Create basic mesh
4729        points, vertices, boundary = rectangular(3, 3)
4730
4731        #Create shallow water domain
4732        domain1 = Domain(points, vertices, boundary)
4733
4734        domain1.reduction = mean
4735        domain1.smooth = False  #Exact result
4736
4737        domain1.default_order = 2
4738        domain1.store = True
4739        domain1.set_datadir('.')
4740        domain1.set_name('spatio_temporal_boundary_source' + str(time.time()))
4741
4742        #FIXME: This is extremely important!
4743        #How can we test if they weren't stored?
4744        domain1.quantities_to_be_stored = ['stage', 'xmomentum', 'ymomentum']
4745
4746
4747        #Bed-slope and friction at vertices (and interpolated elsewhere)
4748        domain1.set_quantity('elevation', 0)
4749        domain1.set_quantity('friction', 0)
4750
4751        # Boundary conditions
4752        Br = Reflective_boundary(domain1)
4753        Bd = Dirichlet_boundary([0.3,0,0])
4754        domain1.set_boundary({'left': Bd, 'top': Bd, 'right': Br, 'bottom': Br})
4755        #Initial condition
4756        domain1.set_quantity('stage', 0)
4757        domain1.check_integrity()
4758
4759        finaltime = 5
4760        #Evolution  (full domain - large steps)
4761        for t in domain1.evolve(yieldstep = 0.671, finaltime = finaltime):
4762            pass
4763            #domain1.write_time()
4764
4765        cv1 = domain1.quantities['stage'].centroid_values
4766
4767
4768        #Create a triangle shaped domain (reusing coordinates from domain 1),
4769        #formed from the lower and right hand  boundaries and
4770        #the sw-ne diagonal
4771        #from domain 1. Call it domain2
4772
4773        points = [ [0,0], [1.0/3,0], [1.0/3,1.0/3],
4774                   [2.0/3,0], [2.0/3,1.0/3], [2.0/3,2.0/3],
4775                   [1,0], [1,1.0/3], [1,2.0/3], [1,1]]
4776
4777        vertices = [ [1,2,0], [3,4,1], [2,1,4], [4,5,2],
4778                     [6,7,3], [4,3,7], [7,8,4], [5,4,8], [8,9,5]]
4779
4780        boundary = { (0,1):'bottom', (1,1):'bottom', (4,1): 'bottom',
4781                     (4,2):'right', (6,2):'right', (8,2):'right',
4782                     (0,0):'diagonal', (3,0):'diagonal', (8,0):'diagonal'}
4783
4784        domain2 = Domain(points, vertices, boundary)
4785
4786        domain2.reduction = domain1.reduction
4787        domain2.smooth = False
4788        domain2.default_order = 2
4789
4790        #Bed-slope and friction at vertices (and interpolated elsewhere)
4791        domain2.set_quantity('elevation', 0)
4792        domain2.set_quantity('friction', 0)
4793        domain2.set_quantity('stage', 0)
4794
4795        # Boundary conditions
4796        Br = Reflective_boundary(domain2)
4797        #Bf = Spatio_temporal_boundary(domain1.get_name() + '.' +\
4798        #                              domain1.format, domain2)
4799        Bf = Field_boundary(domain1.get_name() + '.' +\
4800                            domain1.format, domain2)       
4801        domain2.set_boundary({'right':Br, 'bottom':Br, 'diagonal':Bf})
4802        domain2.check_integrity()
4803
4804
4805
4806        #Evolution (small steps)
4807        for t in domain2.evolve(yieldstep = 0.0711, finaltime = finaltime):
4808            pass
4809
4810
4811        #Use output from domain1 as spatio-temporal boundary for domain2
4812        #and verify that results at right hand side are close.
4813
4814        cv2 = domain2.quantities['stage'].centroid_values
4815
4816        #print take(cv1, (12,14,16))  #Right
4817        #print take(cv2, (4,6,8))
4818        #print take(cv1, (0,6,12))  #Bottom
4819        #print take(cv2, (0,1,4))
4820        #print take(cv1, (0,8,16))  #Diag
4821        #print take(cv2, (0,3,8))
4822
4823        assert allclose( take(cv1, (0,8,16)), take(cv2, (0,3,8))) #Diag
4824        assert allclose( take(cv1, (0,6,12)), take(cv2, (0,1,4))) #Bottom
4825        assert allclose( take(cv1, (12,14,16)), take(cv2, (4,6,8))) #RHS
4826
4827        #Cleanup
4828        os.remove(domain1.get_name() + '.' + domain1.format)
4829        os.remove(domain2.get_name() + '.' + domain2.format)       
4830
4831
4832
4833    def test_spatio_temporal_boundary_2(self):
4834        """Test that boundary values can be read from file and interpolated
4835        in both time and space.
4836        This is a more basic test, verifying that boundary object
4837        produces the expected results
4838
4839
4840        """
4841        import time
4842
4843        #Create sww file of simple propagation from left to right
4844        #through rectangular domain
4845
4846        from mesh_factory import rectangular
4847
4848        #Create basic mesh
4849        points, vertices, boundary = rectangular(3, 3)
4850
4851        #Create shallow water domain
4852        domain1 = Domain(points, vertices, boundary)
4853
4854        domain1.reduction = mean
4855        domain1.smooth = True #To mimic MOST output
4856
4857        domain1.default_order = 2
4858        domain1.store = True
4859        domain1.set_datadir('.')
4860        domain1.set_name('spatio_temporal_boundary_source' + str(time.time()))
4861
4862        #FIXME: This is extremely important!
4863        #How can we test if they weren't stored?
4864        domain1.quantities_to_be_stored = ['stage', 'xmomentum', 'ymomentum']
4865
4866
4867        #Bed-slope and friction at vertices (and interpolated elsewhere)
4868        domain1.set_quantity('elevation', 0)
4869        domain1.set_quantity('friction', 0)
4870
4871        # Boundary conditions
4872        Br = Reflective_boundary(domain1)
4873        Bd = Dirichlet_boundary([0.3,0,0])
4874        domain1.set_boundary({'left': Bd, 'top': Bd, 'right': Br, 'bottom': Br})
4875        #Initial condition
4876        domain1.set_quantity('stage', 0)
4877        domain1.check_integrity()
4878
4879        finaltime = 5
4880        #Evolution  (full domain - large steps)
4881        for t in domain1.evolve(yieldstep = 1, finaltime = finaltime):
4882            pass
4883            #domain1.write_time()
4884
4885
4886        #Create an triangle shaped domain (coinciding with some
4887        #coordinates from domain 1),
4888        #formed from the lower and right hand  boundaries and
4889        #the sw-ne diagonal
4890        #from domain 1. Call it domain2
4891
4892        points = [ [0,0], [1.0/3,0], [1.0/3,1.0/3],
4893                   [2.0/3,0], [2.0/3,1.0/3], [2.0/3,2.0/3],
4894                   [1,0], [1,1.0/3], [1,2.0/3], [1,1]]
4895
4896        vertices = [ [1,2,0],
4897                     [3,4,1], [2,1,4], [4,5,2],
4898                     [6,7,3], [4,3,7], [7,8,4], [5,4,8], [8,9,5]]
4899
4900        boundary = { (0,1):'bottom', (1,1):'bottom', (4,1): 'bottom',
4901                     (4,2):'right', (6,2):'right', (8,2):'right',
4902                     (0,0):'diagonal', (3,0):'diagonal', (8,0):'diagonal'}
4903
4904        domain2 = Domain(points, vertices, boundary)
4905
4906        domain2.reduction = domain1.reduction
4907        domain2.smooth = False
4908        domain2.default_order = 2
4909
4910        #Bed-slope and friction at vertices (and interpolated elsewhere)
4911        domain2.set_quantity('elevation', 0)
4912        domain2.set_quantity('friction', 0)
4913        domain2.set_quantity('stage', 0)
4914
4915
4916        #Read results for specific timesteps t=1 and t=2
4917        from Scientific.IO.NetCDF import NetCDFFile
4918        fid = NetCDFFile(domain1.get_name() + '.' + domain1.format)
4919
4920        x = fid.variables['x'][:]
4921        y = fid.variables['y'][:]
4922        s1 = fid.variables['stage'][1,:]
4923        s2 = fid.variables['stage'][2,:]
4924        fid.close()
4925
4926        from Numeric import take, reshape, concatenate
4927        shp = (len(x), 1)
4928        points = concatenate( (reshape(x, shp), reshape(y, shp)), axis=1)
4929        #The diagonal points of domain 1 are 0, 5, 10, 15
4930
4931        #print points[0], points[5], points[10], points[15]
4932        assert allclose( take(points, [0,5,10,15]),
4933                         [[0,0], [1.0/3, 1.0/3], [2.0/3, 2.0/3], [1,1]])
4934
4935
4936        # Boundary conditions
4937        Br = Reflective_boundary(domain2)
4938        #Bf = Spatio_temporal_boundary(domain1.get_name() + '.' + domain1.format,
4939        #                              domain2)
4940        Bf = Field_boundary(domain1.get_name() + '.' + domain1.format,
4941                            domain2, verbose=False)       
4942        domain2.set_boundary({'right':Br, 'bottom':Br, 'diagonal':Bf})
4943        domain2.check_integrity()
4944
4945        #Test that interpolation points are the mid points of the all boundary
4946        #segments
4947
4948        boundary_midpoints = [[1.0/6, 0], [1.0/2, 0], [5.0/6,0],
4949                              [1.0, 1.0/6], [1.0, 1.0/2], [1.0, 5.0/6],
4950                              [1.0/6, 1.0/6], [0.5, 0.5], [5.0/6, 5.0/6]]
4951
4952        boundary_midpoints.sort()
4953        R = Bf.F.interpolation_points.tolist()
4954        R.sort()
4955        assert allclose(boundary_midpoints, R)
4956
4957        #Check spatially interpolated output at time == 1
4958        domain2.time = 1
4959
4960        #First diagonal midpoint
4961        R0 = Bf.evaluate(0,0)
4962        assert allclose(R0[0], (s1[0] + s1[5])/2)
4963
4964        #Second diagonal midpoint
4965        R0 = Bf.evaluate(3,0)
4966        assert allclose(R0[0], (s1[5] + s1[10])/2)
4967
4968        #First diagonal midpoint
4969        R0 = Bf.evaluate(8,0)
4970        assert allclose(R0[0], (s1[10] + s1[15])/2)
4971
4972        #Check spatially interpolated output at time == 2
4973        domain2.time = 2
4974
4975        #First diagonal midpoint
4976        R0 = Bf.evaluate(0,0)
4977        assert allclose(R0[0], (s2[0] + s2[5])/2)
4978
4979        #Second diagonal midpoint
4980        R0 = Bf.evaluate(3,0)
4981        assert allclose(R0[0], (s2[5] + s2[10])/2)
4982
4983        #First diagonal midpoint
4984        R0 = Bf.evaluate(8,0)
4985        assert allclose(R0[0], (s2[10] + s2[15])/2)
4986
4987
4988        #Now check temporal interpolation
4989
4990        domain2.time = 1 + 2.0/3
4991
4992        #First diagonal midpoint
4993        R0 = Bf.evaluate(0,0)
4994        assert allclose(R0[0], ((s1[0] + s1[5])/2 + 2.0*(s2[0] + s2[5])/2)/3)
4995
4996        #Second diagonal midpoint
4997        R0 = Bf.evaluate(3,0)
4998        assert allclose(R0[0], ((s1[5] + s1[10])/2 + 2.0*(s2[5] + s2[10])/2)/3)
4999
5000        #First diagonal midpoint
5001        R0 = Bf.evaluate(8,0)
5002        assert allclose(R0[0], ((s1[10] + s1[15])/2 + 2.0*(s2[10] + s2[15])/2)/3)
5003
5004
5005
5006        #Cleanup
5007        os.remove(domain1.get_name() + '.' + domain1.format)
5008
5009
5010    def test_spatio_temporal_boundary_3(self):
5011        """Test that boundary values can be read from file and interpolated
5012        in both time and space.
5013        This is a more basic test, verifying that boundary object
5014        produces the expected results
5015
5016        This tests adjusting using mean_stage
5017
5018        """
5019
5020        import time
5021
5022        mean_stage = 5.2 # Adjust stage by this amount in boundary
5023
5024        #Create sww file of simple propagation from left to right
5025        #through rectangular domain
5026
5027        from mesh_factory import rectangular
5028
5029        #Create basic mesh
5030        points, vertices, boundary = rectangular(3, 3)
5031
5032        #Create shallow water domain
5033        domain1 = Domain(points, vertices, boundary)
5034
5035        domain1.reduction = mean
5036        domain1.smooth = True #To mimic MOST output
5037
5038        domain1.default_order = 2
5039        domain1.store = True
5040        domain1.set_datadir('.')
5041        domain1.set_name('spatio_temporal_boundary_source' + str(time.time()))
5042
5043        #FIXME: This is extremely important!
5044        #How can we test if they weren't stored?
5045        domain1.quantities_to_be_stored = ['stage', 'xmomentum', 'ymomentum']
5046
5047
5048        #Bed-slope and friction at vertices (and interpolated elsewhere)
5049        domain1.set_quantity('elevation', 0)
5050        domain1.set_quantity('friction', 0)
5051
5052        # Boundary conditions
5053        Br = Reflective_boundary(domain1)
5054        Bd = Dirichlet_boundary([0.3,0,0])
5055        domain1.set_boundary({'left': Bd, 'top': Bd, 'right': Br, 'bottom': Br})
5056        #Initial condition
5057        domain1.set_quantity('stage', 0)
5058        domain1.check_integrity()
5059
5060        finaltime = 5
5061        #Evolution  (full domain - large steps)
5062        for t in domain1.evolve(yieldstep = 1, finaltime = finaltime):
5063            pass
5064            #domain1.write_time()
5065
5066
5067        #Create an triangle shaped domain (coinciding with some
5068        #coordinates from domain 1),
5069        #formed from the lower and right hand  boundaries and
5070        #the sw-ne diagonal
5071        #from domain 1. Call it domain2
5072
5073        points = [ [0,0],
5074                   [1.0/3,0], [1.0/3,1.0/3],
5075                   [2.0/3,0], [2.0/3,1.0/3], [2.0/3,2.0/3],
5076                   [1,0],     [1,1.0/3],     [1,2.0/3],     [1,1]] 
5077                   
5078        vertices = [ [1,2,0],
5079                     [3,4,1], [2,1,4], [4,5,2],
5080                     [6,7,3], [4,3,7], [7,8,4], [5,4,8], [8,9,5]]
5081
5082        boundary = { (0,1):'bottom', (1,1):'bottom', (4,1): 'bottom',
5083                     (4,2):'right', (6,2):'right', (8,2):'right',
5084                     (0,0):'diagonal', (3,0):'diagonal', (8,0):'diagonal'}
5085
5086        domain2 = Domain(points, vertices, boundary)
5087
5088        domain2.reduction = domain1.reduction
5089        domain2.smooth = False
5090        domain2.default_order = 2
5091
5092        #Bed-slope and friction at vertices (and interpolated elsewhere)
5093        domain2.set_quantity('elevation', 0)
5094        domain2.set_quantity('friction', 0)
5095        domain2.set_quantity('stage', 0)
5096
5097
5098        #Read results for specific timesteps t=1 and t=2
5099        from Scientific.IO.NetCDF import NetCDFFile
5100        fid = NetCDFFile(domain1.get_name() + '