source: branches/numpy/anuga/shallow_water/test_shallow_water_domain.py @ 6410

Last change on this file since 6410 was 6410, checked in by rwilson, 15 years ago

numpy changes.

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