source: branches/anuga_1_1/anuga_core/source/anuga/shallow_water_balanced/test_swb_balance.py @ 7875

Last change on this file since 7875 was 7616, checked in by steve, 14 years ago

Updating to relocated repository

File size: 13.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    def test_balance_deep_and_shallow(self):
37        """Test that balanced limiters preserve conserved quantites.
38        This test is using old depth based balanced limiters
39        """
40
41        import copy
42
43        a = [0.0, 0.0]
44        b = [0.0, 2.0]
45        c = [2.0, 0.0]
46        d = [0.0, 4.0]
47        e = [2.0, 2.0]
48        f = [4.0, 0.0]
49
50        points = [a, b, c, d, e, f]
51        #             bac,     bce,     ecf,     dbe
52        elements = [[1,0,2], [1,2,4], [4,2,5], [3,1,4] ]
53
54        domain = Domain(points, elements)
55        domain.check_integrity()
56
57        # Create a deliberate overshoot
58        domain.set_quantity('stage', [[3,0,3], [2,2,6], [5,3,8], [8,3,5]])
59        domain.set_quantity('elevation', 0)    # Flat bed
60        stage = domain.quantities['stage']
61
62        ref_centroid_values = copy.copy(stage.centroid_values[:])    # Copy
63
64        # Limit
65        domain.distribute_to_vertices_and_edges()
66
67        # Assert that quantities are conserved
68        for k in range(len(domain)):
69            assert num.allclose(ref_centroid_values[k],
70                                num.sum(stage.vertex_values[k,:])/3)
71
72        # Now try with a non-flat bed - closely hugging initial stage in places
73        domain.set_quantity('stage', [[3,0,3], [2,2,6], [5,3,8], [8,3,5]])
74        domain.set_quantity('elevation', [[0,0,0],
75                                          [1.8,1.9,5.9],
76                                          [4.6,0,0],
77                                          [0,2,4]])
78        stage = domain.quantities['stage']
79        elevation = domain.quantities['elevation']
80        height = domain.quantities['height']
81
82        ref_centroid_values = copy.copy(stage.centroid_values[:])    # Copy
83        ref_vertex_values = copy.copy(stage.vertex_values[:])        # Copy
84
85        # Limit
86        domain.distribute_to_vertices_and_edges()
87
88        # Assert that all vertex quantities have changed
89        for k in range(len(domain)):
90            assert not num.allclose(ref_vertex_values[k,:],
91                                    stage.vertex_values[k,:])
92        # and assert that quantities are still conserved
93        for k in range(len(domain)):
94            assert num.allclose(ref_centroid_values[k],
95                                num.sum(stage.vertex_values[k,:])/3)
96
97        # Check actual results
98        assert num.allclose(stage.vertex_values,
99                            [[ 2.66666667,  0.66666667,  2.66666667],
100                             [ 3.33333333,  3.33333333,  3.33333333],
101                             [ 3.73333333,  4.93333333,  7.33333333],
102                             [ 7.33333333,  4.93333333,  3.73333333]])
103
104
105
106    def test_balance_deep_and_shallow_tight_SL(self):
107        """Test that balanced limiters preserve conserved quantites.
108        This test is using Tight Slope Limiters
109        """
110
111        import copy
112
113        a = [0.0, 0.0]
114        b = [0.0, 2.0]
115        c = [2.0, 0.0]
116        d = [0.0, 4.0]
117        e = [2.0, 2.0]
118        f = [4.0, 0.0]
119
120        points = [a, b, c, d, e, f]
121        #             bac,     bce,     ecf,     dbe
122        elements = [[1,0,2], [1,2,4], [4,2,5], [3,1,4] ]
123
124        domain = Domain(points, elements)
125        domain.check_integrity()
126
127        # Create a deliberate overshoot
128        domain.set_quantity('stage', [[3,0,3], [2,2,6], [5,3,8], [8,3,5]])
129        domain.set_quantity('elevation', 0)    # Flat bed
130        stage = domain.quantities['stage']
131
132        ref_centroid_values = copy.copy(stage.centroid_values[:])    # Copy
133
134        # Limit
135        domain.tight_slope_limiters = 1
136        domain.distribute_to_vertices_and_edges()
137
138        # Assert that quantities are conserved
139        for k in range(len(domain)):
140            assert num.allclose (ref_centroid_values[k],
141                                 num.sum(stage.vertex_values[k,:])/3)
142
143        # Now try with a non-flat bed - closely hugging initial stage in places
144        # This will create alphas in the range [0, 0.478260, 1]
145        domain.set_quantity('stage', [[3,0,3], [2,2,6], [5,3,8], [8,3,5]])
146        domain.set_quantity('elevation', [[0,0,0],
147                                          [1.8,1.9,5.9],
148                                          [4.6,0,0],
149                                          [0,2,4]])
150        stage = domain.quantities['stage']
151
152        ref_centroid_values = copy.copy(stage.centroid_values[:])    # Copy
153        ref_vertex_values = copy.copy(stage.vertex_values[:])        # Copy
154
155        # Limit
156        domain.tight_slope_limiters = 1
157        domain.distribute_to_vertices_and_edges()
158
159        # Assert that all vertex quantities have changed
160        for k in range(len(domain)):
161            assert not num.allclose(ref_vertex_values[k,:],
162                                    stage.vertex_values[k,:])
163        # and assert that quantities are still conserved
164        for k in range(len(domain)):
165            assert num.allclose(ref_centroid_values[k],
166                                num.sum(stage.vertex_values[k,:])/3)
167
168    def test_balance_deep_and_shallow_Froude(self):
169        """Test that balanced limiters preserve conserved quantites -
170        and also that excessive Froude numbers are dealt with.
171        This test is using tight slope limiters.
172        """
173
174        import copy
175
176        a = [0.0, 0.0]
177        b = [0.0, 2.0]
178        c = [2.0, 0.0]
179        d = [0.0, 4.0]
180        e = [2.0, 2.0]
181        f = [4.0, 0.0]
182
183        points = [a, b, c, d, e, f]
184        #             bac,     bce,     ecf,     dbe
185        elements = [[1,0,2], [1,2,4], [4,2,5], [3,1,4] ]
186
187        domain = Domain(points, elements)
188        domain.check_integrity()
189        domain.tight_slope_limiters = True
190        domain.use_centroid_velocities = True
191
192        # Create non-flat bed - closely hugging initial stage in places
193        # This will create alphas in the range [0, 0.478260, 1]
194        domain.set_quantity('stage', [[3,0,3], [2,2,6], [5,3,8], [8,3,5]])
195        domain.set_quantity('elevation', [[0,0,0],
196                                          [1.8,1.999,5.999],
197                                          [4.6,0,0],
198                                          [0,2,4]])
199
200        # Create small momenta, that nonetheless will generate large speeds
201        # due to shallow depth at isolated vertices
202        domain.set_quantity('xmomentum', -0.0058)
203        domain.set_quantity('ymomentum', 0.0890)
204
205        stage = domain.quantities['stage']
206        elevation = domain.quantities['elevation']
207        xmomentum = domain.quantities['xmomentum']
208        ymomentum = domain.quantities['ymomentum']
209
210        # Setup triangle #1 to mimick real Froude explosion observed
211        # in the Onslow example 13 Nov 2007.
212        stage.vertex_values[1,:] = [1.6385, 1.6361, 1.2953]
213        elevation.vertex_values[1,:] = [1.6375, 1.6336, 0.4647]
214        xmomentum.vertex_values[1,:] = [-0.0058, -0.0050, -0.0066]
215        ymomentum.vertex_values[1,:] = [0.0890, 0.0890, 0.0890]
216
217        xmomentum.interpolate()
218        ymomentum.interpolate()
219        stage.interpolate()
220        elevation.interpolate()
221
222        # Verify interpolation
223        assert num.allclose(stage.centroid_values[1], 1.5233)
224        assert num.allclose(elevation.centroid_values[1], 1.2452667)
225        assert num.allclose(xmomentum.centroid_values[1], -0.0058)
226        assert num.allclose(ymomentum.centroid_values[1], 0.089)
227
228        # Derived quantities
229        depth = stage-elevation
230        u = xmomentum/depth
231        v = ymomentum/depth
232
233        denom = (depth*g)**0.5
234        Fx = u/denom
235        Fy = v/denom
236
237        # Verify against Onslow example (14 Nov 2007)
238        assert num.allclose(depth.centroid_values[1], 0.278033)
239        assert num.allclose(u.centroid_values[1], -0.0208608)
240        assert num.allclose(v.centroid_values[1], 0.3201055)
241
242        assert num.allclose(denom.centroid_values[1],
243                            num.sqrt(depth.centroid_values[1]*g))
244
245        assert num.allclose(u.centroid_values[1]/denom.centroid_values[1],
246                            -0.012637746977)
247        assert num.allclose(Fx.centroid_values[1],
248                            u.centroid_values[1]/denom.centroid_values[1])
249
250        # Check that Froude numbers are small at centroids.
251        assert num.allclose(Fx.centroid_values[1], -0.012637746977)
252        assert num.allclose(Fy.centroid_values[1], 0.193924048435)
253
254        # But Froude numbers are huge at some vertices and edges
255        assert num.allclose(Fx.vertex_values[1,:], [-5.85888475e+01,
256                                                    -1.27775313e+01,
257                                                    -2.78511420e-03])
258
259        assert num.allclose(Fx.edge_values[1,:], [-6.89150773e-03,
260                                                  -7.38672488e-03,
261                                                  -2.35626238e+01])
262
263        assert num.allclose(Fy.vertex_values[1,:], [8.99035764e+02,
264                                                    2.27440057e+02,
265                                                    3.75568430e-02])
266
267        assert num.allclose(Fy.edge_values[1,:], [1.05748998e-01,
268                                                  1.06035244e-01,
269                                                  3.88346947e+02])
270
271
272        # The task is now to arrange the limiters such that Froude numbers
273        # remain under control whil at the same time obeying the conservation
274        # laws.
275        ref_centroid_values = copy.copy(stage.centroid_values[:])    # Copy
276        ref_vertex_values = copy.copy(stage.vertex_values[:])        # Copy
277
278        # Limit (and invoke balance_deep_and_shallow)
279        domain.tight_slope_limiters = 1
280        domain.distribute_to_vertices_and_edges()
281
282        # Redo derived quantities
283        depth = stage - elevation
284        u = xmomentum/depth
285        v = ymomentum/depth
286
287        # Assert that all vertex velocities stay within one
288        # order of magnitude of centroid velocities.
289        assert num.alltrue(num.absolute(u.vertex_values[1,:]) <=
290                           num.absolute(u.centroid_values[1])*10)
291        assert num.alltrue(num.absolute(v.vertex_values[1,:]) <=
292                           num.absolute(v.centroid_values[1])*10)
293
294        denom = (depth*g)**0.5
295        Fx = u/denom
296        Fy = v/denom
297
298        # Assert that Froude numbers are less than max value (TBA)
299        # at vertices, edges and centroids.
300        from anuga.config import maximum_froude_number
301
302        assert num.alltrue(num.absolute(Fx.vertex_values[1,:]) <
303                           maximum_froude_number)
304        assert num.alltrue(num.absolute(Fy.vertex_values[1,:]) <
305                           maximum_froude_number)
306
307        # Assert that all vertex quantities have changed
308        for k in range(len(domain)):
309            assert not num.allclose(ref_vertex_values[k,:],
310                                    stage.vertex_values[k,:])
311
312        # Assert that quantities are still conserved
313        for k in range(len(domain)):
314            assert num.allclose(ref_centroid_values[k],
315                                num.sum(stage.vertex_values[k,:])/3)
316
317        return
318
319        qwidth = 12
320        for k in [1]:    # range(len(domain)):
321            print 'Triangle %d (C, V, E)' % k
322
323            print ('stage'.ljust(qwidth), stage.centroid_values[k],
324                   stage.vertex_values[k,:], stage.edge_values[k,:])
325            print ('elevation'.ljust(qwidth), elevation.centroid_values[k],
326                   elevation.vertex_values[k,:], elevation.edge_values[k,:])
327            print ('depth'.ljust(qwidth), depth.centroid_values[k],
328                   depth.vertex_values[k,:], depth.edge_values[k,:])
329            print ('xmomentum'.ljust(qwidth), xmomentum.centroid_values[k],
330                   xmomentum.vertex_values[k,:], xmomentum.edge_values[k,:])
331            print ('ymomentum'.ljust(qwidth), ymomentum.centroid_values[k],
332                   ymomentum.vertex_values[k,:], ymomentum.edge_values[k,:])
333            print ('u'.ljust(qwidth),u.centroid_values[k],
334                   u.vertex_values[k,:], u.edge_values[k,:])
335            print ('v'.ljust(qwidth), v.centroid_values[k],
336                   v.vertex_values[k,:], v.edge_values[k,:])
337            print ('Fx'.ljust(qwidth), Fx.centroid_values[k],
338                   Fx.vertex_values[k,:], Fx.edge_values[k,:])
339            print ('Fy'.ljust(qwidth), Fy.centroid_values[k],
340                   Fy.vertex_values[k,:], Fy.edge_values[k,:])
341
342
343#################################################################################
344
345if __name__ == "__main__":
346    suite = unittest.makeSuite(Test_swb_clean, 'test')
347    runner = unittest.TextTestRunner(verbosity=1)
348    runner.run(suite)
Note: See TracBrowser for help on using the repository browser.