source: trunk/anuga_work/development/2010-projects/anuga_1d/sww/test_sww_domain.py @ 7884

Last change on this file since 7884 was 7884, checked in by steve, 12 years ago

Moving 2010 project

File size: 6.2 KB
Line 
1#!/usr/bin/env python
2
3import unittest
4from math import sqrt, pi
5
6
7from sww_domain import *
8import sww_python
9
10
11import numpy
12
13
14class Test_Shallow_Water(unittest.TestCase):
15    def setUp(self):
16        self.points = [0.0, 1.0, 2.0, 3.0]
17        self.vertex_values = [[1.0,2.0],[4.0,5.0],[-1.0,2.0]]
18        self.points2 = [0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0]
19       
20    def tearDown(self):
21        pass
22        #print "  Tearing down"
23
24
25    def test_creation(self):
26        domain = Domain(self.points)
27        assert numpy.allclose(domain.centroids, [0.5, 1.5, 2.5])
28
29    def test_compute_fluxes(self):
30        """
31        Compare shallow_water_domain flux calculation against a previous
32        Python implementation (defined in this file)
33        """
34        domain = Domain(self.points)
35        domain.set_quantity('stage',2.0)
36        domain.set_boundary({'exterior' : Reflective_boundary(domain)})
37
38        stage_ud, xmom_ud = sww_python.compute_fluxes(domain)
39
40        domain.distribute_to_vertices_and_edges()
41        domain.compute_fluxes()
42
43        assert numpy.allclose( domain.quantities['stage'].explicit_update, stage_ud )
44        assert numpy.allclose( domain.quantities['xmomentum'].explicit_update, xmom_ud )
45
46
47    def test_sww_python_flux_function(self):
48        normal = 1.0
49        ql = numpy.array([1.0, 2.0],numpy.float)
50        qr = numpy.array([1.0, 2.0],numpy.float)
51        zl = 0.0
52        zr = 0.0
53
54        #This assumes h0 = 1.0e-3!!
55        edgeflux, maxspeed = sww_python.flux_function(normal, ql,qr,zl,zr)
56        #print maxspeed
57        #print edgeflux
58       
59        assert numpy.allclose([2.0, 8.9], edgeflux, rtol=1.0e-005)
60        assert numpy.allclose(5.1305, maxspeed, rtol=1.0e-005)
61
62        normal = -1.0
63        ql = numpy.array([1.0, 2.0],numpy.float)
64        qr = numpy.array([1.0, 2.0],numpy.float)
65        zl = 0.0
66        zr = 0.0
67
68        edgeflux, maxspeed = sww_python.flux_function(normal, ql,qr,zl,zr)
69
70
71        #print maxspeed
72        #print edgeflux       
73       
74        assert numpy.allclose([-2.0, -8.9], edgeflux, rtol=1.0e-005)
75        assert numpy.allclose(5.1305, maxspeed, rtol=1.0e-005)
76
77
78
79    def test_gravity_for_loops(self):
80        """
81        Compare shallow_water_domain gravity calculation
82        """
83
84        def slope_one(x):
85            return x
86
87        domain = Domain(self.points)
88        domain.set_quantity('stage',4.0)
89        domain.set_quantity('elevation',slope_one)
90        domain.set_boundary({'exterior' : Reflective_boundary(domain)})
91
92        domain.set_spatial_order(1)
93        domain.distribute_to_vertices_and_edges()
94
95        gravity_for_loops(domain)     
96
97        assert numpy.allclose( [-34.3, -24.5, -14.7], domain.quantities['xmomentum'].explicit_update )
98
99
100    def test_gravity(self):
101        """
102        Compare shallow_water_domain gravity calculation
103        """
104
105        def slope_one(x):
106            return x
107
108        domain = Domain(self.points)
109        domain.set_quantity('stage',4.0)
110        domain.set_quantity('elevation',slope_one)
111        domain.set_boundary({'exterior' : Reflective_boundary(domain)})
112
113
114        domain.set_spatial_order(1)
115        domain.distribute_to_vertices_and_edges()
116
117        gravity(domain)
118
119        assert numpy.allclose( [-34.3, -24.5, -14.7], domain.quantities['xmomentum'].explicit_update )
120
121
122    def test_evolve_first_order(self):
123        """
124        Compare still lake solution for various versions of shallow_water_domain
125        """
126       
127        def slope_square(x):
128            return numpy.maximum(4.0-(x-5.0)*(x-5.0), 0.0)
129
130        domain = Domain(self.points2)
131        domain.set_quantity('stage',10.0)
132        domain.set_quantity('elevation',slope_square)
133        domain.set_boundary({'exterior' : Reflective_boundary(domain)})
134
135        domain.default_order = 1
136        domain.set_timestepping_method('euler')
137       
138        yieldstep=1.0
139        finaltime=1.0
140
141        for t in domain.evolve(yieldstep=yieldstep, finaltime=finaltime):
142            pass
143       
144
145        assert numpy.allclose( 10.0*numpy.ones(10), domain.quantities['stage'].centroid_values )
146        assert numpy.allclose( numpy.zeros(10), domain.quantities['xmomentum'].centroid_values )
147
148
149    def test_evolve_euler_second_order_space(self):
150        """
151        Compare still lake solution for various versions of shallow_water_domain
152        """
153
154        def slope_square(x):
155            return numpy.maximum(4.0-(x-5.0)*(x-5.0), 0.0)
156
157        domain = Domain(self.points2)
158        domain.set_quantity('stage',10.0)
159        domain.set_quantity('elevation',slope_square)
160        domain.set_boundary({'exterior' : Reflective_boundary(domain)})
161
162        domain.default_order = 2
163        domain.set_timestepping_method('rk2')
164        yieldstep=1.0
165        finaltime=1.0
166
167        for t in domain.evolve(yieldstep=yieldstep, finaltime=finaltime):
168            pass
169       
170        assert numpy.allclose( 10.0*numpy.ones(10), domain.quantities['stage'].centroid_values )
171        assert numpy.allclose( numpy.zeros(10), domain.quantities['xmomentum'].centroid_values )
172
173    def test_evolve_second_order_space_time(self):
174        """
175        Compare still lake solution for various versions of shallow_water_domain
176        """
177
178        def slope_square(x):
179            return numpy.maximum(4.0-(x-5.0)*(x-5.0), 0.0)
180
181        domain = Domain(self.points2)
182        domain.set_quantity('stage',10.0)
183        domain.set_quantity('elevation',slope_square)
184        domain.set_boundary({'exterior' : Reflective_boundary(domain)})
185
186        domain.default_order = 2
187        domain.set_timestepping_method('rk3')
188
189        yieldstep=1.0
190        finaltime=1.0
191
192        for t in domain.evolve(yieldstep=yieldstep, finaltime=finaltime):
193            pass
194       
195        assert numpy.allclose( 10.0*numpy.ones(10), domain.quantities['stage'].centroid_values )
196        assert numpy.allclose( numpy.zeros(10), domain.quantities['xmomentum'].centroid_values )
197
198
199
200#==============================================================================
201
202if __name__ == "__main__":
203    suite = unittest.makeSuite(Test_Shallow_Water, 'test_gravity')
204    #suite = unittest.makeSuite(Test_Shallow_Water, 'test_evolve_first_order')
205
206
207    runner = unittest.TextTestRunner()
208    runner.run(suite)
Note: See TracBrowser for help on using the repository browser.