source: anuga_core/source/anuga/shallow_water_balanced/test_swb_distribute.py @ 7559

Last change on this file since 7559 was 7559, checked in by steve, 10 years ago

Broke up test_swb_domain so that we can concentrate on new functionality

File size: 14.4 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.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
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_clean(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._order_ = 1
150        domain.distribute_to_vertices_and_edges()
151        assert num.allclose(L[1], val1)
152
153        # Second order
154        domain._order_ = 2
155        domain.beta_w = 0.9
156        domain.beta_w_dry = 0.9
157        domain.beta_uh = 0.9
158        domain.beta_uh_dry = 0.9
159        domain.beta_vh = 0.9
160        domain.beta_vh_dry = 0.9
161        domain.distribute_to_vertices_and_edges()
162        assert num.allclose(L[1], [2.2, 4.9, 4.9])
163
164    def test_distribute_away_from_bed(self):
165        #Using test data generated by abstract_2d_finite_volumes-2
166        #Assuming no friction and flat bed (0.0)
167
168        a = [0.0, 0.0]
169        b = [0.0, 2.0]
170        c = [2.0, 0.0]
171        d = [0.0, 4.0]
172        e = [2.0, 2.0]
173        f = [4.0, 0.0]
174
175        points = [a, b, c, d, e, f]
176        #             bac,     bce,     ecf,     dbe
177        vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
178
179        domain = Domain(points, vertices)
180        L = domain.quantities['stage'].vertex_values
181
182        def stage(x, y):
183            return x**2
184
185        domain.set_quantity('stage', stage, location='centroids')
186
187        domain.quantities['stage'].compute_gradients()
188
189        a, b = domain.quantities['stage'].get_gradients()
190
191        assert num.allclose(a[1], 3.33333334)
192        assert num.allclose(b[1], 0.0)
193
194        domain._order_ = 1
195        domain.distribute_to_vertices_and_edges()
196        assert num.allclose(L[1], 1.77777778)
197
198        domain._order_ = 2
199        domain.beta_w = 0.9
200        domain.beta_w_dry = 0.9
201        domain.beta_uh = 0.9
202        domain.beta_uh_dry = 0.9
203        domain.beta_vh = 0.9
204        domain.beta_vh_dry = 0.9
205        domain.distribute_to_vertices_and_edges()
206        assert num.allclose(L[1], [0.57777777, 2.37777778, 2.37777778])
207
208    def test_distribute_away_from_bed1(self):
209        #Using test data generated by abstract_2d_finite_volumes-2
210        #Assuming no friction and flat bed (0.0)
211
212        a = [0.0, 0.0]
213        b = [0.0, 2.0]
214        c = [2.0, 0.0]
215        d = [0.0, 4.0]
216        e = [2.0, 2.0]
217        f = [4.0, 0.0]
218
219        points = [a, b, c, d, e, f]
220        #             bac,     bce,     ecf,     dbe
221        vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
222
223        domain = Domain(points, vertices)
224        L = domain.quantities['stage'].vertex_values
225
226        def stage(x, y):
227            return x**4 + y**2
228
229        domain.set_quantity('stage', stage, location='centroids')
230
231        domain.quantities['stage'].compute_gradients()
232        a, b = domain.quantities['stage'].get_gradients()
233        assert num.allclose(a[1], 25.18518519)
234        assert num.allclose(b[1], 3.33333333)
235
236        domain._order_ = 1
237        domain.distribute_to_vertices_and_edges()
238        assert num.allclose(L[1], 4.9382716)
239
240        domain._order_ = 2
241        domain.beta_w = 0.9
242        domain.beta_w_dry = 0.9
243        domain.beta_uh = 0.9
244        domain.beta_uh_dry = 0.9
245        domain.beta_vh = 0.9
246        domain.beta_vh_dry = 0.9
247        domain.distribute_to_vertices_and_edges()
248        assert num.allclose(L[1], [1.07160494, 6.46058131, 7.28262855])
249
250    def test_distribute_near_bed(self):
251        a = [0.0, 0.0]
252        b = [0.0, 2.0]
253        c = [2.0, 0.0]
254        d = [0.0, 4.0]
255        e = [2.0, 2.0]
256        f = [4.0, 0.0]
257
258        points = [a, b, c, d, e, f]
259        #             bac,     bce,     ecf,     dbe
260        vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
261
262        domain = Domain(points, vertices)
263
264        # Set up for a gradient of (10,0) at mid triangle (bce)
265        def slope(x, y):
266            return 10*x
267
268        h = 0.1
269        def stage(x, y):
270            return slope(x, y) + h
271
272        domain.set_quantity('elevation', slope)
273        domain.set_quantity('stage', stage, location='centroids')
274
275        E = domain.quantities['elevation'].vertex_values
276        L = domain.quantities['stage'].vertex_values
277
278        # Get reference values
279        volumes = []
280        for i in range(len(L)):
281            volumes.append(num.sum(L[i])/3)
282            assert num.allclose(volumes[i],
283                                domain.quantities['stage'].centroid_values[i])
284
285        domain._order_ = 1
286
287        domain.tight_slope_limiters = 0
288        domain.distribute_to_vertices_and_edges()
289        assert num.allclose(L[1], [0.1, 20.1, 20.1])
290        for i in range(len(L)):
291            assert num.allclose(volumes[i], num.sum(L[i])/3)
292
293        # Allow triangle to be flatter (closer to bed)
294        domain.tight_slope_limiters = 1
295
296        domain.distribute_to_vertices_and_edges()
297        assert num.allclose(L[1], [0.298, 20.001, 20.001])
298        for i in range(len(L)):
299            assert num.allclose(volumes[i], num.sum(L[i])/3)
300
301        domain._order_ = 2
302
303        domain.tight_slope_limiters = 0
304        domain.distribute_to_vertices_and_edges()
305        assert num.allclose(L[1], [0.1, 20.1, 20.1])
306        for i in range(len(L)):
307            assert num.allclose(volumes[i], num.sum(L[i])/3)
308
309        # Allow triangle to be flatter (closer to bed)
310        domain.tight_slope_limiters = 1
311
312        domain.distribute_to_vertices_and_edges()
313        assert num.allclose(L[1], [0.298, 20.001, 20.001])
314        for i in range(len(L)):
315            assert num.allclose(volumes[i], num.sum(L[i])/3)
316
317    def test_distribute_near_bed1(self):
318        a = [0.0, 0.0]
319        b = [0.0, 2.0]
320        c = [2.0, 0.0]
321        d = [0.0, 4.0]
322        e = [2.0, 2.0]
323        f = [4.0, 0.0]
324
325        points = [a, b, c, d, e, f]
326        #             bac,     bce,     ecf,     dbe
327        vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
328
329        domain = Domain(points, vertices)
330
331        # Set up for a gradient of (8,2) at mid triangle (bce)
332        def slope(x, y):
333            return x**4 + y**2
334
335        h = 0.1
336        def stage(x, y):
337            return slope(x, y) + h
338
339        domain.set_quantity('elevation', slope)
340        domain.set_quantity('stage', stage)
341
342        E = domain.quantities['elevation'].vertex_values
343        L = domain.quantities['stage'].vertex_values
344
345        # Get reference values
346        volumes = []
347        for i in range(len(L)):
348            volumes.append(num.sum(L[i])/3)
349            assert num.allclose(volumes[i],
350                                domain.quantities['stage'].centroid_values[i])
351
352        domain._order_ = 1
353
354        domain.tight_slope_limiters = 0
355        domain.distribute_to_vertices_and_edges()
356        assert num.allclose(L[1], [4.1, 16.1, 20.1])       
357        for i in range(len(L)):
358            assert num.allclose(volumes[i], num.sum(L[i])/3)
359       
360               
361        domain.tight_slope_limiters = 1 # Allow triangle to be flatter (closer to bed)
362        domain.distribute_to_vertices_and_edges()
363        assert num.allclose(L[1], [4.2386, 16.0604, 20.001])
364        for i in range(len(L)):
365            assert num.allclose(volumes[i], num.sum(L[i])/3)   
366       
367
368        domain._order_ = 2
369       
370        domain.tight_slope_limiters = 0   
371        domain.distribute_to_vertices_and_edges()
372        assert num.allclose(L[1], [4.1, 16.1, 20.1])
373        for i in range(len(L)):
374            assert num.allclose(volumes[i], num.sum(L[i])/3)
375
376        # Allow triangle to be flatter (closer to bed)
377        domain.tight_slope_limiters = 1
378
379        domain.distribute_to_vertices_and_edges()
380        assert num.allclose(L[1], [4.23370103, 16.06529897, 20.001]) or\
381               num.allclose(L[1], [4.18944138, 16.10955862, 20.001]) or\
382               num.allclose(L[1], [4.19351461, 16.10548539, 20.001]) # old limiters
383       
384        for i in range(len(L)):
385            assert num.allclose(volumes[i], num.sum(L[i])/3)
386
387
388    def test_second_order_distribute_real_data(self):
389        #Using test data generated by abstract_2d_finite_volumes-2
390        #Assuming no friction and flat bed (0.0)
391
392        a = [0.0, 0.0]
393        b = [0.0, 1.0/5]
394        c = [0.0, 2.0/5]
395        d = [1.0/5, 0.0]
396        e = [1.0/5, 1.0/5]
397        f = [1.0/5, 2.0/5]
398        g = [2.0/5, 2.0/5]
399
400        points = [a, b, c, d, e, f, g]
401        #             bae,     efb,     cbf,     feg
402        vertices = [[1,0,4], [4,5,1], [2,1,5], [5,4,6]]
403
404        domain = Domain(points, vertices)
405
406        def slope(x, y):
407            return -x/3
408
409        domain.set_quantity('elevation', slope)
410        domain.set_quantity('stage',
411                            [0.01298164, 0.00365611,
412                             0.01440365, -0.0381856437096],
413                            location='centroids')
414        domain.set_quantity('xmomentum',
415                            [0.00670439, 0.01263789,
416                             0.00647805, 0.0178180740668],
417                            location='centroids')
418        domain.set_quantity('ymomentum',
419                            [-7.23510980e-004, -6.30413883e-005,
420                             6.30413883e-005, 0.000200907255866],
421                            location='centroids')
422
423        E = domain.quantities['elevation'].vertex_values
424        L = domain.quantities['stage'].vertex_values
425        X = domain.quantities['xmomentum'].vertex_values
426        Y = domain.quantities['ymomentum'].vertex_values
427
428        domain._order_ = 2
429        domain.beta_w = 0.9
430        domain.beta_w_dry = 0.9
431        domain.beta_uh = 0.9
432        domain.beta_uh_dry = 0.9
433        domain.beta_vh = 0.9
434        domain.beta_vh_dry = 0.9
435
436        # FIXME (Ole): Need tests where this is commented out
437        domain.tight_slope_limiters = 0       # Backwards compatibility (14/4/7)
438        domain.use_centroid_velocities = 0    # Backwards compatibility (7/5/8)
439
440        domain.distribute_to_vertices_and_edges()
441
442        assert num.allclose(L[1,:], [-0.00825735775384,
443                                     -0.00801881482869,
444                                     0.0272445025825])
445        assert num.allclose(X[1,:], [0.0143507718962,
446                                     0.0142502147066,
447                                     0.00931268339717])
448        assert num.allclose(Y[1,:], [-0.000117062180693,
449                                     7.94434448109e-005,
450                                     -0.000151505429018])
451
452
453
454#################################################################################
455
456if __name__ == "__main__":
457    suite = unittest.makeSuite(Test_swb_clean, 'test')
458    runner = unittest.TextTestRunner(verbosity=1)
459    runner.run(suite)
Note: See TracBrowser for help on using the repository browser.