source: anuga_core/source/anuga/shallow_water_balanced/test_swb_cross_sections.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: 10.6 KB
RevLine 
[7559]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    def test_get_flow_through_cross_section_with_geo(self):
36        """test_get_flow_through_cross_section(self):
37
38        Test that the total flow through a cross section can be
39        correctly obtained at run-time from the ANUGA domain.
40
41        This test creates a flat bed with a known flow through it and tests
42        that the function correctly returns the expected flow.
43
44        The specifics are
45        e = -1 m
46        u = 2 m/s
47        h = 2 m
48        w = 3 m (width of channel)
49
50        q = u*h*w = 12 m^3/s
51
52        This run tries it with georeferencing and with elevation = -1
53        """
54
55        import time, os
56        from Scientific.IO.NetCDF import NetCDFFile
57        from mesh_factory import rectangular
58
59        # Create basic mesh (20m x 3m)
60        width = 3
61        length = 20
62        t_end = 1
63        points, vertices, boundary = rectangular(length, width, length, width)
64
65        # Create shallow water domain
66        domain = Domain(points, vertices, boundary,
67                        geo_reference=Geo_reference(56, 308500, 6189000))
68
69        domain.default_order = 2
70        domain.set_quantities_to_be_stored(None)
71
72        e = -1.0
73        w = 1.0
74        h = w-e
75        u = 2.0
76        uh = u*h
77
78        Br = Reflective_boundary(domain)     # Side walls
79        Bd = Dirichlet_boundary([w, uh, 0])  # 2 m/s across the 3 m inlet:
80
81
82        # Initial conditions
83        domain.set_quantity('elevation', e)
84        domain.set_quantity('stage', w)
85        domain.set_quantity('xmomentum', uh)
86        domain.set_boundary({'left': Bd, 'right': Bd, 'top': Br, 'bottom': Br})
87
88        # Interpolation points down the middle
89        I = [[0, width/2.],
90             [length/2., width/2.],
91             [length, width/2.]]
92        interpolation_points = domain.geo_reference.get_absolute(I)
93
94        # Shortcuts to quantites
95        stage = domain.get_quantity('stage')
96        xmomentum = domain.get_quantity('xmomentum')
97        ymomentum = domain.get_quantity('ymomentum')
98
99        for t in domain.evolve(yieldstep=0.1, finaltime=t_end):
100            # Check that quantities are they should be in the interior
101            w_t = stage.get_values(interpolation_points)
102            uh_t = xmomentum.get_values(interpolation_points)
103            vh_t = ymomentum.get_values(interpolation_points)
104
105            assert num.allclose(w_t, w)
106            assert num.allclose(uh_t, uh)
107            assert num.allclose(vh_t, 0.0, atol=1.0e-6)
108
109            # Check flows through the middle
110            for i in range(5):
111                x = length/2. + i*0.23674563    # Arbitrary
112                cross_section = [[x, 0], [x, width]]
113
114                cross_section = domain.geo_reference.get_absolute(cross_section)
115                Q = domain.get_flow_through_cross_section(cross_section,
116                                                          verbose=False)
117
118                assert num.allclose(Q, uh*width)
119
120    def test_get_energy_through_cross_section_with_geo(self):
121        """test_get_energy_through_cross_section(self):
122
123        Test that the total and specific energy through a cross section can be
124        correctly obtained at run-time from the ANUGA domain.
125
126        This test creates a flat bed with a known flow through it and tests
127        that the function correctly returns the expected energies.
128
129        The specifics are
130        e = -1 m
131        u = 2 m/s
132        h = 2 m
133        w = 3 m (width of channel)
134
135        q = u*h*w = 12 m^3/s
136
137        This run tries it with georeferencing and with elevation = -1
138        """
139
140        import time, os
141        from Scientific.IO.NetCDF import NetCDFFile
142        from mesh_factory import rectangular
143
144        # Create basic mesh (20m x 3m)
145        width = 3
146        length = 20
147        t_end = 1
148        points, vertices, boundary = rectangular(length, width, length, width)
149
150        # Create shallow water domain
151        domain = Domain(points, vertices, boundary,
152                        geo_reference=Geo_reference(56, 308500, 6189000))
153
154        domain.default_order = 2
155        domain.set_quantities_to_be_stored(None)
156
157        e = -1.0
158        w = 1.0
159        h = w-e
160        u = 2.0
161        uh = u*h
162
163        Br = Reflective_boundary(domain)       # Side walls
164        Bd = Dirichlet_boundary([w, uh, 0])    # 2 m/s across the 3 m inlet:
165
166        # Initial conditions
167        domain.set_quantity('elevation', e)
168        domain.set_quantity('stage', w)
169        domain.set_quantity('xmomentum', uh)
170        domain.set_boundary({'left': Bd, 'right': Bd, 'top': Br, 'bottom': Br})
171
172        # Interpolation points down the middle
173        I = [[0, width/2.],
174             [length/2., width/2.],
175             [length, width/2.]]
176        interpolation_points = domain.geo_reference.get_absolute(I)
177
178        # Shortcuts to quantites
179        stage = domain.get_quantity('stage')
180        xmomentum = domain.get_quantity('xmomentum')
181        ymomentum = domain.get_quantity('ymomentum')
182
183        for t in domain.evolve(yieldstep=0.1, finaltime=t_end):
184            # Check that quantities are they should be in the interior
185            w_t = stage.get_values(interpolation_points)
186            uh_t = xmomentum.get_values(interpolation_points)
187            vh_t = ymomentum.get_values(interpolation_points)
188
189            assert num.allclose(w_t, w)
190            assert num.allclose(uh_t, uh)
191            assert num.allclose(vh_t, 0.0, atol=1.0e-6)
192
193            # Check energies through the middle
194            for i in range(5):
195                x = length/2. + i*0.23674563    # Arbitrary
196                cross_section = [[x, 0], [x, width]]
197
198                cross_section = domain.geo_reference.get_absolute(cross_section)
199                Es = domain.get_energy_through_cross_section(cross_section,
200                                                             kind='specific',
201                                                             verbose=False)
202
203                assert num.allclose(Es, h + 0.5*u*u/g)
204
205                Et = domain.get_energy_through_cross_section(cross_section,
206                                                             kind='total',
207                                                             verbose=False)
208                assert num.allclose(Et, w + 0.5*u*u/g)
209
210
211    def test_cross_section_class(self):
212        """test_cross_section_class(self):
213
214        Test that the total and specific energy through a cross section can be
215        correctly obtained at run-time from the ANUGA cross section class.
216
217        This test creates a flat bed with a known flow through it, creates a cross
218        section and tests that the correct flow and energies are calculated
219
220        The specifics are
221        e = -1 m
222        u = 2 m/s
223        h = 2 m
224        w = 3 m (width of channel)
225
226        q = u*h*w = 12 m^3/s
227
228        This run tries it with georeferencing and with elevation = -1
229        """
230
231        import time, os
232        from Scientific.IO.NetCDF import NetCDFFile
233        from mesh_factory import rectangular
234
235        # Create basic mesh (20m x 3m)
236        width = 3
237        length = 20
238        t_end = 1
239        points, vertices, boundary = rectangular(length, width, length, width)
240
241        # Create shallow water domain
242        domain = Domain(points, vertices, boundary,
243                        geo_reference=Geo_reference(56, 308500, 6189000))
244
245        domain.default_order = 2
246        domain.set_quantities_to_be_stored(None)
247
248        e = -1.0
249        w = 1.0
250        h = w-e
251        u = 2.0
252        uh = u*h
253
254        Br = Reflective_boundary(domain)       # Side walls
255        Bd = Dirichlet_boundary([w, uh, 0])    # 2 m/s across the 3 m inlet:
256
257        # Initial conditions
258        domain.set_quantity('elevation', e)
259        domain.set_quantity('stage', w)
260        domain.set_quantity('xmomentum', uh)
261        domain.set_boundary({'left': Bd, 'right': Bd, 'top': Br, 'bottom': Br})
262
263        # Interpolation points down the middle
264        I = [[0, width/2.],
265             [length/2., width/2.],
266             [length, width/2.]]
267        interpolation_points = domain.geo_reference.get_absolute(I)
268
269        # Shortcuts to quantites
270        stage = domain.get_quantity('stage')
271        xmomentum = domain.get_quantity('xmomentum')
272        ymomentum = domain.get_quantity('ymomentum')
273
274
275        # Create some cross sections
276        cross_sections = []
277        for i in range(5):
278            x = length/2. + i*0.23674563    # Arbitrary
279            polyline = [[x, 0], [x, width]]
280
281            polyline = domain.geo_reference.get_absolute(polyline)
282       
283            cross_sections.append(Cross_section(domain,polyline))
284
285
286 
287        for t in domain.evolve(yieldstep=0.1, finaltime=t_end):
288            # Check that quantities are they should be in the interior
289            w_t = stage.get_values(interpolation_points)
290            uh_t = xmomentum.get_values(interpolation_points)
291            vh_t = ymomentum.get_values(interpolation_points)
292
293            assert num.allclose(w_t, w)
294            assert num.allclose(uh_t, uh)
295            assert num.allclose(vh_t, 0.0, atol=1.0e-6)
296
297
298            # Check flows and energies through the middle
299            for cross_section in cross_sections:
300               
301                Q = cross_section.get_flow_through_cross_section()
302
303                assert num.allclose(Q, uh*width)
304
305                Es = cross_section.get_energy_through_cross_section(kind='specific')
306
307                assert num.allclose(Es, h + 0.5*u*u/g)
308
309                Et = cross_section.get_energy_through_cross_section(kind='total')
310               
311                assert num.allclose(Et, w + 0.5*u*u/g)
312
313
314
315#################################################################################
316
317if __name__ == "__main__":
318    suite = unittest.makeSuite(Test_swb_clean, 'test')
319    runner = unittest.TextTestRunner(verbosity=1)
320    runner.run(suite)
Note: See TracBrowser for help on using the repository browser.