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

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

Broke up test_swb_domain so that we can concentrate on new functionality

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