source: trunk/anuga_core/source/anuga/shallow_water_balanced/test_swb_distribute.py @ 8050

Last change on this file since 8050 was 7866, checked in by hudson, 14 years ago

More swb tests passing. Cleaned up some pylint errors.

File size: 12.5 KB
Line 
1#!/usr/bin/env python
2
3import unittest, os
4import os.path
5from math import pi, sqrt
6import tempfile
7
8from anuga.config import g, epsilon
9from anuga.config import netcdf_mode_r, netcdf_mode_w, netcdf_mode_a
10from anuga.utilities.numerical_tools import mean
11from anuga.geometry.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
16
17from anuga.utilities.system_tools import get_pathname_from_package
18from swb_domain import *
19
20import numpy as num
21
22# Get gateway to C implementation of flux function for direct testing
23from shallow_water_ext import flux_function_central as flux_function
24
25
26
27
28class Test_swb_distribute(unittest.TestCase):
29    def setUp(self):
30        pass
31
32    def tearDown(self):
33        pass
34
35
36
37
38    def test_first_order_extrapolator_const_z(self):
39        a = [0.0, 0.0]
40        b = [0.0, 2.0]
41        c = [2.0, 0.0]
42        d = [0.0, 4.0]
43        e = [2.0, 2.0]
44        f = [4.0, 0.0]
45
46        points = [a, b, c, d, e, f]
47        #             bac,     bce,     ecf,     dbe
48        vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
49
50        domain = Domain(points, vertices)
51        val0 = 2. + 2.0/3
52        val1 = 4. + 4.0/3
53        val2 = 8. + 2.0/3
54        val3 = 2. + 8.0/3
55
56        zl = zr = -3.75    # Assume constant bed (must be less than stage)
57        domain.set_quantity('elevation', zl*num.ones((4, 3), num.int)) #array default#
58        domain.set_quantity('stage', [[val0, val0-1, val0-2],
59                                      [val1, val1+1, val1],
60                                      [val2, val2-2, val2],
61                                      [val3-0.5, val3, val3]])
62
63        domain._order_ = 1
64        domain.distribute_to_vertices_and_edges()
65
66        #Check that centroid values were distributed to vertices
67        C = domain.quantities['stage'].centroid_values
68        for i in range(3):
69            assert num.allclose(domain.quantities['stage'].vertex_values[:,i],
70                                C)
71
72
73
74
75    def test_first_order_limiter_variable_z(self):
76        '''Check that first order limiter follows bed_slope'''
77
78        from anuga.config import epsilon
79
80        a = [0.0, 0.0]
81        b = [0.0, 2.0]
82        c = [2.0,0.0]
83        d = [0.0, 4.0]
84        e = [2.0, 2.0]
85        f = [4.0,0.0]
86
87        points = [a, b, c, d, e, f]
88        #bac, bce, ecf, dbe
89        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
90
91        domain = Domain(points, vertices)
92        val0 = 2.+2.0/3
93        val1 = 4.+4.0/3
94        val2 = 8.+2.0/3
95        val3 = 2.+8.0/3
96
97        domain.set_quantity('elevation', [[0,0,0], [6,0,0],
98                                          [6,6,6], [6,6,6]])
99        domain.set_quantity('stage', [[val0, val0, val0],
100                                      [val1, val1, val1],
101                                      [val2, val2, val2],
102                                      [val3, val3, val3]])
103
104        E = domain.quantities['elevation'].vertex_values
105        L = domain.quantities['stage'].vertex_values
106
107
108        #Check that some stages are not above elevation (within eps)
109        #- so that the limiter has something to work with
110        assert not num.alltrue(num.alltrue(num.greater_equal(L,E-epsilon)))
111
112        domain._order_ = 1
113        domain.distribute_to_vertices_and_edges()
114
115        #Check that all stages are above elevation (within eps)
116        assert num.alltrue(num.alltrue(num.greater_equal(L,E-epsilon)))
117
118
119
120
121
122    def test_distribute_basic(self):
123        #Using test data generated by abstract_2d_finite_volumes-2
124        #Assuming no friction and flat bed (0.0)
125
126        a = [0.0, 0.0]
127        b = [0.0, 2.0]
128        c = [2.0, 0.0]
129        d = [0.0, 4.0]
130        e = [2.0, 2.0]
131        f = [4.0, 0.0]
132
133        points = [a, b, c, d, e, f]
134        #             bac,     bce,     ecf,     dbe
135        vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
136
137        domain = Domain(points, vertices)
138
139        val0 = 2.
140        val1 = 4.
141        val2 = 8.
142        val3 = 2.
143
144        domain.set_quantity('stage', [val0, val1, val2, val3],
145                            location='centroids')
146        L = domain.quantities['stage'].vertex_values
147
148        # First order
149        domain.set_default_order(1)
150        domain.distribute_to_vertices_and_edges()
151       
152        assert num.allclose(L[1], val1)
153
154        # Second order
155        domain.set_default_order(2)
156        domain.distribute_to_vertices_and_edges()
157
158        assert num.allclose(L[1], [0.0,   6.0,  6.0], atol=2.0e-3 )
159
160        assert num.allclose(val1, num.sum(L[1])/3)
161
162    def test_distribute_away_from_bed(self):
163        #Using test data generated by abstract_2d_finite_volumes-2
164        #Assuming no friction and flat bed (0.0)
165
166        a = [0.0, 0.0]
167        b = [0.0, 2.0]
168        c = [2.0, 0.0]
169        d = [0.0, 4.0]
170        e = [2.0, 2.0]
171        f = [4.0, 0.0]
172
173        points = [a, b, c, d, e, f]
174        #             bac,     bce,     ecf,     dbe
175        vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
176
177        domain = Domain(points, vertices)
178        L = domain.quantities['stage'].vertex_values
179
180        def stage(x, y):
181            return x**2
182
183        domain.set_quantity('stage', stage, location='centroids')
184        domain.set_quantity('elevation',-3.0)
185
186        domain.quantities['stage'].compute_gradients()
187
188        a, b = domain.quantities['stage'].get_gradients()
189
190        assert num.allclose(a[1], 3.33333334)
191        assert num.allclose(b[1], 0.0)
192
193        domain.set_default_order(1)
194        domain.distribute_to_vertices_and_edges()
195
196        f1 = stage(4.0/3.0, 4.0/3.0)
197        assert num.allclose(L[1], f1)
198
199        domain.set_default_order(2)
200        domain.distribute_to_vertices_and_edges()
201
202
203        fv0 = f1 - a[1]*4.0/3.0 + b[1]*2.0/3.0
204        fv1 = f1 + a[1]*2.0/3.0 - b[1]*4.0/3.0
205        fv2 = f1 + a[1]*2.0/3.0 + b[1]*2.0/3.0
206
207        assert num.allclose(L[1], [fv0,fv1,fv2])
208
209        assert num.allclose(f1, num.sum(L[1])/3)
210
211    def test_distribute_away_from_bed1(self):
212        #Using test data generated by abstract_2d_finite_volumes-2
213        #Assuming no friction and flat bed (0.0)
214
215        a = [0.0, 0.0]
216        b = [0.0, 2.0]
217        c = [2.0, 0.0]
218        d = [0.0, 4.0]
219        e = [2.0, 2.0]
220        f = [4.0, 0.0]
221
222        points = [a, b, c, d, e, f]
223        #             bac,     bce,     ecf,     dbe
224        vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
225
226        domain = Domain(points, vertices)
227        L = domain.quantities['stage'].vertex_values
228
229        def stage(x, y):
230            return x**4 + y**2
231
232        domain.set_quantity('stage', stage, location='centroids')
233        domain.set_quantity('elevation', -10.0)
234
235        domain.quantities['stage'].compute_gradients()
236        a, b = domain.quantities['stage'].get_gradients()
237        assert num.allclose(a[1], 25.18518519)
238        assert num.allclose(b[1], 3.33333333)
239
240        domain.set_default_order(1)
241        domain.distribute_to_vertices_and_edges()
242
243        f1 = stage(4.0/3.0, 4.0/3.0)
244        assert num.allclose(L[1], f1)
245
246        domain.set_default_order(2)
247        domain.distribute_to_vertices_and_edges()
248
249
250        fv0 = f1 - a[1]*4.0/3.0 + b[1]*2.0/3.0
251        fv1 = f1 + a[1]*2.0/3.0 - b[1]*4.0/3.0
252        fv2 = f1 + a[1]*2.0/3.0 + b[1]*2.0/3.0
253
254       
255        assert num.allclose(L[1], [ fv0, fv1, fv2]) or \
256               num.allclose(L[1], [ -9.23392657,  10.51787718,  13.5308642 ]) 
257
258
259    def test_distribute_near_bed(self):
260        a = [0.0, 0.0]
261        b = [0.0, 2.0]
262        c = [1.0, 1.0]
263        d = [2.0, 0.0]
264        e = [2.0, 2.0]
265   
266
267        points = [a, b, c, d, e]
268
269        vertices = [[0,3,2], [0,2,1], [2,3,4], [1,2,4]]
270
271        domain = Domain(points, vertices)
272
273        # Set up for a gradient of (10,0) at mid triangle (bce)
274        def slope(x, y):
275            return 10*x
276
277        h = 0.1
278        def stage(x, y):
279            return slope(x, y) + h
280
281        domain.set_quantity('elevation', slope, location='centroids')
282        domain.set_quantity('stage', stage, location='centroids')
283
284
285
286        E = domain.quantities['elevation']
287        L = domain.quantities['stage']
288        Z = domain.quantities['elevation']
289
290        E_V = E.vertex_values
291        L_V = L.vertex_values
292
293        E_E = E.edge_values
294        L_E = L.edge_values       
295        E_C = E.centroid_values
296        L_C = L.centroid_values
297
298
299        domain.set_default_order(1)
300        domain.distribute_to_vertices_and_edges()
301
302        assert num.allclose(L_V,[[ 10.1,         10.1,         10.1       ],
303                                 [  3.43333333,   3.43333333,   3.43333333],
304                                 [ 16.76666667,  16.76666667,  16.76666667],
305                                 [ 10.1,         10.1,         10.1       ]])
306
307        assert num.allclose(E_V,[[ 10.,          10.,          10.,        ],
308                                 [  3.33333333,   3.33333333,   3.33333333],
309                                 [ 16.66666667,  16.66666667,  16.66666667],
310                                 [ 10.,          10.,          10.        ]])
311
312        domain.set_default_order(2)
313
314        # Setup the elevation to be pw linear (actually linear)
315        Z.extrapolate_second_order_and_limit_by_edge()
316
317
318        domain.distribute_to_vertices_and_edges()
319
320
321
322        assert num.allclose(L_V,[[  0.1,  20.1,  10.1],
323                                 [  0.1,  10.1,   0.1],
324                                 [ 10.1,  20.1,  20.1],
325                                 [  0.1,  10.1,  20.1]]) or \
326               num.allclose(L_V,[[  0.1,         20.1,         10.1,       ],
327                                 [  3.43333333,   3.43333333,   3.43333333],
328                                 [ 16.76666667,  16.76666667,  16.76666667],
329                                 [  0.1,         10.1,         20.1       ]])
330       
331
332        assert num.allclose(E_V,[[  0.,  20.,  10.],
333                                 [  0.,  10.,   0.],
334                                 [ 10.,  20.,  20.],
335                                 [  0.,  10.,  20.]]) or \
336               num.allclose(E_V,[[  0.,          20.,          10.,        ],
337                                 [  3.33333333,   3.33333333,   3.33333333],
338                                 [ 16.66666667,  16.66666667,  16.66666667],
339                                 [  0.,          10.,          20.        ]])
340
341                                 
342
343       
344
345    def test_second_order_distribute_real_data(self):
346        #Using test data generated by abstract_2d_finite_volumes-2
347        #Assuming no friction and flat bed (0.0)
348
349        a = [0.0, 0.0]
350        b = [0.0, 1.0/5]
351        c = [0.0, 2.0/5]
352        d = [1.0/5, 0.0]
353        e = [1.0/5, 1.0/5]
354        f = [1.0/5, 2.0/5]
355        g = [2.0/5, 2.0/5]
356
357        points = [a, b, c, d, e, f, g]
358        #             bae,     efb,     cbf,     feg
359        vertices = [[1,0,4], [4,5,1], [2,1,5], [5,4,6]]
360
361        domain = Domain(points, vertices)
362
363        def slope(x, y):
364            return -x/3
365
366        domain.set_quantity('elevation', slope)
367        domain.set_quantity('stage',
368                            [0.01298164, 0.00365611,
369                             0.01440365, -0.0381856437096],
370                            location='centroids')
371        domain.set_quantity('xmomentum',
372                            [0.00670439, 0.01263789,
373                             0.00647805, 0.0178180740668],
374                            location='centroids')
375        domain.set_quantity('ymomentum',
376                            [-7.23510980e-004, -6.30413883e-005,
377                             6.30413883e-005, 0.000200907255866],
378                            location='centroids')
379
380        E = domain.quantities['elevation'].vertex_values
381        L = domain.quantities['stage'].vertex_values
382        X = domain.quantities['xmomentum'].vertex_values
383        Y = domain.quantities['ymomentum'].vertex_values
384
385        domain.set_default_order(2)
386
387        domain.distribute_to_vertices_and_edges()
388
389
390        assert num.allclose(L[1,:],
391                            [-0.01434766, -0.01292565, 0.03824164], atol=1.0e-2)
392
393        assert num.allclose(X[1,:],
394                            [ 0.01702702, 0.01676034,  0.0057706 ], atol=1.0e-2)
395
396        assert num.allclose(Y[1,:],
397                            [-0.00041792,  0.00076771, -0.00039118], atol=1.0e-4)
398
399
400
401#################################################################################
402
403if __name__ == "__main__":
404    suite = unittest.makeSuite(Test_swb_distribute, 'test')
405    runner = unittest.TextTestRunner(verbosity=1)
406    runner.run(suite)
Note: See TracBrowser for help on using the repository browser.