source: trunk/anuga_core/source/anuga/shallow_water/test_loadsave.py @ 7841

Last change on this file since 7841 was 7841, checked in by hudson, 14 years ago

Refactorings to allow tests to pass.

File size: 7.1 KB
RevLine 
[7841]1#!/usr/bin/env python
2
3import unittest, os, time
4import os.path
5from math import pi, sqrt
6import tempfile
7
8from Scientific.IO.NetCDF import NetCDFFile
9from anuga.file.sww import extent_sww
10
11from anuga.config import g, epsilon
12from anuga.config import netcdf_mode_r, netcdf_mode_w, netcdf_mode_a
13from anuga.utilities.numerical_tools import mean, ensure_numeric
14from anuga.geometry.polygon import is_inside_polygon
15from anuga.coordinate_transforms.geo_reference import Geo_reference
16from anuga.geospatial_data.geospatial_data import Geospatial_data
17from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular_cross, \
18                                            rectangular
19from anuga.abstract_2d_finite_volumes.quantity import Quantity
20from anuga.shallow_water.forcing import Inflow, Cross_section
21from anuga.geospatial_data.geospatial_data import ensure_geospatial
22
23from anuga.utilities.system_tools import get_pathname_from_package
24
25from anuga.abstract_2d_finite_volumes.generic_boundary_conditions \
26        import Dirichlet_boundary
27from anuga.shallow_water.forcing import Rainfall, Wind_stress
28from anuga.shallow_water.forcing import Inflow, Cross_section
29from anuga.shallow_water.sww_interrogate import get_flow_through_cross_section
30
31from shallow_water_domain import Domain
32
33# boundary functions
34from boundaries import Reflective_boundary, \
35            Field_boundary, Transmissive_momentum_set_stage_boundary, \
36            Transmissive_stage_zero_momentum_boundary
37from anuga.abstract_2d_finite_volumes.generic_boundary_conditions\
38     import Transmissive_boundary, Dirichlet_boundary, \
39            Time_boundary, File_boundary, AWI_boundary
40
41import numpy as num
42from anuga.config import g
43
44# Get gateway to C implementation of flux function for direct testing
45from shallow_water_ext import flux_function_central as flux_function
46from shallow_water_ext import rotate
47
48
49def set_bottom_friction(tag, elements, domain):
50    if tag == "bottom":
51        domain.set_quantity('friction', 0.09, indices = elements)
52
53def set_top_friction(tag, elements, domain):
54    if tag == "top":
55        domain.set_quantity('friction', 1., indices = elements)
56
57
58def set_all_friction(tag, elements, domain):
59    if tag == 'all':
60        new_values = domain.get_quantity('friction').get_values(indices = elements) + 10.0
61
62        domain.set_quantity('friction', new_values, indices = elements)
63
64
65# For test_fitting_using_shallow_water_domain example
66def linear_function(point):
67    point = num.array(point)
68    return point[:,0]+point[:,1]
69
70
71
72def scalar_func(t, x, y):
73    """Function that returns a scalar.
74
75    Used to test error message when numeric array is expected
76    """
77
78    return 17.7
79
80
81
82class Test_LoadSave(unittest.TestCase):
83    def setUp(self):
84        pass
85
86    def tearDown(self):
87        pass
88       
89    def test_get_flow_through_cross_section_with_geo(self):
90        """test_get_flow_through_cross_section(self):
91
92        Test that the total flow through a cross section can be
93        correctly obtained at run-time from the ANUGA domain.
94
95        This test creates a flat bed with a known flow through it and tests
96        that the function correctly returns the expected flow.
97
98        The specifics are
99        e = -1 m
100        u = 2 m/s
101        h = 2 m
102        w = 3 m (width of channel)
103
104        q = u*h*w = 12 m^3/s
105
106        This run tries it with georeferencing and with elevation = -1
107        """
108
109        # Create basic mesh (20m x 3m)
110        width = 3
111        length = 20
112        t_end = 1
113        points, vertices, boundary = rectangular(length, width, length, width)
114
115        # Create shallow water domain
116        domain = Domain(points, vertices, boundary,
117                        geo_reference=Geo_reference(56, 308500, 6189000))
118
119        domain.default_order = 2
120        domain.set_quantities_to_be_stored(None)
121
122        e = -1.0
123        w = 1.0
124        h = w-e
125        u = 2.0
126        uh = u*h
127
128        Br = Reflective_boundary(domain)     # Side walls
129        Bd = Dirichlet_boundary([w, uh, 0])  # 2 m/s across the 3 m inlet:
130
131
132        # Initial conditions
133        domain.set_quantity('elevation', e)
134        domain.set_quantity('stage', w)
135        domain.set_quantity('xmomentum', uh)
136        domain.set_boundary({'left': Bd, 'right': Bd, 'top': Br, 'bottom': Br})
137
138        # Interpolation points down the middle
139        I = [[0, width/2.],
140             [length/2., width/2.],
141             [length, width/2.]]
142        interpolation_points = domain.geo_reference.get_absolute(I)
143
144        for t in domain.evolve(yieldstep=0.1, finaltime=0.5):
145            # Shortcuts to quantites
146            stage = domain.get_quantity('stage')
147            xmomentum = domain.get_quantity('xmomentum')
148            ymomentum = domain.get_quantity('ymomentum')
149
150            # Check that quantities are they should be in the interior
151            w_t = stage.get_values(interpolation_points)
152            uh_t = xmomentum.get_values(interpolation_points)
153            vh_t = ymomentum.get_values(interpolation_points)
154
155            assert num.allclose(w_t, w)
156            assert num.allclose(uh_t, uh)
157            assert num.allclose(vh_t, 0.0, atol=1.0e-6)
158
159            # Check flows through the middle
160            for i in range(5):
161                x = length/2. + i*0.23674563    # Arbitrary
162                cross_section = [[x, 0], [x, width]]
163
164                cross_section = domain.geo_reference.get_absolute(cross_section)
165                Q = domain.get_flow_through_cross_section(cross_section,
166                                                          verbose=False)
167
168                assert num.allclose(Q, uh*width)
169
170        import cPickle       
171        cPickle.dump(domain, open('domain_pickle.txt', 'w'))
172        domain_restored = cPickle.load(open('domain_pickle.txt'))
173
174        for t in domain_restored.evolve(yieldstep=0.1, finaltime=1.0):
175            # Shortcuts to quantites
176            stage = domain_restored.get_quantity('stage')
177            xmomentum = domain_restored.get_quantity('xmomentum')
178            ymomentum = domain_restored.get_quantity('ymomentum')
179           
180            # Check that quantities are they should be in the interior
181            w_t = stage.get_values(interpolation_points)
182            uh_t = xmomentum.get_values(interpolation_points)
183            vh_t = ymomentum.get_values(interpolation_points)
184
185            assert num.allclose(w_t, w)
186            assert num.allclose(uh_t, uh)
187            assert num.allclose(vh_t, 0.0, atol=1.0e-6)
188
189            # Check flows through the middle
190            for i in range(5):
191                x = length/2. + i*0.23674563    # Arbitrary
192                cross_section = [[x, 0], [x, width]]
193
194                cross_section = domain_restored.geo_reference.get_absolute(cross_section)
195                Q = domain_restored.get_flow_through_cross_section(cross_section,
196                                                          verbose=False)
197
198                assert num.allclose(Q, uh*width)
199               
200
201
202#-------------------------------------------------------------
203
204if __name__ == "__main__":
205    suite = unittest.makeSuite(Test_LoadSave, 'test')
206    runner = unittest.TextTestRunner() #verbosity=2)
207    runner.run(suite)   
Note: See TracBrowser for help on using the repository browser.