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

Last change on this file since 5304 was 5304, checked in by ole, 16 years ago

Comments about use_edge_limiter

File size: 209.4 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.2 # 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 13/10/2006)
1492
1493        G0 = [-0.20000000000000001, -0.19999681443389281, -0.1986192343695303, -0.19147413648863046, -0.19132688908678019, -0.17642317476621105, -0.167376262630034, -0.16192452887426961, -0.15609171725778803, -0.15127107084302249, -0.14048864340360018, -0.19296484125327093, -0.19997006390580363, -0.19999999999937063, -0.19999999999937063, -0.19999999999938772, -0.19999999999938772, -0.19999999999938772, -0.19999999999938772, -0.19974288463035494, -0.19951636867991712, -0.19966301435195755, -0.19981082259800226, -0.19978575003960128, -0.19992942471933109, -0.19999999931029933, -0.19999999999906989, -0.19999999999906989, -0.19999999999906989, -0.19999999999906989, -0.19999999999906989, -0.19999999999906989, -0.19999999999906989, -0.19999999999906989, -0.19999999999906989, -0.19999999999906989, -0.19999999999906989, -0.19999999999906989, -0.19999999999906989, -0.19999999999906989, -0.19999999999906989, -0.19999999999906989, -0.19999999999906989, -0.19999999999906989, -0.19999999999906989, -0.19999999999906989, -0.19999999999906989, -0.19999999999906989, -0.19999999999906989, -0.19999999999906989, -0.19999999999906989]
1494
1495
1496        G1 = [-0.29999999999999993, -0.29988962537199199, -0.29293904425532025, -0.28329367722887888, -0.25999146407696289, -0.22613875068011896, -0.21190705052094994, -0.19900707995208217, -0.18876305176191882, -0.18132447501091936, -0.17395459512711151, -0.15562414200985644, -0.16212999953643359, -0.18964422820514618, -0.20871181844346975, -0.21672207791083464, -0.21774940291862779, -0.21482868050219833, -0.21057786776704043, -0.20649663432591045, -0.20294932949211578, -0.19974459897911329, -0.19733648772704043, -0.19641404599824669, -0.19654095699184146, -0.19709942852191994, -0.19780873983410741, -0.19853259125123518, -0.19916495938961168, -0.19965391267799168, -0.19993539587158982, -0.2001383705551133, -0.20029344332295113, -0.20035349748150011, -0.20029886541561631, -0.20015541958920294, -0.19997273066429103, -0.19979879448668514, -0.19966016997024041, -0.19957558009501869, -0.19955725674938532, -0.19958083002853366, -0.19961752462568647, -0.19965296611330258, -0.19968998132634594, -0.19972532942208607, -0.19975372922008239, -0.19977196116929855, -0.19977951443660594, -0.19977792107284789, -0.19976991595502003]
1497
1498        G2 = [-0.40000000000000002, -0.39011996186687281, -0.33359026016903887, -0.29757449757405952, -0.27594124995715791, -0.25970211955309436, -0.24482929492054245, -0.23156757139219822, -0.21956485769139392, -0.20844522129026694, -0.19856327660654355, -0.18962303467030903, -0.17371085465024955, -0.16429840256208336, -0.17793711732368575, -0.19287799702389993, -0.20236271260796762, -0.20700727993623128, -0.20847704371373174, -0.20796895600687262, -0.20653398626186478, -0.20480656169870676, -0.20295863990994492, -0.20100199602968896, -0.19940642689498472, -0.19858371478015749, -0.19838672154605322, -0.19851093923669558, -0.19878191998909323, -0.19910827645394291, -0.19943514333832094, -0.19971231361970535, -0.19992429278849655, -0.20010744405928019, -0.20025927002359642, -0.20034751667523681, -0.20035504591467249, -0.20029401385620157, -0.20019492358237226, -0.20008934249434918, -0.19999808924091636, -0.19993869218976712, -0.19991589568150098, -0.19991815777945968, -0.19993012995477188, -0.19994576118144997, -0.19996497193815974, -0.19998586151236197, -0.20000487253824847, -0.20001903000364174, -0.20002698661385457]
1499
1500        G3 = [-0.45000000000000001, -0.37713945714588398, -0.33029565026933816, -0.30598209033945367, -0.28847101155177313, -0.27211191064563195, -0.25701544058818926, -0.24298945948410997, -0.23010402733784807, -0.21820351802867713, -0.20709938367218383, -0.19719881806182216, -0.18568281604361933, -0.16828653906676322, -0.16977310768235579, -0.1832707289594605, -0.19483524345250974, -0.20233480051649216, -0.20630757214159207, -0.20763927857964531, -0.20724458160595791, -0.20599191745446047, -0.20438329669495012, -0.20256105512496606, -0.20071269486729407, -0.19934403619901719, -0.19866860191898347, -0.19849975056296071, -0.19860870923007437, -0.19885838217851401, -0.19916422433758982, -0.19946861981642039, -0.19972267778871666, -0.19993013816258154, -0.20011063428833351, -0.20024891930311628, -0.20031882555219671, -0.20031326268593497, -0.20024881068472311, -0.20015443214902759, -0.20005669097631221, -0.19997542564643309, -0.19992564006223304, -0.19990746148869892, -0.19990923999172872, -0.19991956416813192, -0.19993484556273733, -0.1999538628054662, -0.19997381636620407, -0.19999130900268777, -0.20000388227457688]
1501       
1502        #FIXME (DSG):This is a hack so the anuga install, not precompiled
1503        # works on DSG's win2000, python 2.3
1504        #The problem is the gauge_values[X] are 52 long, not 51.
1505        #
1506        # This was probably fixed by Stephen in changeset:3804
1507        if len(gauge_values[0]) == 52: gauge_values[0].pop()
1508        if len(gauge_values[1]) == 52: gauge_values[1].pop()
1509        if len(gauge_values[2]) == 52: gauge_values[2].pop()
1510        if len(gauge_values[3]) == 52: gauge_values[3].pop()
1511
1512##         print len(G0), len(gauge_values[0])
1513##         print len(G1), len(gauge_values[1])
1514       
1515        #print array(gauge_values[0])-array(G0)
1516
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        domain.forcing_terms.append( Rainfall(domain, rate=2.0, polygon = [[1,1], [2,1], [2,2], [1,2]]))
2234
2235        domain.compute_forcing_terms()
2236        #print domain.quantities['stage'].explicit_update
2237       
2238        assert allclose(domain.quantities['stage'].explicit_update[1], 2.0/1000)
2239        assert allclose(domain.quantities['stage'].explicit_update[0], 0)
2240        assert allclose(domain.quantities['stage'].explicit_update[2:], 0)       
2241       
2242        # FIXME: Do Time dependent rainfall with poly
2243
2244
2245
2246
2247    def test_inflow_using_circle(self):
2248        from math import pi, cos, sin
2249
2250        a = [0.0, 0.0]
2251        b = [0.0, 2.0]
2252        c = [2.0, 0.0]
2253        d = [0.0, 4.0]
2254        e = [2.0, 2.0]
2255        f = [4.0, 0.0]
2256
2257        points = [a, b, c, d, e, f]
2258        #bac, bce, ecf, dbe
2259        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
2260
2261
2262        domain = Domain(points, vertices)
2263
2264        # Flat surface with 1m of water
2265        domain.set_quantity('elevation', 0)
2266        domain.set_quantity('stage', 1.0)
2267        domain.set_quantity('friction', 0)
2268
2269        Br = Reflective_boundary(domain)
2270        domain.set_boundary({'exterior': Br})
2271
2272        # Setup only one forcing term, constant inflow of 2 m^3/s on a circle affecting triangles #0 and #1 (bac and bce)
2273        domain.forcing_terms = []
2274        domain.forcing_terms.append( Inflow(domain, rate=2.0, center=(1,1), radius=1) )
2275
2276        domain.compute_forcing_terms()
2277        #print domain.quantities['stage'].explicit_update
2278       
2279        assert allclose(domain.quantities['stage'].explicit_update[1], 2.0/pi)
2280        assert allclose(domain.quantities['stage'].explicit_update[0], 2.0/pi)
2281        assert allclose(domain.quantities['stage'].explicit_update[2:], 0)       
2282
2283
2284    def test_inflow_using_circle_function(self):
2285        from math import pi, cos, sin
2286
2287        a = [0.0, 0.0]
2288        b = [0.0, 2.0]
2289        c = [2.0, 0.0]
2290        d = [0.0, 4.0]
2291        e = [2.0, 2.0]
2292        f = [4.0, 0.0]
2293
2294        points = [a, b, c, d, e, f]
2295        #bac, bce, ecf, dbe
2296        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
2297
2298
2299        domain = Domain(points, vertices)
2300
2301        # Flat surface with 1m of water
2302        domain.set_quantity('elevation', 0)
2303        domain.set_quantity('stage', 1.0)
2304        domain.set_quantity('friction', 0)
2305
2306        Br = Reflective_boundary(domain)
2307        domain.set_boundary({'exterior': Br})
2308
2309        # Setup only one forcing term, time dependent inflow of 2 m^3/s on a circle affecting triangles #0 and #1 (bac and bce)
2310        domain.forcing_terms = []
2311        domain.forcing_terms.append( Inflow(domain, rate=lambda t: 2., center=(1,1), radius=1) )
2312
2313        domain.compute_forcing_terms()
2314       
2315        assert allclose(domain.quantities['stage'].explicit_update[1], 2.0/pi)
2316        assert allclose(domain.quantities['stage'].explicit_update[0], 2.0/pi)
2317        assert allclose(domain.quantities['stage'].explicit_update[2:], 0)       
2318       
2319
2320
2321    #####################################################
2322    def test_first_order_extrapolator_const_z(self):
2323
2324        a = [0.0, 0.0]
2325        b = [0.0, 2.0]
2326        c = [2.0, 0.0]
2327        d = [0.0, 4.0]
2328        e = [2.0, 2.0]
2329        f = [4.0, 0.0]
2330
2331        points = [a, b, c, d, e, f]
2332        #bac, bce, ecf, dbe
2333        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
2334
2335        domain = Domain(points, vertices)
2336        val0 = 2.+2.0/3
2337        val1 = 4.+4.0/3
2338        val2 = 8.+2.0/3
2339        val3 = 2.+8.0/3
2340
2341        zl=zr=-3.75 #Assume constant bed (must be less than stage)
2342        domain.set_quantity('elevation', zl*ones( (4,3) ))
2343        domain.set_quantity('stage', [[val0, val0-1, val0-2],
2344                                      [val1, val1+1, val1],
2345                                      [val2, val2-2, val2],
2346                                      [val3-0.5, val3, val3]])
2347
2348
2349
2350        domain._order_ = 1
2351        domain.distribute_to_vertices_and_edges()
2352
2353        #Check that centroid values were distributed to vertices
2354        C = domain.quantities['stage'].centroid_values
2355        for i in range(3):
2356            assert allclose( domain.quantities['stage'].vertex_values[:,i], C)
2357
2358
2359    def test_first_order_limiter_variable_z(self):
2360        #Check that first order limiter follows bed_slope
2361        from Numeric import alltrue, greater_equal
2362        from anuga.config import epsilon
2363
2364        a = [0.0, 0.0]
2365        b = [0.0, 2.0]
2366        c = [2.0,0.0]
2367        d = [0.0, 4.0]
2368        e = [2.0, 2.0]
2369        f = [4.0,0.0]
2370
2371        points = [a, b, c, d, e, f]
2372        #bac, bce, ecf, dbe
2373        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
2374
2375        domain = Domain(points, vertices)
2376        val0 = 2.+2.0/3
2377        val1 = 4.+4.0/3
2378        val2 = 8.+2.0/3
2379        val3 = 2.+8.0/3
2380
2381        domain.set_quantity('elevation', [[0,0,0], [6,0,0],
2382                                          [6,6,6], [6,6,6]])
2383        domain.set_quantity('stage', [[val0, val0, val0],
2384                                      [val1, val1, val1],
2385                                      [val2, val2, val2],
2386                                      [val3, val3, val3]])
2387
2388        E = domain.quantities['elevation'].vertex_values
2389        L = domain.quantities['stage'].vertex_values
2390
2391
2392        #Check that some stages are not above elevation (within eps)
2393        #- so that the limiter has something to work with
2394        assert not alltrue(alltrue(greater_equal(L,E-epsilon)))
2395
2396        domain._order_ = 1
2397        domain.distribute_to_vertices_and_edges()
2398
2399        #Check that all stages are above elevation (within eps)
2400        assert alltrue(alltrue(greater_equal(L,E-epsilon)))
2401
2402
2403    #####################################################
2404    def test_distribute_basic(self):
2405        #Using test data generated by abstract_2d_finite_volumes-2
2406        #Assuming no friction and flat bed (0.0)
2407
2408        a = [0.0, 0.0]
2409        b = [0.0, 2.0]
2410        c = [2.0, 0.0]
2411        d = [0.0, 4.0]
2412        e = [2.0, 2.0]
2413        f = [4.0, 0.0]
2414
2415        points = [a, b, c, d, e, f]
2416        #bac, bce, ecf, dbe
2417        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
2418
2419        domain = Domain(points, vertices)
2420
2421        val0 = 2.
2422        val1 = 4.
2423        val2 = 8.
2424        val3 = 2.
2425
2426        domain.set_quantity('stage', [val0, val1, val2, val3],
2427                            location='centroids')
2428        L = domain.quantities['stage'].vertex_values
2429
2430        #First order
2431        domain._order_ = 1
2432        domain.distribute_to_vertices_and_edges()
2433        assert allclose(L[1], val1)
2434
2435        #Second order
2436        domain._order_ = 2
2437        domain.beta_w      = 0.9
2438        domain.beta_w_dry  = 0.9
2439        domain.beta_uh     = 0.9
2440        domain.beta_uh_dry = 0.9
2441        domain.beta_vh     = 0.9
2442        domain.beta_vh_dry = 0.9
2443        domain.distribute_to_vertices_and_edges()
2444        assert allclose(L[1], [2.2, 4.9, 4.9])
2445
2446
2447
2448    def test_distribute_away_from_bed(self):
2449        #Using test data generated by abstract_2d_finite_volumes-2
2450        #Assuming no friction and flat bed (0.0)
2451
2452        a = [0.0, 0.0]
2453        b = [0.0, 2.0]
2454        c = [2.0, 0.0]
2455        d = [0.0, 4.0]
2456        e = [2.0, 2.0]
2457        f = [4.0, 0.0]
2458
2459        points = [a, b, c, d, e, f]
2460        #bac, bce, ecf, dbe
2461        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
2462
2463        domain = Domain(points, vertices)
2464        L = domain.quantities['stage'].vertex_values
2465
2466        def stage(x,y):
2467            return x**2
2468
2469        domain.set_quantity('stage', stage, location='centroids')
2470
2471        domain.quantities['stage'].compute_gradients()
2472
2473        a, b = domain.quantities['stage'].get_gradients()
2474               
2475        assert allclose(a[1], 3.33333334)
2476        assert allclose(b[1], 0.0)
2477
2478        domain._order_ = 1
2479        domain.distribute_to_vertices_and_edges()
2480        assert allclose(L[1], 1.77777778)
2481
2482        domain._order_ = 2
2483        domain.beta_w      = 0.9
2484        domain.beta_w_dry  = 0.9
2485        domain.beta_uh     = 0.9
2486        domain.beta_uh_dry = 0.9
2487        domain.beta_vh     = 0.9
2488        domain.beta_vh_dry = 0.9
2489        domain.distribute_to_vertices_and_edges()
2490        assert allclose(L[1], [0.57777777, 2.37777778, 2.37777778])
2491
2492
2493
2494    def test_distribute_away_from_bed1(self):
2495        #Using test data generated by abstract_2d_finite_volumes-2
2496        #Assuming no friction and flat bed (0.0)
2497
2498        a = [0.0, 0.0]
2499        b = [0.0, 2.0]
2500        c = [2.0, 0.0]
2501        d = [0.0, 4.0]
2502        e = [2.0, 2.0]
2503        f = [4.0, 0.0]
2504
2505        points = [a, b, c, d, e, f]
2506        #bac, bce, ecf, dbe
2507        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
2508
2509        domain = Domain(points, vertices)
2510        L = domain.quantities['stage'].vertex_values
2511
2512        def stage(x,y):
2513            return x**4+y**2
2514
2515        domain.set_quantity('stage', stage, location='centroids')
2516        #print domain.quantities['stage'].centroid_values
2517
2518        domain.quantities['stage'].compute_gradients()
2519        a, b = domain.quantities['stage'].get_gradients()
2520        assert allclose(a[1], 25.18518519)
2521        assert allclose(b[1], 3.33333333)
2522
2523        domain._order_ = 1
2524        domain.distribute_to_vertices_and_edges()
2525        assert allclose(L[1], 4.9382716)
2526
2527        domain._order_ = 2
2528        domain.beta_w      = 0.9
2529        domain.beta_w_dry  = 0.9
2530        domain.beta_uh     = 0.9
2531        domain.beta_uh_dry = 0.9
2532        domain.beta_vh     = 0.9
2533        domain.beta_vh_dry = 0.9
2534        domain.distribute_to_vertices_and_edges()
2535        assert allclose(L[1], [1.07160494, 6.46058131, 7.28262855])
2536
2537
2538
2539    def test_distribute_near_bed(self):
2540
2541        a = [0.0, 0.0]
2542        b = [0.0, 2.0]
2543        c = [2.0, 0.0]
2544        d = [0.0, 4.0]
2545        e = [2.0, 2.0]
2546        f = [4.0, 0.0]
2547
2548        points = [a, b, c, d, e, f]
2549        #bac, bce, ecf, dbe
2550        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
2551
2552        domain = Domain(points, vertices)
2553
2554
2555        #Set up for a gradient of (10,0) at mid triangle (bce)
2556        def slope(x, y):
2557            return 10*x
2558
2559        h = 0.1
2560        def stage(x, y):
2561            return slope(x, y) + h
2562
2563        domain.set_quantity('elevation', slope)
2564        domain.set_quantity('stage', stage, location='centroids')
2565
2566        #print domain.quantities['elevation'].centroid_values
2567        #print domain.quantities['stage'].centroid_values
2568
2569        E = domain.quantities['elevation'].vertex_values
2570        L = domain.quantities['stage'].vertex_values
2571
2572        # Get reference values
2573        volumes = []
2574        for i in range(len(L)):
2575            volumes.append(sum(L[i])/3)
2576            assert allclose(volumes[i], domain.quantities['stage'].centroid_values[i]) 
2577       
2578       
2579        domain._order_ = 1
2580       
2581        domain.tight_slope_limiters = 0
2582        domain.distribute_to_vertices_and_edges()
2583        assert allclose(L[1], [0.1, 20.1, 20.1])
2584        for i in range(len(L)):
2585            assert allclose(volumes[i], sum(L[i])/3)                   
2586       
2587        domain.tight_slope_limiters = 1 # Allow triangle to be flatter (closer to bed)
2588        domain.distribute_to_vertices_and_edges()
2589        assert allclose(L[1], [0.298, 20.001, 20.001])
2590        for i in range(len(L)):
2591            assert allclose(volumes[i], sum(L[i])/3)   
2592
2593        domain._order_ = 2
2594       
2595        domain.tight_slope_limiters = 0
2596        domain.distribute_to_vertices_and_edges()
2597        assert allclose(L[1], [0.1, 20.1, 20.1])       
2598        for i in range(len(L)):
2599            assert allclose(volumes[i], sum(L[i])/3)           
2600       
2601        domain.tight_slope_limiters = 1 # Allow triangle to be flatter (closer to bed)
2602        domain.distribute_to_vertices_and_edges()
2603        assert allclose(L[1], [0.298, 20.001, 20.001])
2604        for i in range(len(L)):
2605            assert allclose(volumes[i], sum(L[i])/3)   
2606       
2607
2608
2609    def test_distribute_near_bed1(self):
2610
2611        a = [0.0, 0.0]
2612        b = [0.0, 2.0]
2613        c = [2.0, 0.0]
2614        d = [0.0, 4.0]
2615        e = [2.0, 2.0]
2616        f = [4.0, 0.0]
2617
2618        points = [a, b, c, d, e, f]
2619        #bac, bce, ecf, dbe
2620        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
2621
2622        domain = Domain(points, vertices)
2623
2624
2625        #Set up for a gradient of (8,2) at mid triangle (bce)
2626        def slope(x, y):
2627            return x**4+y**2
2628
2629        h = 0.1
2630        def stage(x,y):
2631            return slope(x,y)+h
2632
2633        domain.set_quantity('elevation', slope)
2634        domain.set_quantity('stage', stage)
2635
2636        #print domain.quantities['elevation'].centroid_values
2637        #print domain.quantities['stage'].centroid_values
2638
2639        E = domain.quantities['elevation'].vertex_values
2640        L = domain.quantities['stage'].vertex_values
2641
2642        # Get reference values
2643        volumes = []
2644        for i in range(len(L)):
2645            volumes.append(sum(L[i])/3)
2646            assert allclose(volumes[i], domain.quantities['stage'].centroid_values[i]) 
2647       
2648        #print E
2649        domain._order_ = 1
2650       
2651        domain.tight_slope_limiters = 0
2652        domain.distribute_to_vertices_and_edges()
2653        assert allclose(L[1], [4.1, 16.1, 20.1])       
2654        for i in range(len(L)):
2655            assert allclose(volumes[i], sum(L[i])/3)
2656       
2657               
2658        domain.tight_slope_limiters = 1 # Allow triangle to be flatter (closer to bed)
2659        domain.distribute_to_vertices_and_edges()
2660        assert allclose(L[1], [4.2386, 16.0604, 20.001])
2661        for i in range(len(L)):
2662            assert allclose(volumes[i], sum(L[i])/3)   
2663       
2664
2665        domain._order_ = 2
2666       
2667        domain.tight_slope_limiters = 0   
2668        domain.distribute_to_vertices_and_edges()
2669        assert allclose(L[1], [4.1, 16.1, 20.1])
2670        for i in range(len(L)):
2671            assert allclose(volumes[i], sum(L[i])/3)   
2672       
2673        domain.tight_slope_limiters = 1 # Allow triangle to be flatter (closer to bed)
2674        domain.distribute_to_vertices_and_edges()
2675        #print L[1]
2676        assert allclose(L[1], [4.23370103, 16.06529897, 20.001]) or\
2677               allclose(L[1], [4.18944138, 16.10955862, 20.001]) or\
2678               allclose(L[1], [4.19351461, 16.10548539, 20.001]) # old limiters
2679       
2680        for i in range(len(L)):
2681            assert allclose(volumes[i], sum(L[i])/3)
2682
2683
2684    def test_second_order_distribute_real_data(self):
2685        #Using test data generated by abstract_2d_finite_volumes-2
2686        #Assuming no friction and flat bed (0.0)
2687
2688        a = [0.0, 0.0]
2689        b = [0.0, 1.0/5]
2690        c = [0.0, 2.0/5]
2691        d = [1.0/5, 0.0]
2692        e = [1.0/5, 1.0/5]
2693        f = [1.0/5, 2.0/5]
2694        g = [2.0/5, 2.0/5]
2695
2696        points = [a, b, c, d, e, f, g]
2697        #bae, efb, cbf, feg
2698        vertices = [ [1,0,4], [4,5,1], [2,1,5], [5,4,6]]
2699
2700        domain = Domain(points, vertices)
2701
2702        def slope(x, y):
2703            return -x/3
2704
2705        domain.set_quantity('elevation', slope)
2706        domain.set_quantity('stage',
2707                            [0.01298164, 0.00365611,
2708                             0.01440365, -0.0381856437096],
2709                            location='centroids')
2710        domain.set_quantity('xmomentum',
2711                            [0.00670439, 0.01263789,
2712                             0.00647805, 0.0178180740668],
2713                            location='centroids')
2714        domain.set_quantity('ymomentum',
2715                            [-7.23510980e-004, -6.30413883e-005,
2716                             6.30413883e-005, 0.000200907255866],
2717                            location='centroids')
2718
2719        E = domain.quantities['elevation'].vertex_values
2720        L = domain.quantities['stage'].vertex_values
2721        X = domain.quantities['xmomentum'].vertex_values
2722        Y = domain.quantities['ymomentum'].vertex_values
2723
2724        #print E
2725        domain._order_ = 2
2726        domain.beta_w      = 0.9
2727        domain.beta_w_dry  = 0.9
2728        domain.beta_uh     = 0.9
2729        domain.beta_uh_dry = 0.9
2730        domain.beta_vh     = 0.9
2731        domain.beta_vh_dry = 0.9
2732        domain.beta_h = 0.0 #Use first order in h-limiter
2733       
2734        # FIXME (Ole): Need tests where this is commented out
2735        domain.tight_slope_limiters = 0 # Backwards compatibility (14/4/7)                 
2736        domain.use_centroid_velocities = 0 # Backwards compatibility (7/5/8)
2737       
2738               
2739        domain.distribute_to_vertices_and_edges()
2740
2741        #print L[1,:]
2742        #print X[1,:]
2743        #print Y[1,:]
2744
2745        assert allclose(L[1,:], [-0.00825735775384,
2746                                 -0.00801881482869,
2747                                 0.0272445025825])
2748        assert allclose(X[1,:], [0.0143507718962,
2749                                 0.0142502147066,
2750                                 0.00931268339717])
2751        assert allclose(Y[1,:], [-0.000117062180693,
2752                                 7.94434448109e-005,
2753                                 -0.000151505429018])
2754
2755
2756
2757    def test_balance_deep_and_shallow(self):
2758        """Test that balanced limiters preserve conserved quantites.
2759        This test is using old depth based balanced limiters
2760        """
2761        import copy
2762
2763        a = [0.0, 0.0]
2764        b = [0.0, 2.0]
2765        c = [2.0, 0.0]
2766        d = [0.0, 4.0]
2767        e = [2.0, 2.0]
2768        f = [4.0, 0.0]
2769
2770        points = [a, b, c, d, e, f]
2771
2772        #bac, bce, ecf, dbe
2773        elements = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4] ]
2774
2775        domain = Domain(points, elements)
2776        domain.check_integrity()
2777
2778        #Create a deliberate overshoot
2779        domain.set_quantity('stage', [[3,0,3], [2,2,6], [5,3,8], [8,3,5]])
2780        domain.set_quantity('elevation', 0) #Flat bed
2781        stage = domain.quantities['stage']
2782
2783        ref_centroid_values = copy.copy(stage.centroid_values[:]) #Copy
2784
2785        #Limit
2786        domain.tight_slope_limiters = 0               
2787        domain.distribute_to_vertices_and_edges()
2788
2789        #Assert that quantities are conserved
2790        from Numeric import sum
2791        for k in range(len(domain)):
2792            assert allclose (ref_centroid_values[k],
2793                             sum(stage.vertex_values[k,:])/3)
2794
2795
2796        #Now try with a non-flat bed - closely hugging initial stage in places
2797        #This will create alphas in the range [0, 0.478260, 1]
2798        domain.set_quantity('stage', [[3,0,3], [2,2,6], [5,3,8], [8,3,5]])
2799        domain.set_quantity('elevation', [[0,0,0],
2800                                        [1.8,1.9,5.9],
2801                                        [4.6,0,0],
2802                                        [0,2,4]])
2803        stage = domain.quantities['stage']
2804
2805        ref_centroid_values = copy.copy(stage.centroid_values[:]) #Copy
2806        ref_vertex_values = copy.copy(stage.vertex_values[:]) #Copy
2807
2808        #Limit
2809        domain.tight_slope_limiters = 0       
2810        domain.distribute_to_vertices_and_edges()
2811
2812
2813        #Assert that all vertex quantities have changed
2814        for k in range(len(domain)):
2815            #print ref_vertex_values[k,:], stage.vertex_values[k,:]
2816            assert not allclose (ref_vertex_values[k,:], stage.vertex_values[k,:])
2817        #and assert that quantities are still conserved
2818        from Numeric import sum
2819        for k in range(len(domain)):
2820            assert allclose (ref_centroid_values[k],
2821                             sum(stage.vertex_values[k,:])/3)
2822
2823
2824        # Check actual results
2825        assert allclose (stage.vertex_values,
2826                         [[2,2,2],
2827                          [1.93333333, 2.03333333, 6.03333333],
2828                          [6.93333333, 4.53333333, 4.53333333],
2829                          [5.33333333, 5.33333333, 5.33333333]])
2830
2831
2832    def test_balance_deep_and_shallow_tight_SL(self):
2833        """Test that balanced limiters preserve conserved quantites.
2834        This test is using Tight Slope Limiters
2835        """
2836        import copy
2837
2838        a = [0.0, 0.0]
2839        b = [0.0, 2.0]
2840        c = [2.0, 0.0]
2841        d = [0.0, 4.0]
2842        e = [2.0, 2.0]
2843        f = [4.0, 0.0]
2844
2845        points = [a, b, c, d, e, f]
2846
2847        #bac, bce, ecf, dbe
2848        elements = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4] ]
2849
2850        domain = Domain(points, elements)
2851        domain.check_integrity()
2852
2853        #Create a deliberate overshoot
2854        domain.set_quantity('stage', [[3,0,3], [2,2,6], [5,3,8], [8,3,5]])
2855        domain.set_quantity('elevation', 0) #Flat bed
2856        stage = domain.quantities['stage']
2857
2858        ref_centroid_values = copy.copy(stage.centroid_values[:]) #Copy
2859
2860        #Limit
2861        domain.tight_slope_limiters = 1               
2862        domain.distribute_to_vertices_and_edges()
2863
2864        #Assert that quantities are conserved
2865        from Numeric import sum
2866        for k in range(len(domain)):
2867            assert allclose (ref_centroid_values[k],
2868                             sum(stage.vertex_values[k,:])/3)
2869
2870
2871        #Now try with a non-flat bed - closely hugging initial stage in places
2872        #This will create alphas in the range [0, 0.478260, 1]
2873        domain.set_quantity('stage', [[3,0,3], [2,2,6], [5,3,8], [8,3,5]])
2874        domain.set_quantity('elevation', [[0,0,0],
2875                                        [1.8,1.9,5.9],
2876                                        [4.6,0,0],
2877                                        [0,2,4]])
2878        stage = domain.quantities['stage']
2879
2880        ref_centroid_values = copy.copy(stage.centroid_values[:]) #Copy
2881        ref_vertex_values = copy.copy(stage.vertex_values[:]) #Copy
2882
2883        #Limit
2884        domain.tight_slope_limiters = 1       
2885        domain.distribute_to_vertices_and_edges()
2886
2887
2888        #Assert that all vertex quantities have changed
2889        for k in range(len(domain)):
2890            #print ref_vertex_values[k,:], stage.vertex_values[k,:]
2891            assert not allclose (ref_vertex_values[k,:], stage.vertex_values[k,:])
2892        #and assert that quantities are still conserved
2893        from Numeric import sum
2894        for k in range(len(domain)):
2895            assert allclose (ref_centroid_values[k],
2896                             sum(stage.vertex_values[k,:])/3)
2897
2898
2899        #Also check that Python and C version produce the same
2900        # No longer applicable if tight_slope_limiters == 1
2901        #print stage.vertex_values
2902        #assert allclose (stage.vertex_values,
2903        #                 [[2,2,2],
2904        #                  [1.93333333, 2.03333333, 6.03333333],
2905        #                  [6.93333333, 4.53333333, 4.53333333],
2906        #                  [5.33333333, 5.33333333, 5.33333333]])
2907
2908
2909
2910    def test_balance_deep_and_shallow_Froude(self):
2911        """Test that balanced limiters preserve conserved quantites -
2912        and also that excessive Froude numbers are dealt with.
2913        This test is using tight slope limiters.
2914        """
2915        import copy
2916        from Numeric import sqrt, absolute
2917
2918        a = [0.0, 0.0]
2919        b = [0.0, 2.0]
2920        c = [2.0, 0.0]
2921        d = [0.0, 4.0]
2922        e = [2.0, 2.0]
2923        f = [4.0, 0.0]
2924
2925        points = [a, b, c, d, e, f]
2926
2927        # bac, bce, ecf, dbe
2928        elements = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4] ]
2929
2930        domain = Domain(points, elements)
2931        domain.check_integrity()
2932        domain.tight_slope_limiters = True
2933        domain.use_centroid_velocities = True               
2934
2935        # Create non-flat bed - closely hugging initial stage in places
2936        # This will create alphas in the range [0, 0.478260, 1]
2937        domain.set_quantity('stage', [[3,0,3], [2,2,6], [5,3,8], [8,3,5]])
2938        domain.set_quantity('elevation', [[0,0,0],
2939                                        [1.8,1.999,5.999],
2940                                        [4.6,0,0],
2941                                        [0,2,4]])
2942
2943        # Create small momenta, that nonetheless will generate large speeds
2944        # due to shallow depth at isolated vertices
2945        domain.set_quantity('xmomentum', -0.0058)
2946        domain.set_quantity('ymomentum', 0.0890)
2947
2948
2949
2950       
2951        stage = domain.quantities['stage']
2952        elevation = domain.quantities['elevation']
2953        xmomentum = domain.quantities['xmomentum']
2954        ymomentum = domain.quantities['ymomentum']       
2955
2956        # Setup triangle #1 to mimick real Froude explosion observed
2957        # in the Onslow example 13 Nov 2007.
2958
2959        stage.vertex_values[1,:] = [1.6385, 1.6361, 1.2953]
2960        elevation.vertex_values[1,:] = [1.6375, 1.6336, 0.4647]       
2961        xmomentum.vertex_values[1,:] = [-0.0058, -0.0050, -0.0066]
2962        ymomentum.vertex_values[1,:] = [0.0890, 0.0890, 0.0890]
2963
2964        xmomentum.interpolate()
2965        ymomentum.interpolate()       
2966        stage.interpolate()       
2967        elevation.interpolate()
2968
2969        # Verify interpolation
2970        assert allclose(stage.centroid_values[1], 1.5233)
2971        assert allclose(elevation.centroid_values[1], 1.2452667)
2972        assert allclose(xmomentum.centroid_values[1], -0.0058)
2973        assert allclose(ymomentum.centroid_values[1], 0.089)
2974
2975        # Derived quantities
2976        depth = stage-elevation
2977        u = xmomentum/depth
2978        v = ymomentum/depth
2979
2980        denom = (depth*g)**0.5 
2981        Fx = u/denom
2982        Fy = v/denom
2983       
2984   
2985        # Verify against Onslow example (14 Nov 2007)
2986        assert allclose(depth.centroid_values[1], 0.278033)
2987        assert allclose(u.centroid_values[1], -0.0208608)
2988        assert allclose(v.centroid_values[1], 0.3201055)
2989
2990        assert allclose(denom.centroid_values[1],
2991                        sqrt(depth.centroid_values[1]*g))
2992
2993        assert allclose(u.centroid_values[1]/denom.centroid_values[1],
2994                        -0.012637746977)
2995        assert allclose(Fx.centroid_values[1],
2996                        u.centroid_values[1]/denom.centroid_values[1])
2997
2998        # Check that Froude numbers are small at centroids.
2999        assert allclose(Fx.centroid_values[1], -0.012637746977)
3000        assert allclose(Fy.centroid_values[1], 0.193924048435)
3001
3002
3003        # But Froude numbers are huge at some vertices and edges
3004        assert allclose(Fx.vertex_values[1,:], [-5.85888475e+01,
3005                                                -1.27775313e+01,
3006                                                -2.78511420e-03])
3007
3008        assert allclose(Fx.edge_values[1,:], [-6.89150773e-03,
3009                                              -7.38672488e-03,
3010                                              -2.35626238e+01])
3011
3012        assert allclose(Fy.vertex_values[1,:], [8.99035764e+02,
3013                                                2.27440057e+02,
3014                                                3.75568430e-02])
3015
3016        assert allclose(Fy.edge_values[1,:], [1.05748998e-01,
3017                                              1.06035244e-01,
3018                                              3.88346947e+02])
3019       
3020       
3021        # The task is now to arrange the limiters such that Froude numbers
3022        # remain under control whil at the same time obeying the conservation
3023        # laws.
3024
3025       
3026        ref_centroid_values = copy.copy(stage.centroid_values[:]) #Copy
3027        ref_vertex_values = copy.copy(stage.vertex_values[:]) #Copy
3028
3029        # Limit (and invoke balance_deep_and_shallow)
3030        domain.tight_slope_limiters = 1
3031        domain.distribute_to_vertices_and_edges()
3032
3033        # Redo derived quantities
3034        depth = stage-elevation
3035        u = xmomentum/depth
3036        v = ymomentum/depth
3037
3038        # Assert that all vertex velocities stay within one
3039        # order of magnitude of centroid velocities.
3040        #print u.vertex_values[1,:]
3041        #print u.centroid_values[1]
3042       
3043        assert alltrue(absolute(u.vertex_values[1,:]) <= absolute(u.centroid_values[1])*10)
3044        assert alltrue(absolute(v.vertex_values[1,:]) <= absolute(v.centroid_values[1])*10) 
3045
3046        denom = (depth*g)**0.5 
3047        Fx = u/denom
3048        Fy = v/denom
3049
3050
3051        # Assert that Froude numbers are less than max value (TBA)
3052        # at vertices, edges and centroids.
3053        from anuga.config import maximum_froude_number
3054        assert alltrue(absolute(Fx.vertex_values[1,:]) < maximum_froude_number)
3055        assert alltrue(absolute(Fy.vertex_values[1,:]) < maximum_froude_number)
3056
3057
3058        # Assert that all vertex quantities have changed
3059        for k in range(len(domain)):
3060            #print ref_vertex_values[k,:], stage.vertex_values[k,:]
3061            assert not allclose (ref_vertex_values[k,:],
3062                                 stage.vertex_values[k,:])
3063           
3064        # Assert that quantities are still conserved
3065        from Numeric import sum
3066        for k in range(len(domain)):
3067            assert allclose (ref_centroid_values[k],
3068                             sum(stage.vertex_values[k,:])/3)
3069
3070
3071       
3072        return
3073   
3074        qwidth = 12
3075        for k in [1]: #range(len(domain)):
3076            print 'Triangle %d (C, V, E)' %k
3077           
3078            print 'stage'.ljust(qwidth), stage.centroid_values[k],\
3079                  stage.vertex_values[k,:], stage.edge_values[k,:]
3080            print 'elevation'.ljust(qwidth), elevation.centroid_values[k],\
3081                  elevation.vertex_values[k,:], elevation.edge_values[k,:]
3082            print 'depth'.ljust(qwidth), depth.centroid_values[k],\
3083                  depth.vertex_values[k,:], depth.edge_values[k,:]
3084            print 'xmomentum'.ljust(qwidth), xmomentum.centroid_values[k],\
3085                  xmomentum.vertex_values[k,:], xmomentum.edge_values[k,:]
3086            print 'ymomentum'.ljust(qwidth), ymomentum.centroid_values[k],\
3087                  ymomentum.vertex_values[k,:], ymomentum.edge_values[k,:]
3088            print 'u'.ljust(qwidth),u.centroid_values[k],\
3089                  u.vertex_values[k,:], u.edge_values[k,:]
3090            print 'v'.ljust(qwidth), v.centroid_values[k],\
3091                  v.vertex_values[k,:], v.edge_values[k,:]
3092            print 'Fx'.ljust(qwidth), Fx.centroid_values[k],\
3093                  Fx.vertex_values[k,:], Fx.edge_values[k,:]
3094            print 'Fy'.ljust(qwidth), Fy.centroid_values[k],\
3095                  Fy.vertex_values[k,:], Fy.edge_values[k,:]
3096           
3097           
3098       
3099
3100
3101
3102    def test_conservation_1(self):
3103        """Test that stage is conserved globally
3104
3105        This one uses a flat bed, reflective bdries and a suitable
3106        initial condition
3107        """
3108        from mesh_factory import rectangular
3109        from Numeric import array
3110
3111        #Create basic mesh
3112        points, vertices, boundary = rectangular(6, 6)
3113
3114        #Create shallow water domain
3115        domain = Domain(points, vertices, boundary)
3116        domain.smooth = False
3117        domain.default_order=2
3118
3119        #IC
3120        def x_slope(x, y):
3121            return x/3
3122
3123        domain.set_quantity('elevation', 0)
3124        domain.set_quantity('friction', 0)
3125        domain.set_quantity('stage', x_slope)
3126
3127        # Boundary conditions (reflective everywhere)
3128        Br = Reflective_boundary(domain)
3129        domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
3130
3131        domain.check_integrity()
3132
3133        initial_volume = domain.quantities['stage'].get_integral()
3134        initial_xmom = domain.quantities['xmomentum'].get_integral()
3135
3136        #print initial_xmom
3137
3138        #Evolution
3139        for t in domain.evolve(yieldstep = 0.05, finaltime = 5.0):
3140            volume =  domain.quantities['stage'].get_integral()
3141            assert allclose (volume, initial_volume)
3142
3143            #I don't believe that the total momentum should be the same
3144            #It starts with zero and ends with zero though
3145            #xmom = domain.quantities['xmomentum'].get_integral()
3146            #print xmom
3147            #assert allclose (xmom, initial_xmom)
3148
3149        os.remove(domain.get_name() + '.sww')
3150
3151
3152    def test_conservation_2(self):
3153        """Test that stage is conserved globally
3154
3155        This one uses a slopy bed, reflective bdries and a suitable
3156        initial condition
3157        """
3158        from mesh_factory import rectangular
3159        from Numeric import array
3160
3161        #Create basic mesh
3162        points, vertices, boundary = rectangular(6, 6)
3163
3164        #Create shallow water domain
3165        domain = Domain(points, vertices, boundary)
3166        domain.smooth = False
3167        domain.default_order=2
3168
3169        #IC
3170        def x_slope(x, y):
3171            return x/3
3172
3173        domain.set_quantity('elevation', x_slope)
3174        domain.set_quantity('friction', 0)
3175        domain.set_quantity('stage', 0.4) #Steady
3176
3177        # Boundary conditions (reflective everywhere)
3178        Br = Reflective_boundary(domain)
3179        domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
3180
3181        domain.check_integrity()
3182
3183        initial_volume = domain.quantities['stage'].get_integral()
3184        initial_xmom = domain.quantities['xmomentum'].get_integral()
3185
3186        #print initial_xmom
3187
3188        #Evolution
3189        for t in domain.evolve(yieldstep = 0.05, finaltime = 5.0):
3190            volume =  domain.quantities['stage'].get_integral()
3191            assert allclose (volume, initial_volume)
3192
3193            #FIXME: What would we expect from momentum
3194            #xmom = domain.quantities['xmomentum'].get_integral()
3195            #print xmom
3196            #assert allclose (xmom, initial_xmom)
3197
3198        os.remove(domain.get_name() + '.sww')
3199
3200    def test_conservation_3(self):
3201        """Test that stage is conserved globally
3202
3203        This one uses a larger grid, convoluted bed, reflective bdries and a suitable
3204        initial condition
3205        """
3206        from mesh_factory import rectangular
3207        from Numeric import array
3208
3209        #Create basic mesh
3210        points, vertices, boundary = rectangular(2, 1)
3211
3212        #Create shallow water domain
3213        domain = Domain(points, vertices, boundary)
3214        domain.smooth = False
3215        domain.default_order = 2
3216        domain.beta_h = 0.2
3217        domain.set_quantities_to_be_stored(['stage', 'xmomentum', 'ymomentum'])
3218
3219        #IC
3220        def x_slope(x, y):
3221            z = 0*x
3222            for i in range(len(x)):
3223                if x[i] < 0.3:
3224                    z[i] = x[i]/3
3225                if 0.3 <= x[i] < 0.5:
3226                    z[i] = -0.5
3227                if 0.5 <= x[i] < 0.7:
3228                    z[i] = 0.39
3229                if 0.7 <= x[i]:
3230                    z[i] = x[i]/3
3231            return z
3232
3233
3234
3235        domain.set_quantity('elevation', x_slope)
3236        domain.set_quantity('friction', 0)
3237        domain.set_quantity('stage', 0.4) #Steady
3238
3239        # Boundary conditions (reflective everywhere)
3240        Br = Reflective_boundary(domain)
3241        domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
3242
3243        domain.check_integrity()
3244
3245        initial_volume = domain.quantities['stage'].get_integral()
3246        initial_xmom = domain.quantities['xmomentum'].get_integral()
3247
3248        import copy
3249        ref_centroid_values =\
3250                 copy.copy(domain.quantities['stage'].centroid_values)
3251
3252        #print 'ORG', domain.quantities['stage'].centroid_values
3253        domain.distribute_to_vertices_and_edges()
3254
3255
3256        #print domain.quantities['stage'].centroid_values
3257        assert allclose(domain.quantities['stage'].centroid_values,
3258                        ref_centroid_values)
3259
3260
3261        #Check that initial limiter doesn't violate cons quan
3262        assert allclose (domain.quantities['stage'].get_integral(),
3263                         initial_volume)
3264
3265        #Evolution
3266        for t in domain.evolve(yieldstep = 0.05, finaltime = 10):
3267            volume =  domain.quantities['stage'].get_integral()
3268            #print t, volume, initial_volume
3269            assert allclose (volume, initial_volume)
3270
3271        os.remove(domain.get_name() + '.sww')
3272
3273    def test_conservation_4(self):
3274        """Test that stage is conserved globally
3275
3276        This one uses a larger grid, convoluted bed, reflective bdries and a suitable
3277        initial condition
3278        """
3279        from mesh_factory import rectangular
3280        from Numeric import array
3281
3282        #Create basic mesh
3283        points, vertices, boundary = rectangular(6, 6)
3284
3285        #Create shallow water domain
3286        domain = Domain(points, vertices, boundary)
3287        domain.smooth = False
3288        domain.default_order=2
3289        domain.beta_h = 0.0
3290
3291        #IC
3292        def x_slope(x, y):
3293            z = 0*x
3294            for i in range(len(x)):
3295                if x[i] < 0.3:
3296                    z[i] = x[i]/3
3297                if 0.3 <= x[i] < 0.5:
3298                    z[i] = -0.5
3299                if 0.5 <= x[i] < 0.7:
3300                    #z[i] = 0.3 #OK with beta == 0.2
3301                    z[i] = 0.34 #OK with beta == 0.0
3302                    #z[i] = 0.35#Fails after 80 timesteps with an error
3303                                #of the order 1.0e-5
3304                if 0.7 <= x[i]:
3305                    z[i] = x[i]/3
3306            return z
3307
3308
3309
3310        domain.set_quantity('elevation', x_slope)
3311        domain.set_quantity('friction', 0)
3312        domain.set_quantity('stage', 0.4) #Steady
3313
3314        # Boundary conditions (reflective everywhere)
3315        Br = Reflective_boundary(domain)
3316        domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
3317
3318        domain.check_integrity()
3319
3320        initial_volume = domain.quantities['stage'].get_integral()
3321        initial_xmom = domain.quantities['xmomentum'].get_integral()
3322
3323        import copy
3324        ref_centroid_values =\
3325                 copy.copy(domain.quantities['stage'].centroid_values)
3326
3327        #Test limiter by itself
3328        domain.distribute_to_vertices_and_edges()
3329
3330        #Check that initial limiter doesn't violate cons quan
3331        assert allclose (domain.quantities['stage'].get_integral(),
3332                         initial_volume)
3333        #NOTE: This would fail if any initial stage was less than the
3334        #corresponding bed elevation - but that is reasonable.
3335
3336
3337        #Evolution
3338        for t in domain.evolve(yieldstep = 0.05, finaltime = 10.0):
3339            volume =  domain.quantities['stage'].get_integral()
3340
3341            #print t, volume, initial_volume
3342
3343            assert allclose (volume, initial_volume)
3344
3345
3346        os.remove(domain.get_name() + '.sww')
3347
3348
3349    def test_conservation_5(self):
3350        """Test that momentum is conserved globally in
3351        steady state scenario
3352
3353        This one uses a slopy bed, dirichlet and reflective bdries
3354        """
3355        from mesh_factory import rectangular
3356        from Numeric import array
3357
3358        # Create basic mesh
3359        points, vertices, boundary = rectangular(6, 6)
3360
3361        # Create shallow water domain
3362        domain = Domain(points, vertices, boundary)
3363        domain.smooth = False
3364        domain.default_order = 2
3365
3366        # IC
3367        def x_slope(x, y):
3368            return x/3
3369
3370        domain.set_quantity('elevation', x_slope)
3371        domain.set_quantity('friction', 0)
3372        domain.set_quantity('stage', 0.4) # Steady
3373
3374        # Boundary conditions (reflective everywhere)
3375        Br = Reflective_boundary(domain)
3376        Bleft = Dirichlet_boundary([0.5,0,0])
3377        Bright = Dirichlet_boundary([0.1,0,0])
3378        domain.set_boundary({'left': Bleft, 'right': Bright,
3379                             'top': Br, 'bottom': Br})
3380
3381        domain.check_integrity()
3382
3383        initial_volume = domain.quantities['stage'].get_integral()
3384        initial_xmom = domain.quantities['xmomentum'].get_integral()
3385
3386
3387        # Evolution
3388        for t in domain.evolve(yieldstep = 0.05, finaltime = 15.0):
3389            stage =  domain.quantities['stage'].get_integral()
3390            xmom = domain.quantities['xmomentum'].get_integral()
3391            ymom = domain.quantities['ymomentum'].get_integral()
3392
3393            if allclose(t, 6):  # Steady state reached
3394                steady_xmom = domain.quantities['xmomentum'].get_integral()
3395                steady_ymom = domain.quantities['ymomentum'].get_integral()
3396                steady_stage = domain.quantities['stage'].get_integral()
3397
3398            if t > 6:
3399                #print '%.2f %14.8f %14.8f' %(t, ymom, steady_ymom)
3400                msg = 'xmom=%.2f, steady_xmom=%.2f' %(xmom, steady_xmom)
3401                assert allclose(xmom, steady_xmom), msg
3402                assert allclose(ymom, steady_ymom)
3403                assert allclose(stage, steady_stage)
3404
3405
3406        os.remove(domain.get_name() + '.sww')
3407
3408
3409
3410
3411
3412    def test_conservation_real(self):
3413        """Test that momentum is conserved globally
3414       
3415        Stephen finally made a test that revealed the problem.
3416        This test failed with code prior to 25 July 2005
3417        """
3418
3419        yieldstep = 0.01
3420        finaltime = 0.05
3421        min_depth = 1.0e-2
3422
3423
3424        import sys
3425        from os import sep; sys.path.append('..'+sep+'abstract_2d_finite_volumes')
3426        from mesh_factory import rectangular
3427
3428
3429        #Create shallow water domain
3430        points, vertices, boundary = rectangular(10, 10, len1=500, len2=500)
3431        domain = Domain(points, vertices, boundary)
3432        domain.smooth = False
3433        domain.default_order = 1
3434        domain.minimum_allowed_height = min_depth
3435
3436        # Set initial condition
3437        class Set_IC:
3438            """Set an initial condition with a constant value, for x0<x<x1
3439            """
3440
3441            def __init__(self, x0=0.25, x1=0.5, h=1.0):
3442                self.x0 = x0
3443                self.x1 = x1
3444                self.= h
3445
3446            def __call__(self, x, y):
3447                return self.h*((x>self.x0)&(x<self.x1))
3448
3449
3450        domain.set_quantity('stage', Set_IC(200.0,300.0,5.0))
3451
3452
3453        #Boundaries
3454        R = Reflective_boundary(domain)
3455        domain.set_boundary( {'left': R, 'right': R, 'top':R, 'bottom': R})
3456
3457        ref = domain.quantities['stage'].get_integral()
3458
3459        # Evolution
3460        for t in domain.evolve(yieldstep = yieldstep, finaltime = finaltime):
3461            pass
3462            #print 'Integral stage = ',\
3463            #      domain.quantities['stage'].get_integral(),\
3464            #      ' Time = ',domain.time
3465
3466
3467        now = domain.quantities['stage'].get_integral()
3468
3469        msg = 'Stage not conserved: was %f, now %f' %(ref, now)
3470        assert allclose(ref, now), msg
3471
3472        os.remove(domain.get_name() + '.sww')
3473
3474    def test_second_order_flat_bed_onestep(self):
3475
3476        from mesh_factory import rectangular
3477        from Numeric import array
3478
3479        #Create basic mesh
3480        points, vertices, boundary = rectangular(6, 6)
3481
3482        #Create shallow water domain
3483        domain = Domain(points, vertices, boundary)
3484        domain.smooth = False
3485        domain.default_order = 2
3486        domain.beta_w      = 0.9
3487        domain.beta_w_dry  = 0.9
3488        domain.beta_uh     = 0.9
3489        domain.beta_uh_dry = 0.9
3490        domain.beta_vh     = 0.9
3491        domain.beta_vh_dry = 0.9
3492        domain.H0 = 0 # Backwards compatibility (6/2/7)       
3493       
3494        # Boundary conditions
3495        Br = Reflective_boundary(domain)
3496        Bd = Dirichlet_boundary([0.1, 0., 0.])
3497        domain.set_boundary({'left': Bd, 'right': Br, 'top': Br, 'bottom': Br})
3498
3499        domain.check_integrity()
3500
3501        # Evolution
3502        for t in domain.evolve(yieldstep = 0.05, finaltime = 0.05):
3503            pass# domain.write_time()
3504
3505        # Data from earlier version of abstract_2d_finite_volumes
3506        assert allclose(domain.min_timestep, 0.0396825396825)
3507        assert allclose(domain.max_timestep, 0.0396825396825)
3508
3509        assert allclose(domain.quantities['stage'].centroid_values[:12],
3510                        [0.00171396, 0.02656103, 0.00241523, 0.02656103,
3511                        0.00241523, 0.02656103, 0.00241523, 0.02656103,
3512                        0.00241523, 0.02656103, 0.00241523, 0.0272623])
3513
3514        domain.distribute_to_vertices_and_edges()
3515
3516        assert allclose(domain.quantities['stage'].vertex_values[:12,0],
3517                        [0.0001714, 0.02656103, 0.00024152,
3518                        0.02656103, 0.00024152, 0.02656103,
3519                        0.00024152, 0.02656103, 0.00024152,
3520                        0.02656103, 0.00024152, 0.0272623])
3521
3522        assert allclose(domain.quantities['stage'].vertex_values[:12,1],
3523                        [0.00315012, 0.02656103, 0.00024152, 0.02656103,
3524                         0.00024152, 0.02656103, 0.00024152, 0.02656103,
3525                         0.00024152, 0.02656103, 0.00040506, 0.0272623])
3526
3527        assert allclose(domain.quantities['stage'].vertex_values[:12,2],
3528                        [0.00182037, 0.02656103, 0.00676264,
3529                         0.02656103, 0.00676264, 0.02656103,
3530                         0.00676264, 0.02656103, 0.00676264,
3531                         0.02656103, 0.0065991, 0.0272623])
3532
3533        assert allclose(domain.quantities['xmomentum'].centroid_values[:12],
3534                        [0.00113961, 0.01302432, 0.00148672,
3535                         0.01302432, 0.00148672, 0.01302432,
3536                         0.00148672, 0.01302432, 0.00148672 ,
3537                         0.01302432, 0.00148672, 0.01337143])
3538
3539        assert allclose(domain.quantities['ymomentum'].centroid_values[:12],
3540                        [-2.91240050e-004 , 1.22721531e-004,
3541                         -1.22721531e-004, 1.22721531e-004 ,
3542                         -1.22721531e-004, 1.22721531e-004,
3543                         -1.22721531e-004 , 1.22721531e-004,
3544                         -1.22721531e-004, 1.22721531e-004,
3545                         -1.22721531e-004, -4.57969873e-005])
3546
3547        os.remove(domain.get_name() + '.sww')
3548
3549
3550    def test_second_order_flat_bed_moresteps(self):
3551
3552        from mesh_factory import rectangular
3553        from Numeric import array
3554
3555        #Create basic mesh
3556        points, vertices, boundary = rectangular(6, 6)
3557
3558        #Create shallow water domain
3559        domain = Domain(points, vertices, boundary)
3560        domain.smooth = False
3561        domain.default_order=2
3562
3563        # Boundary conditions
3564        Br = Reflective_boundary(domain)
3565        Bd = Dirichlet_boundary([0.1, 0., 0.])
3566        domain.set_boundary({'left': Bd, 'right': Br, 'top': Br, 'bottom': Br})
3567
3568        domain.check_integrity()
3569
3570        #Evolution
3571        for t in domain.evolve(yieldstep = 0.05, finaltime = 0.1):
3572            pass
3573
3574        #Data from earlier version of abstract_2d_finite_volumes
3575        #assert allclose(domain.min_timestep, 0.0396825396825)
3576        #assert allclose(domain.max_timestep, 0.0396825396825)
3577        #print domain.quantities['stage'].centroid_values
3578
3579        os.remove(domain.get_name() + '.sww')       
3580
3581
3582    def test_flatbed_first_order(self):
3583        from mesh_factory import rectangular
3584        from Numeric import array
3585
3586        #Create basic mesh
3587        N = 8
3588        points, vertices, boundary = rectangular(N, N)
3589
3590        #Create shallow water domain
3591        domain = Domain(points, vertices, boundary)
3592        domain.smooth = False
3593        domain.default_order=1
3594        domain.H0 = 0 # Backwards compatibility (6/2/7)       
3595
3596        # Boundary conditions
3597        Br = Reflective_boundary(domain)
3598        Bd = Dirichlet_boundary([0.2,0.,0.])
3599
3600        domain.set_boundary({'left': Bd, 'right': Br, 'top': Br, 'bottom': Br})
3601        domain.check_integrity()
3602
3603
3604        #Evolution
3605        for t in domain.evolve(yieldstep = 0.02, finaltime = 0.5):
3606            pass
3607            #domain.write_time()
3608
3609        #FIXME: These numbers were from version before 25/10
3610        #assert allclose(domain.min_timestep, 0.0140413643926)
3611        #assert allclose(domain.max_timestep, 0.0140947355753)
3612
3613        for i in range(3):
3614            #assert allclose(domain.quantities['stage'].edge_values[:4,i],
3615            #                [0.10730244,0.12337617,0.11200126,0.12605666])
3616
3617            assert allclose(domain.quantities['xmomentum'].edge_values[:4,i],
3618                            [0.07610894,0.06901572,0.07284461,0.06819712])
3619
3620            #assert allclose(domain.quantities['ymomentum'].edge_values[:4,i],
3621            #                [-0.0060238, -0.00157404, -0.00309633, -0.0001637])
3622
3623
3624        os.remove(domain.get_name() + '.sww')           
3625
3626    def test_flatbed_second_order(self):
3627        from mesh_factory import rectangular
3628        from Numeric import array
3629
3630        #Create basic mesh
3631        N = 8
3632        points, vertices, boundary = rectangular(N, N)
3633
3634        #Create shallow water domain
3635        domain = Domain(points, vertices, boundary)
3636        domain.smooth = False
3637        domain.default_order=2
3638        domain.beta_w      = 0.9
3639        domain.beta_w_dry  = 0.9
3640        domain.beta_uh     = 0.9
3641        domain.beta_uh_dry = 0.9
3642        domain.beta_vh     = 0.9
3643        domain.beta_vh_dry = 0.9       
3644        #domain.minimum_allowed_height = 0.0 #Makes it like the 'oldstyle' balance
3645        domain.H0 = 0 # Backwards compatibility (6/2/7)
3646        domain.use_centroid_velocities = False # Backwards compatibility (8/5/8)
3647        domain.set_maximum_allowed_speed(1.0)       
3648
3649        # Boundary conditions
3650        Br = Reflective_boundary(domain)
3651        Bd = Dirichlet_boundary([0.2,0.,0.])
3652
3653        domain.set_boundary({'left': Bd, 'right': Br, 'top': Br, 'bottom': Br})
3654        domain.check_integrity()
3655
3656        # Evolution
3657        for t in domain.evolve(yieldstep = 0.01, finaltime = 0.03):
3658            pass
3659
3660        msg = 'min step was %f instead of %f' %(domain.min_timestep,
3661                                                0.0210448446782) 
3662
3663        assert allclose(domain.min_timestep, 0.0210448446782), msg
3664        assert allclose(domain.max_timestep, 0.0210448446782)
3665
3666        #print domain.quantities['stage'].vertex_values[:4,0]
3667        #print domain.quantities['xmomentum'].vertex_values[:4,0]
3668        #print domain.quantities['ymomentum'].vertex_values[:4,0]
3669
3670        #FIXME: These numbers were from version before 25/10
3671        #assert allclose(domain.quantities['stage'].vertex_values[:4,0],
3672        #                [0.00101913,0.05352143,0.00104852,0.05354394])
3673
3674        #FIXME: These numbers were from version before 21/3/6 -
3675        #could be recreated by setting maximum_allowed_speed to 0 maybe
3676        #assert allclose(domain.quantities['xmomentum'].vertex_values[:4,0],
3677        #                [ 0.00064835, 0.03685719, 0.00085073, 0.03687313])
3678       
3679        assert allclose(domain.quantities['xmomentum'].vertex_values[:4,0],
3680                        [ 0.00090581, 0.03685719, 0.00088303, 0.03687313])
3681                       
3682                       
3683
3684        #assert allclose(domain.quantities['xmomentum'].vertex_values[:4,0],
3685        #                [0.00090581,0.03685719,0.00088303,0.03687313])
3686
3687        assert allclose(domain.quantities['ymomentum'].vertex_values[:4,0],
3688                        [-0.00139463,0.0006156,-0.00060364,0.00061827])
3689
3690
3691        os.remove(domain.get_name() + '.sww')
3692
3693       
3694    def test_flatbed_second_order_vmax_0(self):
3695        from mesh_factory import rectangular
3696        from Numeric import array
3697
3698        #Create basic mesh
3699        N = 8
3700        points, vertices, boundary = rectangular(N, N)
3701
3702        #Create shallow water domain
3703        domain = Domain(points, vertices, boundary)
3704        domain.smooth = False
3705        domain.default_order=2
3706        domain.beta_w      = 0.9
3707        domain.beta_w_dry  = 0.9
3708        domain.beta_uh     = 0.9
3709        domain.beta_uh_dry = 0.9
3710        domain.beta_vh     = 0.9
3711        domain.beta_vh_dry = 0.9       
3712        domain.maximum_allowed_speed = 0.0 #Makes it like the 'oldstyle'
3713        domain.H0 = 0 # Backwards compatibility (6/2/7)
3714        domain.use_centroid_velocities = False # Backwards compatibility (8/5/8)       
3715
3716        # Boundary conditions
3717        Br = Reflective_boundary(domain)
3718        Bd = Dirichlet_boundary([0.2,0.,0.])
3719
3720        domain.set_boundary({'left': Bd, 'right': Br, 'top': Br, 'bottom': Br})
3721        domain.check_integrity()
3722
3723        #Evolution
3724        for t in domain.evolve(yieldstep = 0.01, finaltime = 0.03):
3725            pass
3726
3727
3728        assert allclose(domain.min_timestep, 0.0210448446782)
3729        assert allclose(domain.max_timestep, 0.0210448446782)
3730
3731        #FIXME: These numbers were from version before 21/3/6 -
3732        #could be recreated by setting maximum_allowed_speed to 0 maybe
3733        assert allclose(domain.quantities['xmomentum'].vertex_values[:4,0],
3734                        [ 0.00064835, 0.03685719, 0.00085073, 0.03687313]) 
3735       
3736
3737        assert allclose(domain.quantities['ymomentum'].vertex_values[:4,0],
3738                        [-0.00139463,0.0006156,-0.00060364,0.00061827])
3739
3740
3741        os.remove(domain.get_name() + '.sww')
3742
3743       
3744
3745    def test_flatbed_second_order_distribute(self):
3746        #Use real data from anuga.abstract_2d_finite_volumes 2
3747        #painfully setup and extracted.
3748        from mesh_factory import rectangular
3749        from Numeric import array
3750
3751        #Create basic mesh
3752        N = 8
3753        points, vertices, boundary = rectangular(N, N)
3754
3755        #Create shallow water domain
3756        domain = Domain(points, vertices, boundary)
3757        domain.smooth = False
3758        domain.default_order=domain._order_=2
3759        domain.beta_w      = 0.9
3760        domain.beta_w_dry  = 0.9
3761        domain.beta_uh     = 0.9
3762        domain.beta_uh_dry = 0.9
3763        domain.beta_vh     = 0.9
3764        domain.beta_vh_dry = 0.9
3765        domain.H0 = 0 # Backwards compatibility (6/2/7)
3766        domain.use_centroid_velocities = False # Backwards compatibility (8/5/8)       
3767        domain.set_maximum_allowed_speed(1.0)       
3768
3769        # Boundary conditions
3770        Br = Reflective_boundary(domain)
3771        Bd = Dirichlet_boundary([0.2,0.,0.])
3772
3773        domain.set_boundary({'left': Bd, 'right': Br, 'top': Br, 'bottom': Br})
3774        domain.check_integrity()
3775
3776
3777
3778        for V in [False, True]:
3779            if V:
3780                #Set centroids as if system had been evolved
3781                L = zeros(2*N*N, Float)
3782                L[:32] = [7.21205592e-003, 5.35214298e-002, 1.00910824e-002,
3783                          5.35439433e-002, 1.00910824e-002, 5.35439433e-002,
3784                          1.00910824e-002, 5.35439433e-002, 1.00910824e-002,
3785                          5.35439433e-002, 1.00910824e-002, 5.35439433e-002,
3786                          1.00910824e-002, 5.35393928e-002, 1.02344264e-002,
3787                          5.59605058e-002, 0.00000000e+000, 3.31027800e-004,
3788                          0.00000000e+000, 4.37962142e-005, 0.00000000e+000,
3789                          4.37962142e-005, 0.00000000e+000, 4.37962142e-005,
3790                          0.00000000e+000, 4.37962142e-005, 0.00000000e+000,
3791                          4.37962142e-005, 0.00000000e+000, 4.37962142e-005,
3792                          0.00000000e+000, 5.57305948e-005]
3793
3794                X = zeros(2*N*N, Float)
3795                X[:32] = [6.48351607e-003, 3.68571894e-002, 8.50733285e-003,
3796                          3.68731327e-002, 8.50733285e-003, 3.68731327e-002,
3797                          8.50733285e-003, 3.68731327e-002, 8.50733285e-003,
3798                          3.68731327e-002, 8.50733285e-003, 3.68731327e-002,
3799                          8.50733285e-003, 3.68693861e-002, 8.65220973e-003,
3800                          3.85055387e-002, 0.00000000e+000, 2.86060840e-004,
3801                          0.00000000e+000, 3.58905503e-005, 0.00000000e+000,
3802                          3.58905503e-005, 0.00000000e+000, 3.58905503e-005,
3803                          0.00000000e+000, 3.58905503e-005, 0.00000000e+000,
3804                          3.58905503e-005, 0.00000000e+000, 3.58905503e-005,
3805                          0.00000000e+000, 4.57662812e-005]
3806
3807                Y = zeros(2*N*N, Float)
3808                Y[:32]=[-1.39463104e-003, 6.15600298e-004, -6.03637382e-004,
3809                        6.18272251e-004, -6.03637382e-004, 6.18272251e-004,
3810                        -6.03637382e-004, 6.18272251e-004, -6.03637382e-004,
3811                        6.18272251e-004, -6.03637382e-004, 6.18272251e-004,
3812                        -6.03637382e-004, 6.18599320e-004, -6.74622797e-004,
3813                        -1.48934756e-004, 0.00000000e+000, -5.35079969e-005,
3814                        0.00000000e+000, -2.57264987e-005, 0.00000000e+000,
3815                        -2.57264987e-005, 0.00000000e+000, -2.57264987e-005,
3816                        0.00000000e+000, -2.57264987e-005, 0.00000000e+000,
3817                        -2.57264987e-005, 0.00000000e+000, -2.57264987e-005,
3818                        0.00000000e+000, -2.57635178e-005]
3819
3820
3821                domain.set_quantity('stage', L, location='centroids')
3822                domain.set_quantity('xmomentum', X, location='centroids')
3823                domain.set_quantity('ymomentum', Y, location='centroids')
3824
3825                domain.check_integrity()
3826            else:
3827                #Evolution
3828                for t in domain.evolve(yieldstep = 0.01, finaltime = 0.03):
3829                    pass
3830                assert allclose(domain.min_timestep, 0.0210448446782)
3831                assert allclose(domain.max_timestep, 0.0210448446782)
3832
3833
3834            #Centroids were correct but not vertices.
3835            #Hence the check of distribute below.
3836            assert allclose(domain.quantities['stage'].centroid_values[:4],
3837                            [0.00721206,0.05352143,0.01009108,0.05354394])
3838
3839            assert allclose(domain.quantities['xmomentum'].centroid_values[:4],
3840                            [0.00648352,0.03685719,0.00850733,0.03687313])
3841
3842            assert allclose(domain.quantities['ymomentum'].centroid_values[:4],
3843                            [-0.00139463,0.0006156,-0.00060364,0.00061827])
3844
3845            #print 'C17=', domain.quantities['xmomentum'].centroid_values[17]
3846            #print 'C19=', domain.quantities['xmomentum'].centroid_values[19]
3847
3848            #assert allclose(domain.quantities['xmomentum'].centroid_values[17],0.00028606084)
3849            ##print domain.quantities['xmomentum'].centroid_values[17], V
3850            ##print
3851            if not V:
3852                #FIXME: These numbers were from version before 21/3/6 -
3853                #could be recreated by setting maximum_allowed_speed to 0 maybe           
3854               
3855                #assert allclose(domain.quantities['xmomentum'].centroid_values[17], 0.0)               
3856                assert allclose(domain.quantities['xmomentum'].centroid_values[17], 0.000286060839592)                         
3857
3858            else:
3859                assert allclose(domain.quantities['xmomentum'].centroid_values[17], 0.00028606084)
3860
3861            import copy
3862            XX = copy.copy(domain.quantities['xmomentum'].centroid_values)
3863            assert allclose(domain.quantities['xmomentum'].centroid_values, XX)
3864
3865            domain.distribute_to_vertices_and_edges()
3866
3867            #assert allclose(domain.quantities['xmomentum'].centroid_values, XX)
3868
3869            #assert allclose(domain.quantities['xmomentum'].centroid_values[17],
3870            #                0.0)
3871            assert allclose(domain.quantities['xmomentum'].centroid_values[17], 0.000286060839592)                                                 
3872
3873
3874            #FIXME: These numbers were from version before 25/10
3875            #assert allclose(domain.quantities['stage'].vertex_values[:4,0],
3876            #                [0.00101913,0.05352143,0.00104852,0.05354394])
3877
3878            assert allclose(domain.quantities['ymomentum'].vertex_values[:4,0],
3879                            [-0.00139463,0.0006156,-0.00060364,0.00061827])
3880
3881
3882            assert allclose(domain.quantities['xmomentum'].vertex_values[:4,0],
3883                            [0.00090581,0.03685719,0.00088303,0.03687313])
3884
3885
3886            #NB NO longer relvant:
3887
3888            #This was the culprit. First triangles vertex 0 had an
3889            #x-momentum of 0.0064835 instead of 0.00090581 and
3890            #third triangle had 0.00850733 instead of 0.00088303
3891            #print domain.quantities['xmomentum'].vertex_values[:4,0]
3892
3893            #print domain.quantities['xmomentum'].vertex_values[:4,0]
3894            #assert allclose(domain.quantities['xmomentum'].vertex_values[:4,0],
3895            #                [0.00090581,0.03685719,0.00088303,0.03687313])
3896
3897        os.remove(domain.get_name() + '.sww')
3898
3899
3900
3901    def test_bedslope_problem_first_order(self):
3902
3903        from mesh_factory import rectangular
3904        from Numeric import array
3905
3906        #Create basic mesh
3907        points, vertices, boundary = rectangular(6, 6)
3908
3909        #Create shallow water domain
3910        domain = Domain(points, vertices, boundary)
3911        domain.smooth = False
3912        domain.default_order = 1
3913
3914        #Bed-slope and friction
3915        def x_slope(x, y):
3916            return -x/3
3917
3918        domain.set_quantity('elevation', x_slope)
3919
3920        # Boundary conditions
3921        Br = Reflective_boundary(domain)
3922        domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
3923
3924        #Initial condition
3925        #domain.set_quantity('stage', Constant_height(x_slope, 0.05))
3926        domain.set_quantity('stage', expression='elevation+0.05')
3927        domain.check_integrity()
3928
3929        #Evolution
3930        for t in domain.evolve(yieldstep = 0.05, finaltime = 0.05):
3931            pass# domain.write_time()
3932
3933        # FIXME (Ole): Need some other assertion here!
3934        #print domain.min_timestep, domain.max_timestep   
3935        #assert allclose(domain.min_timestep, 0.050010003001)
3936        #assert allclose(domain.max_timestep, 0.050010003001)
3937
3938
3939        os.remove(domain.get_name() + '.sww')
3940       
3941    def test_bedslope_problem_first_order_moresteps(self):
3942
3943        from mesh_factory import rectangular
3944        from Numeric import array
3945
3946        #Create basic mesh
3947        points, vertices, boundary = rectangular(6, 6)
3948
3949        #Create shallow water domain
3950        domain = Domain(points, vertices, boundary)
3951        domain.smooth = False
3952        domain.default_order = 1
3953        domain.beta_h = 0.0 # Use first order in h-limiter
3954       
3955        # FIXME (Ole): Need tests where these two are commented out
3956        domain.H0 = 0        # Backwards compatibility (6/2/7)       
3957        domain.tight_slope_limiters = 0 # Backwards compatibility (14/4/7)
3958        domain.use_centroid_velocities = 0 # Backwards compatibility (7/5/8)                 
3959
3960        #Bed-slope and friction
3961        def x_slope(x, y):
3962            return -x/3
3963
3964        domain.set_quantity('elevation', x_slope)
3965
3966        # Boundary conditions
3967        Br = Reflective_boundary(domain)
3968        domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
3969
3970        #Initial condition
3971        domain.set_quantity('stage', expression='elevation+0.05')
3972        domain.check_integrity()
3973
3974        #Evolution
3975        for t in domain.evolve(yieldstep = 0.05, finaltime = 0.5):
3976            pass# domain.write_time()
3977
3978        #Data from earlier version of abstract_2d_finite_volumes
3979        #print domain.quantities['stage'].centroid_values
3980
3981        assert allclose(domain.quantities['stage'].centroid_values,
3982                        [-0.02998628, -0.01520652, -0.03043492,
3983                        -0.0149132, -0.03004706, -0.01476251,
3984                        -0.0298215, -0.01467976, -0.02988158,
3985                        -0.01474662, -0.03036161, -0.01442995,
3986                        -0.07624583, -0.06297061, -0.07733792,
3987                        -0.06342237, -0.07695439, -0.06289595,
3988                        -0.07635559, -0.0626065, -0.07633628,
3989                        -0.06280072, -0.07739632, -0.06386738,
3990                        -0.12161738, -0.11028239, -0.1223796,
3991                        -0.11095953, -0.12189744, -0.11048616,
3992                        -0.12074535, -0.10987605, -0.12014311,
3993                        -0.10976691, -0.12096859, -0.11087692,
3994                        -0.16868259, -0.15868061, -0.16801135,
3995                        -0.1588003, -0.16674343, -0.15813323,
3996                        -0.16457595, -0.15693826, -0.16281096,
3997                        -0.15585154, -0.16283873, -0.15540068,
3998                        -0.17450362, -0.19919913, -0.18062882,
3999                        -0.19764131, -0.17783111, -0.19407213,
4000                        -0.1736915, -0.19053624, -0.17228678,
4001                        -0.19105634, -0.17920133, -0.1968828,
4002                        -0.14244395, -0.14604641, -0.14473537,
4003                        -0.1506107, -0.14510055, -0.14919522,
4004                        -0.14175896, -0.14560798, -0.13911658,
4005                        -0.14439383, -0.13924047, -0.14829043])
4006
4007        os.remove(domain.get_name() + '.sww')
4008       
4009    def test_bedslope_problem_second_order_one_step(self):
4010
4011        from mesh_factory import rectangular
4012        from Numeric import array
4013
4014        #Create basic mesh
4015        points, vertices, boundary = rectangular(6, 6)
4016
4017        #Create shallow water domain
4018        domain = Domain(points, vertices, boundary)
4019        domain.smooth = False
4020        domain.default_order=2
4021        domain.beta_w      = 0.9
4022        domain.beta_w_dry  = 0.9
4023        domain.beta_uh     = 0.9
4024        domain.beta_uh_dry = 0.9
4025        domain.beta_vh     = 0.9
4026        domain.beta_vh_dry = 0.9
4027
4028       
4029        # FIXME (Ole): Need tests where this is commented out
4030        domain.tight_slope_limiters = 0 # Backwards compatibility (14/4/7)
4031        domain.use_centroid_velocities = 0 # Backwards compatibility (7/5/8)     
4032       
4033        #Bed-slope and friction at vertices (and interpolated elsewhere)
4034        def x_slope(x, y):
4035            return -x/3
4036
4037        domain.set_quantity('elevation', x_slope)
4038
4039        # Boundary conditions
4040        Br = Reflective_boundary(domain)
4041        domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
4042
4043        #Initial condition
4044        domain.set_quantity('stage', expression='elevation+0.05')
4045        domain.check_integrity()
4046
4047        assert allclose(domain.quantities['stage'].centroid_values,
4048                        [0.01296296, 0.03148148, 0.01296296,
4049                        0.03148148, 0.01296296, 0.03148148,
4050                        0.01296296, 0.03148148, 0.01296296,
4051                        0.03148148, 0.01296296, 0.03148148,
4052                        -0.04259259, -0.02407407, -0.04259259,
4053                        -0.02407407, -0.04259259, -0.02407407,
4054                        -0.04259259, -0.02407407, -0.04259259,
4055                        -0.02407407, -0.04259259, -0.02407407,
4056                        -0.09814815, -0.07962963, -0.09814815,
4057                        -0.07962963, -0.09814815, -0.07962963,
4058                        -0.09814815, -0.07962963, -0.09814815,
4059                        -0.07962963, -0.09814815, -0.07962963,
4060                        -0.1537037 , -0.13518519, -0.1537037,
4061                        -0.13518519, -0.1537037, -0.13518519,
4062                        -0.1537037 , -0.13518519, -0.1537037,
4063                        -0.13518519, -0.1537037, -0.13518519,
4064                        -0.20925926, -0.19074074, -0.20925926,
4065                        -0.19074074, -0.20925926, -0.19074074,
4066                        -0.20925926, -0.19074074, -0.20925926,
4067                        -0.19074074, -0.20925926, -0.19074074,
4068                        -0.26481481, -0.2462963, -0.26481481,
4069                        -0.2462963, -0.26481481, -0.2462963,
4070                        -0.26481481, -0.2462963, -0.26481481,
4071                        -0.2462963, -0.26481481, -0.2462963])
4072
4073
4074        #print domain.quantities['stage'].extrapolate_second_order()
4075        #domain.distribute_to_vertices_and_edges()
4076        #print domain.quantities['stage'].vertex_values[:,0]
4077
4078        #Evolution
4079        for t in domain.evolve(yieldstep = 0.05, finaltime = 0.05):
4080            #domain.write_time()
4081            pass
4082
4083
4084        #print domain.quantities['stage'].centroid_values
4085        assert allclose(domain.quantities['stage'].centroid_values,
4086                 [0.01290985, 0.02356019, 0.01619096, 0.02356019, 0.01619096,
4087                  0.02356019, 0.01619096, 0.02356019, 0.01619096, 0.02356019,
4088                  0.01619096, 0.0268413, -0.04411074, -0.0248011, -0.04186556,
4089                  -0.0248011, -0.04186556, -0.0248011, -0.04186556, -0.0248011,
4090                  -0.04186556, -0.0248011, -0.04186556, -0.02255593,
4091                  -0.09966629, -0.08035666, -0.09742112, -0.08035666,
4092                  -0.09742112, -0.08035666, -0.09742112, -0.08035666,
4093                  -0.09742112, -0.08035666, -0.09742112, -0.07811149,
4094                  -0.15522185, -0.13591222, -0.15297667, -0.13591222,
4095                  -0.15297667, -0.13591222, -0.15297667, -0.13591222,
4096                  -0.15297667, -0.13591222, -0.15297667, -0.13366704,
4097                  -0.2107774, -0.19146777, -0.20853223, -0.19146777,
4098                  -0.20853223, -0.19146777, -0.20853223, -0.19146777,
4099                  -0.20853223, -0.19146777, -0.20853223, -0.1892226,
4100                  -0.26120669, -0.24776246, -0.25865535, -0.24776246,
4101                  -0.25865535, -0.24776246, -0.25865535, -0.24776246,
4102                  -0.25865535, -0.24776246, -0.25865535, -0.24521113])
4103
4104        os.remove(domain.get_name() + '.sww')
4105
4106    def test_bedslope_problem_second_order_two_steps(self):
4107
4108        from mesh_factory import rectangular
4109        from Numeric import array
4110
4111        #Create basic mesh
4112        points, vertices, boundary = rectangular(6, 6)
4113
4114        #Create shallow water domain
4115        domain = Domain(points, vertices, boundary)
4116        domain.smooth = False
4117        domain.default_order=2
4118        domain.beta_w      = 0.9
4119        domain.beta_w_dry  = 0.9
4120        domain.beta_uh     = 0.9
4121        domain.beta_uh_dry = 0.9
4122        domain.beta_vh     = 0.9
4123        domain.beta_vh_dry = 0.9
4124        domain.beta_h = 0.0 #Use first order in h-limiter
4125       
4126        # FIXME (Ole): Need tests where this is commented out
4127        domain.tight_slope_limiters = 0 # Backwards compatibility (14/4/7)                 
4128        domain.H0 = 0 # Backwards compatibility (6/2/7)
4129        domain.use_centroid_velocities = 0 # Backwards compatibility (7/5/8)
4130       
4131
4132        #Bed-slope and friction at vertices (and interpolated elsewhere)
4133        def x_slope(x, y):
4134            return -x/3
4135
4136        domain.set_quantity('elevation', x_slope)
4137
4138        # Boundary conditions
4139        Br = Reflective_boundary(domain)
4140        domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
4141
4142        #Initial condition
4143        domain.set_quantity('stage', expression='elevation+0.05')
4144        domain.check_integrity()
4145
4146        assert allclose(domain.quantities['stage'].centroid_values,
4147                        [0.01296296, 0.03148148, 0.01296296,
4148                        0.03148148, 0.01296296, 0.03148148,
4149                        0.01296296, 0.03148148, 0.01296296,
4150                        0.03148148, 0.01296296, 0.03148148,
4151                        -0.04259259, -0.02407407, -0.04259259,
4152                        -0.02407407, -0.04259259, -0.02407407,
4153                        -0.04259259, -0.02407407, -0.04259259,
4154                        -0.02407407, -0.04259259, -0.02407407,
4155                        -0.09814815, -0.07962963, -0.09814815,
4156                        -0.07962963, -0.09814815, -0.07962963,
4157                        -0.09814815, -0.07962963, -0.09814815,
4158                        -0.07962963, -0.09814815, -0.07962963,
4159                        -0.1537037 , -0.13518519, -0.1537037,
4160                        -0.13518519, -0.1537037, -0.13518519,
4161                        -0.1537037 , -0.13518519, -0.1537037,
4162                        -0.13518519, -0.1537037, -0.13518519,
4163                        -0.20925926, -0.19074074, -0.20925926,
4164                        -0.19074074, -0.20925926, -0.19074074,
4165                        -0.20925926, -0.19074074, -0.20925926,
4166                        -0.19074074, -0.20925926, -0.19074074,
4167                        -0.26481481, -0.2462963, -0.26481481,
4168                        -0.2462963, -0.26481481, -0.2462963,
4169                        -0.26481481, -0.2462963, -0.26481481,
4170                        -0.2462963, -0.26481481, -0.2462963])
4171
4172
4173        #print domain.quantities['stage'].extrapolate_second_order()
4174        #domain.distribute_to_vertices_and_edges()
4175        #print domain.quantities['stage'].vertex_values[:,0]
4176
4177        #Evolution
4178        for t in domain.evolve(yieldstep = 0.05, finaltime = 0.1):
4179            pass
4180
4181
4182        #Data from earlier version of abstract_2d_finite_volumes ft=0.1
4183        assert allclose(domain.min_timestep, 0.0376895634803)
4184        assert allclose(domain.max_timestep, 0.0415635655309)
4185
4186
4187        assert allclose(domain.quantities['stage'].centroid_values,
4188                        [0.00855788, 0.01575204, 0.00994606, 0.01717072,
4189                         0.01005985, 0.01716362, 0.01005985, 0.01716299,
4190                         0.01007098, 0.01736248, 0.01216452, 0.02026776,
4191                         -0.04462374, -0.02479045, -0.04199789, -0.0229465,
4192                         -0.04184033, -0.02295693, -0.04184013, -0.02295675,
4193                         -0.04184486, -0.0228168, -0.04028876, -0.02036486,
4194                         -0.10029444, -0.08170809, -0.09772846, -0.08021704,
4195                         -0.09760006, -0.08022143, -0.09759984, -0.08022124,
4196                         -0.09760261, -0.08008893, -0.09603914, -0.07758209,
4197                         -0.15584152, -0.13723138, -0.15327266, -0.13572906,
4198                         -0.15314427, -0.13573349, -0.15314405, -0.13573331,
4199                         -0.15314679, -0.13560104, -0.15158523, -0.13310701,
4200                         -0.21208605, -0.19283913, -0.20955631, -0.19134189,
4201                         -0.20942821, -0.19134598, -0.20942799, -0.1913458,
4202                         -0.20943005, -0.19120952, -0.20781177, -0.18869401,
4203                         -0.25384082, -0.2463294, -0.25047649, -0.24464654,
4204                         -0.25031159, -0.24464253, -0.25031112, -0.24464253,
4205                         -0.25031463, -0.24454764, -0.24885323, -0.24286438])
4206
4207
4208        os.remove(domain.get_name() + '.sww')
4209
4210    def test_bedslope_problem_second_order_two_yieldsteps(self):
4211
4212        from mesh_factory import rectangular
4213        from Numeric import array
4214
4215        #Create basic mesh
4216        points, vertices, boundary = rectangular(6, 6)
4217
4218        #Create shallow water domain
4219        domain = Domain(points, vertices, boundary)
4220        domain.smooth = False
4221        domain.default_order=2
4222        domain.beta_w      = 0.9
4223        domain.beta_w_dry  = 0.9
4224        domain.beta_uh     = 0.9
4225        domain.beta_uh_dry = 0.9
4226        domain.beta_vh     = 0.9
4227        domain.beta_vh_dry = 0.9
4228        domain.beta_h = 0.0 #Use first order in h-limiter
4229       
4230        # FIXME (Ole): Need tests where this is commented out
4231        domain.tight_slope_limiters = 0 # Backwards compatibility (14/4/7)                 
4232        domain.H0 = 0 # Backwards compatibility (6/2/7)
4233        domain.use_centroid_velocities = 0 # Backwards compatibility (7/5/8)
4234       
4235
4236        #Bed-slope and friction at vertices (and interpolated elsewhere)
4237        def x_slope(x, y):
4238            return -x/3
4239
4240        domain.set_quantity('elevation', x_slope)
4241
4242        # Boundary conditions
4243        Br = Reflective_boundary(domain)
4244        domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
4245
4246        #Initial condition
4247        domain.set_quantity('stage', expression='elevation+0.05')
4248        domain.check_integrity()
4249
4250        assert allclose(domain.quantities['stage'].centroid_values,
4251                        [0.01296296, 0.03148148, 0.01296296,
4252                        0.03148148, 0.01296296, 0.03148148,
4253                        0.01296296, 0.03148148, 0.01296296,
4254                        0.03148148, 0.01296296, 0.03148148,
4255                        -0.04259259, -0.02407407, -0.04259259,
4256                        -0.02407407, -0.04259259, -0.02407407,
4257                        -0.04259259, -0.02407407, -0.04259259,
4258                        -0.02407407, -0.04259259, -0.02407407,
4259                        -0.09814815, -0.07962963, -0.09814815,
4260                        -0.07962963, -0.09814815, -0.07962963,
4261                        -0.09814815, -0.07962963, -0.09814815,
4262                        -0.07962963, -0.09814815, -0.07962963,
4263                        -0.1537037 , -0.13518519, -0.1537037,
4264                        -0.13518519, -0.1537037, -0.13518519,
4265                        -0.1537037 , -0.13518519, -0.1537037,
4266                        -0.13518519, -0.1537037, -0.13518519,
4267                        -0.20925926, -0.19074074, -0.20925926,
4268                        -0.19074074, -0.20925926, -0.19074074,
4269                        -0.20925926, -0.19074074, -0.20925926,
4270                        -0.19074074, -0.20925926, -0.19074074,
4271                        -0.26481481, -0.2462963, -0.26481481,
4272                        -0.2462963, -0.26481481, -0.2462963,
4273                        -0.26481481, -0.2462963, -0.26481481,
4274                        -0.2462963, -0.26481481, -0.2462963])
4275
4276
4277        #print domain.quantities['stage'].extrapolate_second_order()
4278        #domain.distribute_to_vertices_and_edges()
4279        #print domain.quantities['stage'].vertex_values[:,0]
4280
4281        #Evolution
4282        for t in domain.evolve(yieldstep = 0.05, finaltime = 0.1):   #0.05??
4283            #domain.write_time()
4284            pass
4285
4286
4287
4288        assert allclose(domain.quantities['stage'].centroid_values,
4289                 [0.00855788, 0.01575204, 0.00994606, 0.01717072, 0.01005985,
4290                  0.01716362, 0.01005985, 0.01716299, 0.01007098, 0.01736248,
4291                  0.01216452, 0.02026776, -0.04462374, -0.02479045, -0.04199789,
4292                  -0.0229465, -0.04184033, -0.02295693, -0.04184013,
4293                  -0.02295675, -0.04184486, -0.0228168, -0.04028876,
4294                  -0.02036486, -0.10029444, -0.08170809, -0.09772846,
4295                  -0.08021704, -0.09760006, -0.08022143, -0.09759984,
4296                  -0.08022124, -0.09760261, -0.08008893, -0.09603914,
4297                  -0.07758209, -0.15584152, -0.13723138, -0.15327266,
4298                  -0.13572906, -0.15314427, -0.13573349, -0.15314405,
4299                  -0.13573331, -0.15314679, -0.13560104, -0.15158523,
4300                  -0.13310701, -0.21208605, -0.19283913, -0.20955631,
4301                  -0.19134189, -0.20942821, -0.19134598, -0.20942799,
4302                  -0.1913458, -0.20943005, -0.19120952, -0.20781177,
4303                  -0.18869401, -0.25384082, -0.2463294, -0.25047649,
4304                  -0.24464654, -0.25031159, -0.24464253, -0.25031112,
4305                  -0.24464253, -0.25031463, -0.24454764, -0.24885323,
4306                  -0.24286438])
4307
4308        os.remove(domain.get_name() + '.sww')
4309
4310    def test_bedslope_problem_second_order_more_steps(self):
4311
4312        from mesh_factory import rectangular
4313        from Numeric import array
4314
4315        #Create basic mesh
4316        points, vertices, boundary = rectangular(6, 6)
4317
4318        #Create shallow water domain
4319        domain = Domain(points, vertices, boundary)
4320        domain.smooth = False
4321        domain.default_order=2
4322        domain.beta_w      = 0.9
4323        domain.beta_w_dry  = 0.9
4324        domain.beta_uh     = 0.9
4325        domain.beta_uh_dry = 0.9
4326        domain.beta_vh     = 0.9
4327        domain.beta_vh_dry = 0.9
4328        domain.beta_h = 0.0 #Use first order in h-limiter
4329       
4330       
4331        # FIXME (Ole): Need tests where these two are commented out
4332        domain.H0 = 0        # Backwards compatibility (6/2/7)       
4333        domain.tight_slope_limiters = 0 # Backwards compatibility (14/4/7)                 
4334        domain.use_centroid_velocities = 0 # Backwards compatibility (7/5/8)
4335       
4336               
4337
4338        #Bed-slope and friction at vertices (and interpolated elsewhere)
4339        def x_slope(x, y):
4340            return -x/3
4341
4342        domain.set_quantity('elevation', x_slope)
4343
4344        # Boundary conditions
4345        Br = Reflective_boundary(domain)
4346        domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
4347
4348        #Initial condition
4349        domain.set_quantity('stage', expression = 'elevation + 0.05')
4350        domain.check_integrity()
4351
4352        assert allclose(domain.quantities['stage'].centroid_values,
4353                        [0.01296296, 0.03148148, 0.01296296,
4354                        0.03148148, 0.01296296, 0.03148148,
4355                        0.01296296, 0.03148148, 0.01296296,
4356                        0.03148148, 0.01296296, 0.03148148,
4357                        -0.04259259, -0.02407407, -0.04259259,
4358                        -0.02407407, -0.04259259, -0.02407407,
4359                        -0.04259259, -0.02407407, -0.04259259,
4360                        -0.02407407, -0.04259259, -0.02407407,
4361                        -0.09814815, -0.07962963, -0.09814815,
4362                        -0.07962963, -0.09814815, -0.07962963,
4363                        -0.09814815, -0.07962963, -0.09814815,
4364                        -0.07962963, -0.09814815, -0.07962963,
4365                        -0.1537037 , -0.13518519, -0.1537037,
4366                        -0.13518519, -0.1537037, -0.13518519,
4367                        -0.1537037 , -0.13518519, -0.1537037,
4368                        -0.13518519, -0.1537037, -0.13518519,
4369                        -0.20925926, -0.19074074, -0.20925926,
4370                        -0.19074074, -0.20925926, -0.19074074,
4371                        -0.20925926, -0.19074074, -0.20925926,
4372                        -0.19074074, -0.20925926, -0.19074074,
4373                        -0.26481481, -0.2462963, -0.26481481,
4374                        -0.2462963, -0.26481481, -0.2462963,
4375                        -0.26481481, -0.2462963, -0.26481481,
4376                        -0.2462963, -0.26481481, -0.2462963])
4377
4378
4379        #print domain.quantities['stage'].extrapolate_second_order()
4380        #domain.distribute_to_vertices_and_edges()
4381        #print domain.quantities['stage'].vertex_values[:,0]
4382
4383        #Evolution
4384        for t in domain.evolve(yieldstep = 0.05, finaltime = 0.5):
4385
4386            # Check that diagnostics works
4387            msg = domain.timestepping_statistics(track_speeds=True)
4388            #FIXME(Ole): One might check the contents of msg here.
4389
4390
4391
4392        assert allclose(domain.quantities['stage'].centroid_values,
4393     [-0.02907028, -0.01475478, -0.02973417, -0.01447186, -0.02932665, -0.01428285,
4394      -0.02901975, -0.0141361,  -0.02898816, -0.01418135, -0.02961409, -0.01403487,
4395      -0.07597998, -0.06252591, -0.07664854, -0.06312532, -0.07638287, -0.06265139,
4396      -0.07571145, -0.06235231, -0.0756817,  -0.06245309, -0.07652292, -0.06289946,
4397      -0.12367464, -0.11088981, -0.12237277, -0.11115338, -0.1218934,  -0.1107174,
4398      -0.12081485, -0.11000491, -0.12038451, -0.11010335, -0.12102113, -0.11012105,
4399      -0.16909116, -0.15831543, -0.16730214, -0.15786249, -0.1665493,  -0.15697919,
4400      -0.16496618, -0.15559852, -0.16338679, -0.15509088, -0.16364092, -0.15424423,
4401      -0.18771107, -0.19903904, -0.18903759, -0.19858437, -0.18701552, -0.19697797,
4402      -0.1833593,  -0.19505871, -0.1818806,  -0.19418042, -0.18586159, -0.19576946,
4403      -0.13986873, -0.14170053, -0.14132188, -0.14560674, -0.14095617, -0.14373292,
4404      -0.13785933, -0.14033364, -0.13592955, -0.13936356, -0.13596008, -0.14216296])
4405
4406        assert allclose(domain.quantities['xmomentum'].centroid_values,
4407     [ 0.00831121,  0.00317948,  0.00731797,  0.00334939,  0.00764717,  0.00348053,
4408       0.00788729,  0.00356522,  0.00780649,  0.00341919,  0.00693525,  0.00310375,
4409       0.02166196,  0.01421475,  0.02017737,  0.01316839,  0.02037015,  0.01368659,
4410       0.02106,     0.01399161,  0.02074514,  0.01354935,  0.01887407,  0.0123113,
4411       0.03775083,  0.02855197,  0.03689337,  0.02759782,  0.03732848,  0.02812072,
4412       0.03872545,  0.02913348,  0.03880939,  0.02803804,  0.03546499,  0.0260039,
4413       0.0632131,   0.04730634,  0.0576324,   0.04592336,  0.05790921,  0.04690514,
4414       0.05986467,  0.04871165,  0.06170068,  0.04811572,  0.05657041,  0.04416292,
4415       0.08489642,  0.07188097,  0.07835261,  0.06843406,  0.07986412,  0.0698247,
4416       0.08201071,  0.07216756,  0.08378418,  0.07273624,  0.080399,    0.06645841,
4417       0.01631548,  0.04691608,  0.0206632,   0.044409,    0.02115518,  0.04560305,
4418       0.02160608,  0.04663725,  0.02174734,  0.04795559,  0.02281427,  0.05667111])
4419
4420
4421        assert allclose(domain.quantities['ymomentum'].centroid_values,
4422     [ 1.45876601e-004, -3.24627393e-004, -1.57572719e-004, -2.92790187e-004,
4423      -9.90988382e-005, -3.06677335e-004, -1.62493106e-004, -3.71310004e-004,
4424      -1.99445058e-004, -3.28493467e-004,  6.68217349e-005, -8.42042805e-006,
4425       5.05093371e-004, -1.42842214e-004, -6.81454718e-005, -5.02084057e-004,
4426      -8.50583861e-005, -4.65443981e-004, -1.96406564e-004, -5.88889562e-004,
4427      -2.70160173e-004, -5.35485454e-004,  2.60780997e-004,  3.12145471e-005,
4428       5.16189608e-004,  1.07069062e-004,  9.29989252e-005, -3.71211119e-004,
4429       1.16350246e-004, -3.82407830e-004, -1.62077969e-004, -6.30906636e-004,
4430      -4.74025708e-004, -6.94463009e-004,  6.15092843e-005,  2.22106820e-004,
4431      -6.29589294e-004,  2.43611937e-004, -5.88125094e-004, -6.94293192e-005,
4432      -4.17914641e-004,  6.64609019e-005, -7.68334577e-004, -3.40232101e-004,
4433      -1.67424308e-003, -7.39485066e-004, -1.59966988e-003,  5.68262838e-005,
4434      -1.48470633e-003, -1.84554882e-003, -2.27200099e-003, -1.67506848e-003,
4435      -1.95610258e-003, -1.47638801e-003, -1.73779477e-003, -1.85498791e-003,
4436      -2.01357843e-003, -2.17675471e-003, -1.65783870e-003, -1.15818681e-003,
4437      -1.18663036e-003, -2.94229849e-003, -3.59309018e-003, -5.13496584e-003,
4438      -6.17359400e-003, -5.98761937e-003, -6.00540116e-003, -5.01121966e-003,
4439      -4.50964850e-003, -3.06319963e-003,  6.08950810e-004, -4.79537921e-004])
4440
4441        os.remove(domain.get_name() + '.sww')
4442
4443
4444
4445    def NOtest_bedslope_problem_second_order_more_steps_feb_2007(self):
4446        """test_bedslope_problem_second_order_more_steps_feb_2007
4447
4448        Test shallow water finite volumes, using parameters from
4449        feb 2007 rather than backward compatibility ad infinitum
4450       
4451        """
4452        from mesh_factory import rectangular
4453        from Numeric import array
4454
4455        #Create basic mesh
4456        points, vertices, boundary = rectangular(6, 6)
4457
4458        #Create shallow water domain
4459        domain = Domain(points, vertices, boundary)
4460        domain.smooth = False
4461        domain.default_order = 2
4462        domain.beta_w      = 0.9
4463        domain.beta_w_dry  = 0.9
4464        domain.beta_uh     = 0.9
4465        domain.beta_uh_dry = 0.9
4466        domain.beta_vh     = 0.9
4467        domain.beta_vh_dry = 0.9
4468        domain.beta_h = 0.0 #Use first order in h-limiter
4469        domain.H0 = 0.001
4470        domain.tight_slope_limiters = 1
4471
4472        #Bed-slope and friction at vertices (and interpolated elsewhere)
4473        def x_slope(x, y):
4474            return -x/3
4475
4476        domain.set_quantity('elevation', x_slope)
4477
4478        # Boundary conditions
4479        Br = Reflective_boundary(domain)
4480        domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
4481
4482        #Initial condition
4483        domain.set_quantity('stage', expression = 'elevation + 0.05')
4484        domain.check_integrity()
4485
4486        assert allclose(domain.quantities['stage'].centroid_values,
4487                        [0.01296296, 0.03148148, 0.01296296,
4488                        0.03148148, 0.01296296, 0.03148148,
4489                        0.01296296, 0.03148148, 0.01296296,
4490                        0.03148148, 0.01296296, 0.03148148,
4491                        -0.04259259, -0.02407407, -0.04259259,
4492                        -0.02407407, -0.04259259, -0.02407407,
4493                        -0.04259259, -0.02407407, -0.04259259,
4494                        -0.02407407, -0.04259259, -0.02407407,
4495                        -0.09814815, -0.07962963, -0.09814815,
4496                        -0.07962963, -0.09814815, -0.07962963,
4497                        -0.09814815, -0.07962963, -0.09814815,
4498                        -0.07962963, -0.09814815, -0.07962963,
4499                        -0.1537037 , -0.13518519, -0.1537037,
4500                        -0.13518519, -0.1537037, -0.13518519,
4501                        -0.1537037 , -0.13518519, -0.1537037,
4502                        -0.13518519, -0.1537037, -0.13518519,
4503                        -0.20925926, -0.19074074, -0.20925926,
4504                        -0.19074074, -0.20925926, -0.19074074,
4505                        -0.20925926, -0.19074074, -0.20925926,
4506                        -0.19074074, -0.20925926, -0.19074074,
4507                        -0.26481481, -0.2462963, -0.26481481,
4508                        -0.2462963, -0.26481481, -0.2462963,
4509                        -0.26481481, -0.2462963, -0.26481481,
4510                        -0.2462963, -0.26481481, -0.2462963])
4511
4512
4513        #print domain.quantities['stage'].extrapolate_second_order()
4514        #domain.distribute_to_vertices_and_edges()
4515        #print domain.quantities['stage'].vertex_values[:,0]
4516
4517        #Evolution
4518        for t in domain.evolve(yieldstep = 0.05, finaltime = 0.5):
4519            pass
4520
4521
4522        #print domain.quantities['stage'].centroid_values
4523           
4524        assert allclose(domain.quantities['stage'].centroid_values,
4525         [-0.03348416, -0.01749303, -0.03299091, -0.01739241, -0.03246447, -0.01732016,
4526          -0.03205390, -0.01717833, -0.03146383, -0.01699831, -0.03076577, -0.01671795,
4527          -0.07952656, -0.06684763, -0.07721455, -0.06668388, -0.07632976, -0.06600113,
4528          -0.07523678, -0.06546373, -0.07447040, -0.06508861, -0.07438723, -0.06359288,
4529          -0.12526729, -0.11205668, -0.12179433, -0.11068104, -0.12048395, -0.10968948,
4530          -0.11912023, -0.10862628, -0.11784090, -0.10803744, -0.11790629, -0.10742354,
4531          -0.16859613, -0.15427413, -0.16664444, -0.15464452, -0.16570816, -0.15327556,
4532          -0.16409162, -0.15204092, -0.16264608, -0.15102139, -0.16162736, -0.14969205,
4533          -0.18736511, -0.19874036, -0.18811230, -0.19758289, -0.18590182, -0.19580301,
4534          -0.18234588, -0.19423215, -0.18100376, -0.19380116, -0.18509710, -0.19501636,
4535          -0.13982382, -0.14166819, -0.14132775, -0.14528694, -0.14096905, -0.14351126,
4536          -0.13800356, -0.14027920, -0.13613538, -0.13936795, -0.13621902, -0.14204982])
4537
4538                     
4539        assert allclose(domain.quantities['xmomentum'].centroid_values,
4540      [0.00600290,  0.00175780,  0.00591905,  0.00190903,  0.00644462,  0.00203095,
4541       0.00684561,  0.00225089,  0.00708208,  0.00236235,  0.00649095,  0.00222343,
4542       0.02068693,  0.01164034,  0.01983343,  0.01159526,  0.02044611,  0.01233252,
4543       0.02135685,  0.01301289,  0.02161290,  0.01260280,  0.01867612,  0.01133078,
4544       0.04091313,  0.02668283,  0.03634781,  0.02733469,  0.03767692,  0.02836840,
4545       0.03906338,  0.02958073,  0.04025669,  0.02953292,  0.03665616,  0.02583565,
4546       0.06314558,  0.04830935,  0.05663609,  0.04564362,  0.05756200,  0.04739673,
4547       0.05967379,  0.04919083,  0.06124330,  0.04965808,  0.05879240,  0.04629319,
4548       0.08220739,  0.06924725,  0.07713556,  0.06782640,  0.07909499,  0.06992544,
4549       0.08116621,  0.07210181,  0.08281548,  0.07222669,  0.07941059,  0.06755612,
4550       0.01581588,  0.04533609,  0.02017939,  0.04342565,  0.02073232,  0.04476108,
4551       0.02117439,  0.04573358,  0.02129473,  0.04694267,  0.02220398,  0.05533458])
4552
4553       
4554        assert allclose(domain.quantities['ymomentum'].centroid_values,
4555     [-7.65882069e-005, -1.46087080e-004, -1.09630102e-004, -7.80950424e-005,
4556      -1.15922807e-005, -9.09134899e-005, -1.35994542e-004, -1.95673476e-004,
4557      -4.25779199e-004, -2.95890312e-004, -4.00060341e-004, -9.42021290e-005,
4558      -3.41372596e-004, -1.54560195e-004, -2.94810038e-004, -1.08844546e-004,
4559      -6.97240892e-005,  3.50299623e-005, -2.40159184e-004, -2.01805883e-004,
4560      -7.60732405e-004, -5.10897642e-004, -1.00940001e-003, -1.38037759e-004,
4561      -1.06169131e-003, -3.12307760e-004, -9.90602307e-004, -4.21634250e-005,
4562      -6.02424239e-004,  1.52230578e-004, -7.63833035e-004, -1.10273481e-004,
4563      -1.40187071e-003, -5.57831837e-004, -1.63988285e-003, -2.48018092e-004,
4564      -1.83309840e-003, -6.19360836e-004, -1.29955242e-003, -3.76237145e-004,
4565      -1.00613007e-003, -8.63641918e-005, -1.13604124e-003, -3.90589728e-004,
4566      -1.91457355e-003, -9.43783961e-004, -2.28090840e-003, -5.79107025e-004,
4567      -1.54091533e-003, -2.39785792e-003, -2.47947427e-003, -2.02694009e-003,
4568      -2.10441194e-003, -1.82082650e-003, -1.80229336e-003, -2.10418336e-003,
4569      -1.93104408e-003, -2.23200334e-003, -1.57239706e-003, -1.31486358e-003,
4570      -1.17564993e-003, -2.85846494e-003, -3.52956754e-003, -5.12658193e-003,
4571      -6.24238960e-003, -6.01820113e-003, -6.09602201e-003, -5.04787190e-003,
4572      -4.59373845e-003, -3.01393146e-003,  5.08550095e-004, -4.35896549e-004])
4573
4574        os.remove(domain.get_name() + '.sww')
4575
4576
4577    def test_temp_play(self):
4578
4579        from mesh_factory import rectangular
4580        from Numeric import array
4581
4582        #Create basic mesh
4583        points, vertices, boundary = rectangular(5, 5)
4584
4585        #Create shallow water domain
4586        domain = Domain(points, vertices, boundary)
4587        domain.smooth = False
4588        domain.default_order=2
4589        domain.beta_w      = 0.9
4590        domain.beta_w_dry  = 0.9
4591        domain.beta_uh     = 0.9
4592        domain.beta_uh_dry = 0.9
4593        domain.beta_vh     = 0.9
4594        domain.beta_vh_dry = 0.9
4595        domain.beta_h = 0.0 #Use first order in h-limiter
4596       
4597        # FIXME (Ole): Need tests where these two are commented out
4598        domain.H0 = 0        # Backwards compatibility (6/2/7)       
4599        domain.tight_slope_limiters = False # Backwards compatibility (14/4/7)
4600        domain.use_centroid_velocities = False # Backwards compatibility (7/5/8)
4601        domain.use_edge_limiter = False # Backwards compatibility (9/5/8)       
4602       
4603
4604        #Bed-slope and friction at vertices (and interpolated elsewhere)
4605        def x_slope(x, y):
4606            return -x/3
4607
4608        domain.set_quantity('elevation', x_slope)
4609
4610        # Boundary conditions
4611        Br = Reflective_boundary(domain)
4612        domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
4613
4614        #Initial condition
4615        domain.set_quantity('stage', expression='elevation+0.05')
4616        domain.check_integrity()
4617
4618        #Evolution
4619        for t in domain.evolve(yieldstep = 0.05, finaltime = 0.1):
4620            pass
4621
4622        assert allclose(domain.quantities['stage'].centroid_values[:4],
4623                        [0.00206836, 0.01296714, 0.00363415, 0.01438924])
4624        #print domain.quantities['xmomentum'].centroid_values[:4]
4625        assert allclose(domain.quantities['xmomentum'].centroid_values[:4],
4626                        [0.01360154, 0.00671133, 0.01264578, 0.00648503])
4627        assert allclose(domain.quantities['ymomentum'].centroid_values[:4],
4628                        [-1.19201077e-003, -7.23647546e-004,
4629                        -6.39083123e-005, 6.29815168e-005])
4630
4631        os.remove(domain.get_name() + '.sww')
4632
4633    def test_complex_bed(self):
4634        #No friction is tested here
4635
4636        from mesh_factory import rectangular
4637        from Numeric import array
4638
4639        N = 12
4640        points, vertices, boundary = rectangular(N, N/2, len1=1.2,len2=0.6,
4641                                                 origin=(-0.07, 0))
4642
4643
4644        domain = Domain(points, vertices, boundary)
4645        domain.smooth = False
4646        domain.default_order=2
4647
4648
4649        inflow_stage = 0.1
4650        Z = Weir(inflow_stage)
4651        domain.set_quantity('elevation', Z)
4652
4653        Br = Reflective_boundary(domain)
4654        Bd = Dirichlet_boundary([inflow_stage, 0.0, 0.0])
4655        domain.set_boundary({'left': Bd, 'right': Br, 'bottom': Br, 'top': Br})
4656
4657        domain.set_quantity('stage', expression='elevation')
4658
4659        for t in domain.evolve(yieldstep = 0.02, finaltime = 0.2):
4660            pass
4661
4662
4663        #print domain.quantities['stage'].centroid_values
4664
4665        #FIXME: These numbers were from version before 25/10
4666        #assert allclose(domain.quantities['stage'].centroid_values,
4667# [3.95822638e-002,  5.61022588e-002,  4.66437868e-002,  5.73081011e-002,
4668#  4.72394613e-002,  5.74684939e-002,  4.74309483e-002,  5.77458084e-002,
4669#  4.80628177e-002,  5.85656225e-002,  4.90498542e-002,  6.02609831e-002,
4670#  1.18470315e-002,  1.75136443e-002,  1.18035266e-002,  2.15565695e-002,
4671#  1.31620268e-002,  2.14351640e-002,  1.32351076e-002,  2.15450687e-002,
4672#  1.36414028e-002,  2.24274619e-002,  1.51689511e-002,  2.21789655e-002,
4673# -7.54337535e-003, -6.86362021e-004, -7.74146760e-003, -1.83756530e-003,
4674# -8.16773628e-003, -4.49916813e-004, -8.08202599e-003, -3.91118720e-004,
4675# -8.10292716e-003, -3.88584984e-004, -7.35226124e-003,  2.73985295e-004,
4676#  1.86166683e-001,  8.74070369e-002,  1.86166712e-001,  8.74035875e-002,
4677#  6.11666935e-002, -3.76173225e-002, -6.38333276e-002, -3.76147365e-002,
4678#  6.11666725e-002,  8.73846774e-002,  1.86166697e-001,  8.74171550e-002,
4679# -4.83333333e-002,  1.18333333e-001, -4.83333333e-002,  1.18333333e-001,
4680# -4.83333333e-002, -6.66666667e-003, -1.73333333e-001, -1.31666667e-001,
4681# -1.73333333e-001, -6.66666667e-003, -4.83333333e-002,  1.18333333e-001,
4682# -2.48333333e-001, -2.31666667e-001, -2.48333333e-001, -2.31666667e-001,
4683# -2.48333333e-001, -2.31666667e-001, -2.48333333e-001, -2.31666667e-001,
4684# -2.48333333e-001, -2.31666667e-001, -2.48333333e-001, -2.31666667e-001,
4685# -4.65000000e-001, -3.65000000e-001, -4.65000000e-001, -3.65000000e-001,
4686# -4.65000000e-001, -3.65000000e-001, -4.65000000e-001, -3.65000000e-001,
4687# -4.65000000e-001, -3.65000000e-001, -4.65000000e-001, -3.65000000e-001,
4688# -5.98333333e-001, -5.81666667e-001, -5.98333333e-001, -5.81666667e-001,
4689# -5.98333333e-001, -5.81666667e-001, -5.98333333e-001, -5.81666667e-001,
4690# -5.98333333e-001, -5.81666667e-001, -5.98333333e-001, -5.81666667e-001,
4691# -6.48333333e-001, -6.31666667e-001, -6.48333333e-001, -6.31666667e-001,
4692# -6.48333333e-001, -6.31666667e-001, -6.48333333e-001, -6.31666667e-001,
4693# -6.48333333e-001, -6.31666667e-001, -6.48333333e-001, -6.31666667e-001,
4694# -5.31666667e-001, -5.98333333e-001, -5.31666667e-001, -5.98333333e-001,
4695# -5.31666667e-001, -5.98333333e-001, -5.31666667e-001, -5.98333333e-001,
4696# -5.31666667e-001, -5.98333333e-001, -5.31666667e-001, -5.98333333e-001,
4697# -4.98333333e-001, -4.81666667e-001, -4.98333333e-001, -4.81666667e-001,
4698# -4.98333333e-001, -4.81666667e-001, -4.98333333e-001, -4.81666667e-001,
4699# -4.98333333e-001, -4.81666667e-001, -4.98333333e-001, -4.81666667e-001,
4700# -5.48333333e-001, -5.31666667e-001, -5.48333333e-001, -5.31666667e-001,
4701# -5.48333333e-001, -5.31666667e-001, -5.48333333e-001, -5.31666667e-001,
4702# -5.48333333e-001, -5.31666667e-001, -5.48333333e-001, -5.31666667e-001])
4703
4704        os.remove(domain.get_name() + '.sww')
4705
4706    def test_spatio_temporal_boundary_1(self):
4707        """Test that boundary values can be read from file and interpolated
4708        in both time and space.
4709
4710        Verify that the same steady state solution is arrived at and that
4711        time interpolation works.
4712
4713        The full solution history is not exactly the same as
4714        file boundary must read and interpolate from *smoothed* version
4715        as stored in sww.
4716        """
4717        import time
4718
4719        #Create sww file of simple propagation from left to right
4720        #through rectangular domain
4721
4722        from mesh_factory import rectangular
4723
4724        #Create basic mesh
4725        points, vertices, boundary = rectangular(3, 3)
4726
4727        #Create shallow water domain
4728        domain1 = Domain(points, vertices, boundary)
4729
4730        domain1.reduction = mean
4731        domain1.smooth = False  #Exact result
4732
4733        domain1.default_order = 2
4734        domain1.store = True
4735        domain1.set_datadir('.')
4736        domain1.set_name('spatio_temporal_boundary_source' + str(time.time()))
4737
4738        #FIXME: This is extremely important!
4739        #How can we test if they weren't stored?
4740        domain1.quantities_to_be_stored = ['stage', 'xmomentum', 'ymomentum']
4741
4742
4743        #Bed-slope and friction at vertices (and interpolated elsewhere)
4744        domain1.set_quantity('elevation', 0)
4745        domain1.set_quantity('friction', 0)
4746
4747        # Boundary conditions
4748        Br = Reflective_boundary(domain1)
4749        Bd = Dirichlet_boundary([0.3,0,0])
4750        domain1.set_boundary({'left': Bd, 'top': Bd, 'right': Br, 'bottom': Br})
4751        #Initial condition
4752        domain1.set_quantity('stage', 0)
4753        domain1.check_integrity()
4754
4755        finaltime = 5
4756        #Evolution  (full domain - large steps)
4757        for t in domain1.evolve(yieldstep = 0.671, finaltime = finaltime):
4758            pass
4759            #domain1.write_time()
4760
4761        cv1 = domain1.quantities['stage'].centroid_values
4762
4763
4764        #Create a triangle shaped domain (reusing coordinates from domain 1),
4765        #formed from the lower and right hand  boundaries and
4766        #the sw-ne diagonal
4767        #from domain 1. Call it domain2
4768
4769        points = [ [0,0], [1.0/3,0], [1.0/3,1.0/3],
4770                   [2.0/3,0], [2.0/3,1.0/3], [2.0/3,2.0/3],
4771                   [1,0], [1,1.0/3], [1,2.0/3], [1,1]]
4772
4773        vertices = [ [1,2,0], [3,4,1], [2,1,4], [4,5,2],
4774                     [6,7,3], [4,3,7], [7,8,4], [5,4,8], [8,9,5]]
4775
4776        boundary = { (0,1):'bottom', (1,1):'bottom', (4,1): 'bottom',
4777                     (4,2):'right', (6,2):'right', (8,2):'right',
4778                     (0,0):'diagonal', (3,0):'diagonal', (8,0):'diagonal'}
4779
4780        domain2 = Domain(points, vertices, boundary)
4781
4782        domain2.reduction = domain1.reduction
4783        domain2.smooth = False
4784        domain2.default_order = 2
4785
4786        #Bed-slope and friction at vertices (and interpolated elsewhere)
4787        domain2.set_quantity('elevation', 0)
4788        domain2.set_quantity('friction', 0)
4789        domain2.set_quantity('stage', 0)
4790
4791        # Boundary conditions
4792        Br = Reflective_boundary(domain2)
4793        #Bf = Spatio_temporal_boundary(domain1.get_name() + '.' +\
4794        #                              domain1.format, domain2)
4795        Bf = Field_boundary(domain1.get_name() + '.' +\
4796                            domain1.format, domain2)       
4797        domain2.set_boundary({'right':Br, 'bottom':Br, 'diagonal':Bf})
4798        domain2.check_integrity()
4799
4800
4801
4802        #Evolution (small steps)
4803        for t in domain2.evolve(yieldstep = 0.0711, finaltime = finaltime):
4804            pass
4805
4806
4807        #Use output from domain1 as spatio-temporal boundary for domain2
4808        #and verify that results at right hand side are close.
4809
4810        cv2 = domain2.quantities['stage'].centroid_values
4811
4812        #print take(cv1, (12,14,16))  #Right
4813        #print take(cv2, (4,6,8))
4814        #print take(cv1, (0,6,12))  #Bottom
4815        #print take(cv2, (0,1,4))
4816        #print take(cv1, (0,8,16))  #Diag
4817        #print take(cv2, (0,3,8))
4818
4819        assert allclose( take(cv1, (0,8,16)), take(cv2, (0,3,8))) #Diag
4820        assert allclose( take(cv1, (0,6,12)), take(cv2, (0,1,4))) #Bottom
4821        assert allclose( take(cv1, (12,14,16)), take(cv2, (4,6,8))) #RHS
4822
4823        #Cleanup
4824        os.remove(domain1.get_name() + '.' + domain1.format)
4825        os.remove(domain2.get_name() + '.' + domain2.format)       
4826
4827
4828
4829    def test_spatio_temporal_boundary_2(self):
4830        """Test that boundary values can be read from file and interpolated
4831        in both time and space.
4832        This is a more basic test, verifying that boundary object
4833        produces the expected results
4834
4835
4836        """
4837        import time
4838
4839        #Create sww file of simple propagation from left to right
4840        #through rectangular domain
4841
4842        from mesh_factory import rectangular
4843
4844        #Create basic mesh
4845        points, vertices, boundary = rectangular(3, 3)
4846
4847        #Create shallow water domain
4848        domain1 = Domain(points, vertices, boundary)
4849
4850        domain1.reduction = mean
4851        domain1.smooth = True #To mimic MOST output
4852
4853        domain1.default_order = 2
4854        domain1.store = True
4855        domain1.set_datadir('.')
4856        domain1.set_name('spatio_temporal_boundary_source' + str(time.time()))
4857
4858        #FIXME: This is extremely important!
4859        #How can we test if they weren't stored?
4860        domain1.quantities_to_be_stored = ['stage', 'xmomentum', 'ymomentum']
4861
4862
4863        #Bed-slope and friction at vertices (and interpolated elsewhere)
4864        domain1.set_quantity('elevation', 0)
4865        domain1.set_quantity('friction', 0)
4866
4867        # Boundary conditions
4868        Br = Reflective_boundary(domain1)
4869        Bd = Dirichlet_boundary([0.3,0,0])
4870        domain1.set_boundary({'left': Bd, 'top': Bd, 'right': Br, 'bottom': Br})
4871        #Initial condition
4872        domain1.set_quantity('stage', 0)
4873        domain1.check_integrity()
4874
4875        finaltime = 5
4876        #Evolution  (full domain - large steps)
4877        for t in domain1.evolve(yieldstep = 1, finaltime = finaltime):
4878            pass
4879            #domain1.write_time()
4880
4881
4882        #Create an triangle shaped domain (coinciding with some
4883        #coordinates from domain 1),
4884        #formed from the lower and right hand  boundaries and
4885        #the sw-ne diagonal
4886        #from domain 1. Call it domain2
4887
4888        points = [ [0,0], [1.0/3,0], [1.0/3,1.0/3],
4889                   [2.0/3,0], [2.0/3,1.0/3], [2.0/3,2.0/3],
4890                   [1,0], [1,1.0/3], [1,2.0/3], [1,1]]
4891
4892        vertices = [ [1,2,0],
4893                     [3,4,1], [2,1,4], [4,5,2],
4894                     [6,7,3], [4,3,7], [7,8,4], [5,4,8], [8,9,5]]
4895
4896        boundary = { (0,1):'bottom', (1,1):'bottom', (4,1): 'bottom',
4897                     (4,2):'right', (6,2):'right', (8,2):'right',
4898                     (0,0):'diagonal', (3,0):'diagonal', (8,0):'diagonal'}
4899
4900        domain2 = Domain(points, vertices, boundary)
4901
4902        domain2.reduction = domain1.reduction
4903        domain2.smooth = False
4904        domain2.default_order = 2
4905
4906        #Bed-slope and friction at vertices (and interpolated elsewhere)
4907        domain2.set_quantity('elevation', 0)
4908        domain2.set_quantity('friction', 0)
4909        domain2.set_quantity('stage', 0)
4910
4911
4912        #Read results for specific timesteps t=1 and t=2
4913        from Scientific.IO.NetCDF import NetCDFFile
4914        fid = NetCDFFile(domain1.get_name() + '.' + domain1.format)
4915
4916        x = fid.variables['x'][:]
4917        y = fid.variables['y'][:]
4918        s1 = fid.variables['stage'][1,:]
4919        s2 = fid.variables['stage'][2,:]
4920        fid.close()
4921
4922        from Numeric import take, reshape, concatenate
4923        shp = (len(x), 1)
4924        points = concatenate( (reshape(x, shp), reshape(y, shp)), axis=1)
4925        #The diagonal points of domain 1 are 0, 5, 10, 15
4926
4927        #print points[0], points[5], points[10], points[15]
4928        assert allclose( take(points, [0,5,10,15]),
4929                         [[0,0], [1.0/3, 1.0/3], [2.0/3, 2.0/3], [1,1]])
4930
4931
4932        # Boundary conditions
4933        Br = Reflective_boundary(domain2)
4934        #Bf = Spatio_temporal_boundary(domain1.get_name() + '.' + domain1.format,
4935        #                              domain2)
4936        Bf = Field_boundary(domain1.get_name() + '.' + domain1.format,
4937                            domain2, verbose=False)       
4938        domain2.set_boundary({'right':Br, 'bottom':Br, 'diagonal':Bf})
4939        domain2.check_integrity()
4940
4941        #Test that interpolation points are the mid points of the all boundary
4942        #segments
4943
4944        boundary_midpoints = [[1.0/6, 0], [1.0/2, 0], [5.0/6,0],
4945                              [1.0, 1.0/6], [1.0, 1.0/2], [1.0, 5.0/6],
4946                              [1.0/6, 1.0/6], [0.5, 0.5], [5.0/6, 5.0/6]]
4947
4948        boundary_midpoints.sort()
4949        R = Bf.F.interpolation_points.tolist()
4950        R.sort()
4951        assert allclose(boundary_midpoints, R)
4952
4953        #Check spatially interpolated output at time == 1
4954        domain2.time = 1
4955
4956        #First diagonal midpoint
4957        R0 = Bf.evaluate(0,0)
4958        assert allclose(R0[0], (s1[0] + s1[5])/2)
4959
4960        #Second diagonal midpoint
4961        R0 = Bf.evaluate(3,0)
4962        assert allclose(R0[0], (s1[5] + s1[10])/2)
4963
4964        #First diagonal midpoint
4965        R0 = Bf.evaluate(8,0)
4966        assert allclose(R0[0], (s1[10] + s1[15])/2)
4967
4968        #Check spatially interpolated output at time == 2
4969        domain2.time = 2
4970
4971        #First diagonal midpoint
4972        R0 = Bf.evaluate(0,0)
4973        assert allclose(R0[0], (s2[0] + s2[5])/2)
4974
4975        #Second diagonal midpoint
4976        R0 = Bf.evaluate(3,0)
4977        assert allclose(R0[0], (s2[5] + s2[10])/2)
4978
4979        #First diagonal midpoint
4980        R0 = Bf.evaluate(8,0)
4981        assert allclose(R0[0], (s2[10] + s2[15])/2)
4982
4983
4984        #Now check temporal interpolation
4985
4986        domain2.time = 1 + 2.0/3
4987
4988        #First diagonal midpoint
4989        R0 = Bf.evaluate(0,0)
4990        assert allclose(R0[0], ((s1[0] + s1[5])/2 + 2.0*(s2[0] + s2[5])/2)/3)
4991
4992        #Second diagonal midpoint
4993        R0 = Bf.evaluate(3,0)
4994        assert allclose(R0[0], ((s1[5] + s1[10])/2 + 2.0*(s2[5] + s2[10])/2)/3)
4995
4996        #First diagonal midpoint
4997        R0 = Bf.evaluate(8,0)
4998        assert allclose(R0[0], ((s1[10] + s1[15])/2 + 2.0*(s2[10] + s2[15])/2)/3)
4999
5000
5001
5002        #Cleanup
5003        os.remove(domain1.get_name() + '.' + domain1.format)
5004
5005
5006    def test_spatio_temporal_boundary_3(self):
5007        """Test that boundary values can be read from file and interpolated
5008        in both time and space.
5009        This is a more basic test, verifying that boundary object
5010        produces the expected results
5011
5012        This tests adjusting using mean_stage
5013
5014        """
5015
5016        import time
5017
5018        mean_stage = 5.2 # Adjust stage by this amount in boundary
5019
5020        #Create sww file of simple propagation from left to right
5021        #through rectangular domain
5022
5023        from mesh_factory import rectangular
5024
5025        #Create basic mesh
5026        points, vertices, boundary = rectangular(3, 3)
5027
5028        #Create shallow water domain
5029        domain1 = Domain(points, vertices, boundary)
5030
5031        domain1.reduction = mean
5032        domain1.smooth = True #To mimic MOST output
5033
5034        domain1.default_order = 2
5035        domain1.store = True
5036        domain1.set_datadir('.')
5037        domain1.set_name('spatio_temporal_boundary_source' + str(time.time()))
5038
5039        #FIXME: This is extremely important!
5040        #How can we test if they weren't stored?
5041        domain1.quantities_to_be_stored = ['stage', 'xmomentum', 'ymomentum']
5042
5043
5044        #Bed-slope and friction at vertices (and interpolated elsewhere)
5045        domain1.set_quantity('elevation', 0)
5046        domain1.set_quantity('friction', 0)
5047
5048        # Boundary conditions
5049        Br = Reflective_boundary(domain1)
5050        Bd = Dirichlet_boundary([0.3,0,0])
5051        domain1.set_boundary({'left': Bd, 'top': Bd, 'right': Br, 'bottom': Br})
5052        #Initial condition
5053        domain1.set_quantity('stage', 0)
5054        domain1.check_integrity()
5055
5056        finaltime = 5
5057        #Evolution  (full domain - large steps)
5058        for t in domain1.evolve(yieldstep = 1, finaltime = finaltime):
5059            pass
5060            #domain1.write_time()
5061
5062
5063        #Create an triangle shaped domain (coinciding with some
5064        #coordinates from domain 1),
5065        #formed from the lower and right hand  boundaries and
5066        #the sw-ne diagonal
5067        #from domain 1. Call it domain2
5068
5069        points = [ [0,0],
5070                   [1.0/3,0], [1.0/3,1.0/3],
5071                   [2.0/3,0], [2.0/3,1.0/3], [2.0/3,2.0/3],
5072                   [1,0],     [1,1.0/3],     [1,2.0/3],     [1,1]] 
5073                   
5074        vertices = [ [1,2,0],
5075                     [3,4,1], [2,1,4], [4,5,2],
5076                     [6,7,3], [4,3,7], [7,8,4], [5,4,8], [8,9,5]]
5077
5078        boundary = { (0,1):'bottom', (1,1):'bottom', (4,1): 'bottom',
5079                     (4,2):'right', (6,2):'right', (8,2):'right',
5080                     (0,0):'diagonal', (3,0):'diagonal', (8,0):'diagonal'}
5081
5082        domain2 = Domain(points, vertices, boundary)
5083
5084        domain2.reduction = domain1.reduction
5085        domain2.smooth = False
5086        domain2.default_order = 2
5087
5088        #Bed-slope and friction at vertices (and interpolated elsewhere)
5089        domain2.set_quantity('elevation', 0)
5090        domain2.set_quantity('friction', 0)
5091        domain2.set_quantity('stage', 0)
5092
5093
5094        #Read results for specific timesteps t=1 and t=2
5095        from Scientific.IO.NetCDF import NetCDFFile
5096        fid = NetCDFFile(domain1.get_name() + '.' + domain1.format)
5097
5098        x = fid.variables['x'][:]
5099        y = fid.variables['y'][:]
5100        s1 = fid.variables['stage'][1,:]
5101        s2 = fid.variables['stage'][2,:]
5102        fid.close()
5103
5104        from Numeric import take, reshape, concatenate
5105        shp = (len(x), 1)
5106        points = concatenate( (reshape(x, shp), reshape(y, shp)), axis=1)
5107        #The diagonal points of domain 1 are 0, 5, 10, 15
5108
5109        #print points[0], points[5], points[10], points[15]
5110        assert allclose( take(points, [0,5,10,15]),
5111                         [[0,0], [1.0/3, 1.0/3], [2.0/3, 2.0/3], [1,1]])
5112
5113
5114        # Boundary conditions
5115        Br = Reflective_boundary(domain2)
5116        #Bf = Spatio_temporal_boundary(domain1.get_name() + '.' + domain1.format,
5117        #                              domain2)
5118        Bf = Field_boundary(domain1.get_name() + '.' + domain1.format,
5119                            domain2, mean_stage=mean_stage, verbose=False)
5120       
5121        domain2.set_boundary({'right':Br, 'bottom':Br, 'diagonal':Bf})
5122        domain2.check_integrity()
5123
5124        #Test that interpolation points are the mid points of the all boundary
5125        #segments
5126
5127        boundary_midpoints = [[1.0/6, 0], [1.0/2, 0], [5.0/6,0],
5128                              [1.0, 1.0/6], [1.0, 1.0/2], [1.0, 5.0/6],
5129                              [1.0/6, 1.0/6], [0.5, 0.5], [5.0/6, 5.0/6]]
5130
5131        boundary_midpoints.sort()
5132        R = Bf.F.interpolation_points.tolist()
5133        R.sort()
5134        assert allclose(boundary_midpoints, R)
5135
5136        #Check spatially interpolated output at time == 1
5137        domain2.time = 1
5138
5139        #First diagonal midpoint
5140        R0 = Bf.evaluate(0,0)
5141        assert allclose(R0[0], (s1[0] + s1[5])/2 + mean_stage)
5142
5143        #Second diagonal midpoint
5144        R0 = Bf.evaluate(3,0)
5145        assert allclose(R0[0], (s1[5] + s1[10])/2 + mean_stage)
5146
5147        #First diagonal midpoint
5148        R0 = Bf.evaluate(8,0)
5149        assert allclose(R0[0], (s1[10] + s1[15])/2 + mean_stage)
5150
5151        #Check spatially interpolated output at time == 2
5152        domain2.time = 2
5153
5154        #First diagonal midpoint
5155        R0 = Bf.evaluate(0,0)
5156        assert allclose(R0[0], (s2[0] + s2[5])/2 + mean_stage)
5157
5158        #Second diagonal midpoint
5159        R0 = Bf.evaluate(3,0)
5160        assert allclose(R0[0], (s2[5] + s2[10])/2 + mean_stage)
5161
5162        #First diagonal midpoint
5163        R0 = Bf.evaluate(8,0)
5164        assert allclose(R0[0], (s2[10] + s2[15])/2 + mean_stage)
5165
5166
5167        #Now check temporal interpolation
5168
5169        domain2.time = 1 + 2.0/3
5170
5171        #First diagonal midpoint
5172        R0 = Bf.evaluate(0,0)
5173        assert allclose(R0[0], ((s1[0] + s1[5])/2 + 2.0*(s2[0] + s2[5])/2)/3 + mean_stage)
5174
5175        #Second diagonal midpoint
5176        R0 = Bf.evaluate(3,0)
5177        assert allclose(R0[0], ((s1[5] + s1[10])/2 + 2.0*(s2[5] + s2[10])/2)/3 + mean_stage)
5178
5179        #First diagonal midpoint
5180        R0 = Bf.evaluate(8,0)
5181        assert allclose(R0[0], ((s1[10] + s1[15])/2 + 2.0*(s2[10] + s2[15])/2)/3 + mean_stage)
5182
5183
5184        #Cleanup
5185        os.remove(domain1.get_name() + '.' + domain1.format)
5186
5187
5188    def test_spatio_temporal_boundary_outside(self):
5189        """Test that field_boundary catches if a point is outside the sww that defines it
5190        """
5191
5192        import time
5193        #Create sww file of simple propagation from left to right
5194        #through rectangular domain
5195
5196        from mesh_factory import rectangular
5197
5198        #Create basic mesh
5199        points, vertices, boundary = rectangular(3, 3)
5200
5201        #Create shallow water domain
5202        domain1 = Domain(points, vertices, boundary)
5203
5204        domain1.reduction = mean
5205        domain1.smooth = True #To mimic MOST output
5206
5207        domain1.default_order = 2
5208        domain1.store = True
5209        domain1.set_datadir('.')
5210        domain1.set_name('spatio_temporal_boundary_source' + str(time.time()))
5211
5212        #FIXME: This is extremely important!
5213        #How can we test if they weren't stored?
5214        domain1.quantities_to_be_stored = ['stage', 'xmomentum', 'ymomentum']
5215
5216
5217        #Bed-slope and friction at vertices (and interpolated elsewhere)
5218        domain1.set_quantity('elevation', 0)
5219        domain1.set_quantity('friction', 0)
5220
5221        # Boundary conditions
5222        Br = Reflective_boundary(domain1)
5223        Bd = Dirichlet_boundary([0.3,0,0])
5224        domain1.set_boundary({'left': Bd, 'top': Bd, 'right': Br, 'bottom': Br})
5225        #Initial condition
5226        domain1.set_quantity('stage', 0)
5227        domain1.check_integrity()
5228
5229        finaltime = 5
5230        #Evolution  (full domain - large steps)
5231        for t in domain1.evolve(yieldstep = 1, finaltime = finaltime):
5232            pass
5233            #domain1.write_time()
5234
5235
5236        #Create an triangle shaped domain (coinciding with some
5237        #coordinates from domain 1, but one edge outside!),
5238        #formed from the lower and right hand  boundaries and
5239        #the sw-ne diagonal as in the previous test but scaled
5240        #in the x direction by a factor of 2
5241
5242        points = [ [0,0],
5243                   [2.0/3,0], [2.0/3,1.0/3],
5244                   [4.0/3,0], [4.0/3,1.0/3], [4.0/3,2.0/3],
5245                   [2,0],     [2,1.0/3],     [2,2.0/3],     [2,1] 
5246                   ]
5247
5248        vertices = [ [1,2,0],
5249                     [3,4,1], [2,1,4], [4,5,2],
5250                     [6,7,3], [4,3,7], [7,8,4], [5,4,8], [8,9,5]]
5251
5252        boundary = { (0,1):'bottom', (1,1):'bottom', (4,1): 'bottom',
5253                     (4,2):'right', (6,2):'right', (8,2):'right',
5254                     (0,0):'diagonal', (3,0):'diagonal', (8,0):'diagonal'}
5255
5256        domain2 = Domain(points, vertices, boundary)
5257
5258        domain2.reduction = domain1.reduction
5259        domain2.smooth = False
5260        domain2.default_order = 2
5261
5262        #Bed-slope and friction at vertices (and interpolated elsewhere)
5263        domain2.set_quantity('elevation', 0)
5264        domain2.set_quantity('friction', 0)
5265        domain2.set_quantity('stage', 0)
5266
5267
5268        #Read results for specific timesteps t=1 and t=2
5269        from Scientific.IO.NetCDF import NetCDFFile
5270        fid = NetCDFFile(domain1.get_name() + '.' + domain1.format)
5271
5272        x = fid.variables['x'][:]
5273        y = fid.variables['y'][:]
5274        s1 = fid.variables['stage'][1,:]
5275        s2 = fid.variables['stage'][2,:]
5276        fid.close()
5277
5278        from Numeric import take, reshape, concatenate
5279        shp = (len(x), 1)
5280        points = concatenate( (reshape(x, shp), reshape(y, shp)), axis=1)
5281        #The diagonal points of domain 1 are 0, 5, 10, 15
5282
5283        #print points[0], points[5], points[10], points[15]
5284        assert allclose( take(points, [0,5,10,15]),
5285                         [[0,0], [1.0/3, 1.0/3], [2.0/3, 2.0/3], [1,1]])
5286
5287
5288        # Boundary conditions
5289        Br = Reflective_boundary(domain2)
5290        #Bf = Spatio_temporal_boundary(domain1.get_name() + '.' + domain1.format,
5291        #                              domain2)
5292        Bf = Field_boundary(domain1.get_name() + '.' + domain1.format,
5293                            domain2, mean_stage=1, verbose=False)
5294       
5295        domain2.set_boundary({'right':Br, 'bottom':Br, 'diagonal':Bf})
5296        domain2.check_integrity()
5297
5298        try:
5299            for t in domain2.evolve(yieldstep = 1, finaltime = finaltime):
5300                pass
5301        except:
5302            pass
5303        else:
5304            msg = 'This should have caught NAN at boundary'
5305            raise Exception, msg
5306
5307
5308        #Cleanup
5309        os.remove(domain1.get_name() + '.' + domain1.format)
5310
5311
5312
5313
5314    def test_extrema(self):
5315        """Test that extrema of quantities are computed correctly
5316        Extrema are updated at every *internal* timestep
5317        """
5318
5319        from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular_cross
5320
5321        initial_runup_height = -0.4
5322        final_runup_height = -0.3
5323
5324
5325        #--------------------------------------------------------------
5326        # Setup computational domain
5327        #--------------------------------------------------------------
5328        N = 5
5329        points, vertices, boundary = rectangular_cross(N, N) 
5330        domain = Domain(points, vertices, boundary)
5331        domain.set_name('extrema_test')
5332
5333        #--------------------------------------------------------------
5334        # Setup initial conditions
5335        #--------------------------------------------------------------
5336        def topography(x,y):
5337            return -x/2                             # linear bed slope
5338           
5339
5340        domain.set_quantity('elevation', topography)       # Use function for elevation
5341        domain.set_quantity('friction', 0.)                # Zero friction
5342        domain.set_quantity('stage', initial_runup_height) # Constant negative initial stage
5343        domain.set_quantities_to_be_monitored(['stage', 'stage-elevation'],
5344                                              time_interval = [0.5, 2.7],
5345                                              polygon = [[0,0], [0,1], [1,1], [1,0]])
5346       
5347        assert len(domain.quantities_to_be_monitored) == 2
5348        assert domain.quantities_to_be_monitored.has_key('stage')
5349        assert domain.quantities_to_be_monitored.has_key('stage-elevation')
5350        for key in domain.quantities_to_be_monitored['stage'].keys():
5351            assert domain.quantities_to_be_monitored['stage'][key] is None       
5352
5353
5354        #--------------------------------------------------------------
5355        # Setup boundary conditions
5356        #--------------------------------------------------------------
5357        Br = Reflective_boundary(domain)              # Reflective wall
5358        Bd = Dirichlet_boundary([final_runup_height,  # Constant inflow
5359                                 0,
5360                                 0])
5361
5362        # All reflective to begin with (still water)
5363        domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
5364
5365
5366        #--------------------------------------------------------------
5367        # Let triangles adjust and check extrema
5368        #--------------------------------------------------------------
5369        for t in domain.evolve(yieldstep = 0.1, finaltime = 1.0):
5370            domain.quantity_statistics() # Run it silently
5371
5372
5373
5374        #--------------------------------------------------------------
5375        # Test extrema
5376        #--------------------------------------------------------------
5377
5378        stage = domain.quantities_to_be_monitored['stage']
5379        assert stage['min'] <= stage['max']
5380
5381        #print stage['min'], stage['max']
5382        assert allclose(stage['min'], initial_runup_height,
5383                        rtol = 1.0/N) # First order accuracy
5384
5385
5386        depth = domain.quantities_to_be_monitored['stage-elevation']
5387        assert depth['min'] <= depth['max'] 
5388        assert depth['min'] >= 0.0
5389        assert depth['max'] >= 0.0       
5390        ##assert depth[1] <= ?? initial_runup_height       
5391
5392
5393        #--------------------------------------------------------------
5394        # Update boundary to allow inflow
5395        #--------------------------------------------------------------
5396        domain.set_boundary({'right': Bd})
5397
5398       
5399        #--------------------------------------------------------------
5400        # Evolve system through time
5401        #--------------------------------------------------------------
5402        for t in domain.evolve(yieldstep = 0.1, finaltime = 3.0):
5403            #domain.write_time()
5404            domain.quantity_statistics() # Run it silently           
5405           
5406   
5407        #--------------------------------------------------------------
5408        # Test extrema again
5409        #--------------------------------------------------------------
5410
5411        stage = domain.quantities_to_be_monitored['stage']
5412        assert stage['min'] <= stage['max']
5413
5414        assert allclose(stage['min'], initial_runup_height,
5415                        rtol = 1.0/N) # First order accuracy       
5416
5417        depth = domain.quantities_to_be_monitored['stage-elevation']
5418        assert depth['min'] <= depth['max'] 
5419        assert depth['min'] >= 0.0
5420        assert depth['max'] >= 0.0       
5421
5422        #Cleanup
5423        os.remove(domain.get_name() + '.' + domain.format)
5424       
5425
5426
5427    def test_tight_slope_limiters(self):
5428        """Test that new slope limiters (Feb 2007) don't induce extremely
5429        small timesteps. This test actually reveals the problem as it
5430        was in March-April 2007
5431        """
5432
5433        import time, os
5434        from Numeric import array, zeros, allclose, Float, concatenate
5435        from Scientific.IO.NetCDF import NetCDFFile
5436        from data_manager import get_dataobject, extent_sww
5437        from mesh_factory import rectangular
5438
5439       
5440        #Create basic mesh
5441        points, vertices, boundary = rectangular(2, 2)
5442
5443        #Create shallow water domain
5444        domain = Domain(points, vertices, boundary)
5445        domain.default_order = 2
5446
5447        # This will pass
5448        #domain.tight_slope_limiters = 1
5449        #domain.H0 = 0.01
5450       
5451        # This will fail
5452        #domain.tight_slope_limiters = 1
5453        #domain.H0 = 0.001
5454
5455        # This will pass provided C extension implements limiting of
5456        # momentum in _compute_speeds
5457        domain.tight_slope_limiters = 1
5458        domain.H0 = 0.001       
5459        domain.protect_against_isolated_degenerate_timesteps = True       
5460
5461        #Set some field values
5462        domain.set_quantity('elevation', lambda x,y: -x)
5463        domain.set_quantity('friction', 0.03)
5464
5465
5466        ######################
5467        # Boundary conditions
5468        B = Transmissive_boundary(domain)
5469        domain.set_boundary( {'left': B, 'right': B, 'top': B, 'bottom': B})
5470
5471
5472        ######################
5473        #Initial condition - with jumps
5474
5475
5476        bed = domain.quantities['elevation'].vertex_values
5477        stage = zeros(bed.shape, Float)
5478
5479        h = 0.3
5480        for i in range(stage.shape[0]):
5481            if i % 2 == 0:
5482                stage[i,:] = bed[i,:] + h
5483            else:
5484                stage[i,:] = bed[i,:]
5485
5486        domain.set_quantity('stage', stage)
5487
5488
5489        domain.distribute_to_vertices_and_edges()               
5490
5491       
5492
5493        domain.set_name('tight_limiters')
5494        domain.format = 'sww'
5495        domain.smooth = True
5496        domain.reduction = mean
5497        domain.set_datadir('.')
5498        domain.smooth = False
5499        domain.store = True
5500        domain.beta_h = 0
5501       
5502
5503        #Evolution
5504        for t in domain.evolve(yieldstep = 0.1, finaltime = 0.3):
5505           
5506            #domain.write_time(track_speeds=True)
5507            stage = domain.quantities['stage'].vertex_values
5508
5509            #Get NetCDF
5510            fid = NetCDFFile(domain.writer.filename, 'r')
5511            stage_file = fid.variables['stage']
5512           
5513            fid.close()
5514
5515        os.remove(domain.writer.filename)
5516
5517
5518    def test_pmesh2Domain(self):
5519         import os
5520         import tempfile
5521
5522         fileName = tempfile.mktemp(".tsh")
5523         file = open(fileName,"w")
5524         file.write("4 3 # <vertex #> <x> <y> [attributes]\n \
55250 0.0 0.0 0.0 0.0 0.01 \n \
55261 1.0 0.0 10.0 10.0 0.02  \n \
55272 0.0 1.0 0.0 10.0 0.03  \n \
55283 0.5 0.25 8.0 12.0 0.04  \n \
5529# Vert att title  \n \
5530elevation  \n \
5531stage  \n \
5532friction  \n \
55332 # <triangle #> [<vertex #>] [<neigbouring triangle #>]  \n\
55340 0 3 2 -1  -1  1 dsg\n\
55351 0 1 3 -1  0 -1   ole nielsen\n\
55364 # <segment #> <vertex #>  <vertex #> [boundary tag] \n\
55370 1 0 2 \n\
55381 0 2 3 \n\
55392 2 3 \n\
55403 3 1 1 \n\
55413 0 # <x> <y> [attributes] ...Mesh Vertices... \n \
55420 216.0 -86.0 \n \
55431 160.0 -167.0 \n \
55442 114.0 -91.0 \n \
55453 # <vertex #>  <vertex #> [boundary tag] ...Mesh Segments... \n \
55460 0 1 0 \n \
55471 1 2 0 \n \
55482 2 0 0 \n \
55490 # <x> <y> ...Mesh Holes... \n \
55500 # <x> <y> <attribute>...Mesh Regions... \n \
55510 # <x> <y> <attribute>...Mesh Regions, area... \n\
5552#Geo reference \n \
555356 \n \
5554140 \n \
5555120 \n")
5556         file.close()
5557
5558         tags = {}
5559         b1 =  Dirichlet_boundary(conserved_quantities = array([0.0]))
5560         b2 =  Dirichlet_boundary(conserved_quantities = array([1.0]))
5561         b3 =  Dirichlet_boundary(conserved_quantities = array([2.0]))
5562         tags["1"] = b1
5563         tags["2"] = b2
5564         tags["3"] = b3
5565
5566         #from anuga.abstract_2d_finite_volumes.pmesh2domain import pmesh_to_domain_instance
5567         #domain = pmesh_to_domain_instance(fileName, Domain)
5568
5569         domain = Domain(mesh_filename=fileName)
5570                         #verbose=True, use_cache=True)
5571         
5572         #print "domain.tagged_elements", domain.tagged_elements
5573         ## check the quantities
5574         #print domain.quantities['elevation'].vertex_values
5575         answer = [[0., 8., 0.],
5576                   [0., 10., 8.]]
5577         assert allclose(domain.quantities['elevation'].vertex_values,
5578                        answer)
5579
5580         #print domain.quantities['stage'].vertex_values
5581         answer = [[0., 12., 10.],
5582                   [0., 10., 12.]]
5583         assert allclose(domain.quantities['stage'].vertex_values,
5584                        answer)
5585
5586         #print domain.quantities['friction'].vertex_values
5587         answer = [[0.01, 0.04, 0.03],
5588                   [0.01, 0.02, 0.04]]
5589         assert allclose(domain.quantities['friction'].vertex_values,
5590                        answer)
5591
5592         #print domain.quantities['friction'].vertex_values
5593         assert allclose(domain.tagged_elements['dsg'][0],0)
5594         assert allclose(domain.tagged_elements['ole nielsen'][0],1)
5595
5596         self.failUnless( domain.boundary[(1, 0)]  == '1',
5597                          "test_tags_to_boundaries  failed. Single boundary wasn't added.")
5598         self.failUnless( domain.boundary[(1, 2)]  == '2',
5599                          "test_tags_to_boundaries  failed. Single boundary wasn't added.")
5600         self.failUnless( domain.boundary[(0, 1)]  == '3',
5601                          "test_tags_to_boundaries  failed. Single boundary wasn't added.")
5602         self.failUnless( domain.boundary[(0, 0)]  == 'exterior',
5603                          "test_tags_to_boundaries  failed. Single boundary wasn't added.")
5604         #print "domain.boundary",domain.boundary
5605         self.failUnless( len(domain.boundary)  == 4,
5606                          "test_pmesh2Domain Too many boundaries")
5607         #FIXME change to use get_xllcorner
5608         #print "d.geo_reference.xllcorner",domain.geo_reference.xllcorner
5609         self.failUnless(domain.geo_reference.xllcorner  == 140.0,
5610                          "bad geo_referece")
5611
5612
5613         #************
5614
5615   
5616         domain = Domain(fileName)
5617         
5618         #print "domain.tagged_elements", domain.tagged_elements
5619         ## check the quantities
5620         #print domain.quantities['elevation'].vertex_values
5621         answer = [[0., 8., 0.],
5622                   [0., 10., 8.]]
5623         assert allclose(domain.quantities['elevation'].vertex_values,
5624                        answer)
5625
5626         #print domain.quantities['stage'].vertex_values
5627         answer = [[0., 12., 10.],
5628                   [0., 10., 12.]]
5629         assert allclose(domain.quantities['stage'].vertex_values,
5630                        answer)
5631
5632         #print domain.quantities['friction'].vertex_values
5633         answer = [[0.01, 0.04, 0.03],
5634                   [0.01, 0.02, 0.04]]
5635         assert allclose(domain.quantities['friction'].vertex_values,
5636                        answer)
5637
5638         #print domain.quantities['friction'].vertex_values
5639         assert allclose(domain.tagged_elements['dsg'][0],0)
5640         assert allclose(domain.tagged_elements['ole nielsen'][0],1)
5641
5642         self.failUnless( domain.boundary[(1, 0)]  == '1',
5643                          "test_tags_to_boundaries  failed. Single boundary wasn't added.")
5644         self.failUnless( domain.boundary[(1, 2)]  == '2',
5645                          "test_tags_to_boundaries  failed. Single boundary wasn't added.")
5646         self.failUnless( domain.boundary[(0, 1)]  == '3',
5647                          "test_tags_to_boundaries  failed. Single boundary wasn't added.")
5648         self.failUnless( domain.boundary[(0, 0)]  == 'exterior',
5649                          "test_tags_to_boundaries  failed. Single boundary wasn't added.")
5650         #print "domain.boundary",domain.boundary
5651         self.failUnless( len(domain.boundary)  == 4,
5652                          "test_pmesh2Domain Too many boundaries")
5653         #FIXME change to use get_xllcorner
5654         #print "d.geo_reference.xllcorner",domain.geo_reference.xllcorner
5655         self.failUnless(domain.geo_reference.xllcorner  == 140.0,
5656                          "bad geo_referece")
5657         #************
5658         os.remove(fileName)
5659
5660        #-------------------------------------------------------------
5661
5662    def test_get_lone_vertices(self):
5663       
5664        a = [0.0, 0.0]
5665        b = [0.0, 2.0]
5666        c = [2.0,0.0]
5667        d = [0.0, 4.0]
5668        e = [2.0, 2.0]
5669        f = [4.0,0.0]
5670
5671        points = [a, b, c, d, e, f]
5672        #bac, bce, ecf, dbe
5673        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4] ]
5674        boundary = { (0, 0): 'Third',
5675                     (0, 2): 'First',
5676                     (2, 0): 'Second',
5677                     (2, 1): 'Second',
5678                     (3, 1): 'Second',
5679                     (3, 2): 'Third'}
5680
5681
5682        domain = Domain(points, vertices, boundary)
5683        #domain.check_integrity()
5684        domain.get_lone_vertices()
5685
5686       
5687    def test_fitting_using_shallow_water_domain(self):
5688       
5689        #Mesh in zone 56 (absolute coords)
5690
5691        x0 = 314036.58727982
5692        y0 = 6224951.2960092
5693
5694        a = [x0+0.0, y0+0.0]
5695        b = [x0+0.0, y0+2.0]
5696        c = [x0+2.0, y0+0.0]
5697        d = [x0+0.0, y0+4.0]
5698        e = [x0+2.0, y0+2.0]
5699        f = [x0+4.0, y0+0.0]
5700
5701        points = [a, b, c, d, e, f]
5702
5703        #bac, bce, ecf, dbe
5704        elements = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4] ]
5705
5706        #absolute going in ..
5707        mesh4 = Domain(points, elements,
5708                       geo_reference = Geo_reference(56, 0, 0))
5709        mesh4.check_integrity()
5710        quantity = Quantity(mesh4)
5711
5712        #Get (enough) datapoints (relative to georef)
5713        data_points_rel = [[ 0.66666667, 0.66666667],
5714                       [ 1.33333333, 1.33333333],
5715                       [ 2.66666667, 0.66666667],
5716                       [ 0.66666667, 2.66666667],
5717                       [ 0.0, 1.0],
5718                       [ 0.0, 3.0],
5719                       [ 1.0, 0.0],
5720                       [ 1.0, 1.0],
5721                       [ 1.0, 2.0],
5722                       [ 1.0, 3.0],
5723                       [ 2.0, 1.0],
5724                       [ 3.0, 0.0],
5725                       [ 3.0, 1.0]]
5726
5727        data_geo_spatial = Geospatial_data(data_points_rel,
5728                         geo_reference = Geo_reference(56, x0, y0))
5729        data_points_absolute = data_geo_spatial.get_data_points(absolute=True)
5730        attributes = linear_function(data_points_absolute)
5731        att = 'spam_and_eggs'
5732       
5733        #Create .txt file
5734        ptsfile = tempfile.mktemp(".txt")
5735        file = open(ptsfile,"w")
5736        file.write(" x,y," + att + " \n")
5737        for data_point, attribute in map(None, data_points_absolute
5738                                         ,attributes):
5739            row = str(data_point[0]) + ',' + str(data_point[1]) \
5740                  + ',' + str(attribute)
5741            file.write(row + "\n")
5742        file.close()
5743
5744        #file = open(ptsfile, 'r')
5745        #lines = file.readlines()
5746        #file.close()
5747     
5748
5749        #Check that values can be set from file
5750        quantity.set_values(filename = ptsfile,
5751                            attribute_name = att, alpha = 0)
5752        answer = linear_function(quantity.domain.get_vertex_coordinates())
5753
5754        assert allclose(quantity.vertex_values.flat, answer)
5755
5756
5757        #Check that values can be set from file using default attribute
5758        quantity.set_values(filename = ptsfile, alpha = 0)
5759        assert allclose(quantity.vertex_values.flat, answer)
5760
5761        #Cleanup
5762        import os
5763        os.remove(ptsfile)
5764
5765
5766       
5767if __name__ == "__main__":
5768
5769    suite = unittest.makeSuite(Test_Shallow_Water,'test')
5770
5771    #suite = unittest.makeSuite(Test_Shallow_Water,'test_bedslope_problem_first_order_moresteps')   
5772    #suite = unittest.makeSuite(Test_Shallow_Water,'test_fitting_using_shallow_water_domain')   
5773    #suite = unittest.makeSuite(Test_Shallow_Water,'test_tight_slope_limiters')
5774    #suite = unittest.makeSuite(Test_Shallow_Water,'test_get_maximum_inundation_from_sww')
5775    #suite = unittest.makeSuite(Test_Shallow_Water,'test_temp')   
5776   
5777
5778   
5779    runner = unittest.TextTestRunner(verbosity=1)   
5780    runner.run(suite)
Note: See TracBrowser for help on using the repository browser.