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

Last change on this file since 7573 was 7573, checked in by steve, 13 years ago

Committing a version of shallow_water_balanced which passes it unit test
using a version of edge limiting which doesn't limit boundary edges. THis
is useful to allow linear functions to be reconstructed.

Had to play with the transmissive BC and use centroid values which is
set via domain set_centroid_transmissive_bc

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