source: anuga_work/development/pipeflow/test_sww_domain.py @ 7824

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

test_sww_vel_domain is nowworking with numpy

File size: 5.7 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(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        gravity(domain)
93       
94        #print domain.quantities['stage'].vertex_values
95        #print domain.quantities['elevation'].vertex_values
96        #print domain.quantities['xmomentum'].explicit_update       
97
98        assert numpy.allclose( [-34.3, -24.5, -14.7], domain.quantities['xmomentum'].explicit_update )
99
100
101    def test_evolve_first_order(self):
102        """
103        Compare still lake solution for various versions of shallow_water_domain
104        """
105       
106        def slope_square(x):
107            return numpy.maximum(4.0-(x-5.0)*(x-5.0), 0.0)
108
109        domain = Domain(self.points2)
110        domain.set_quantity('stage',10.0)
111        domain.set_quantity('elevation',slope_square)
112        domain.set_boundary({'exterior' : Reflective_boundary(domain)})
113
114        domain.default_order = 1
115        domain.set_timestepping_method('euler')
116       
117        yieldstep=1.0
118        finaltime=1.0
119
120        for t in domain.evolve(yieldstep=yieldstep, finaltime=finaltime):
121            pass
122       
123
124        assert numpy.allclose( 10.0*numpy.ones(10), domain.quantities['stage'].centroid_values )
125        assert numpy.allclose( numpy.zeros(10), domain.quantities['xmomentum'].centroid_values )
126
127
128    def test_evolve_euler_second_order_space(self):
129        """
130        Compare still lake solution for various versions of shallow_water_domain
131        """
132
133        def slope_square(x):
134            return numpy.maximum(4.0-(x-5.0)*(x-5.0), 0.0)
135
136        domain = Domain(self.points2)
137        domain.set_quantity('stage',10.0)
138        domain.set_quantity('elevation',slope_square)
139        domain.set_boundary({'exterior' : Reflective_boundary(domain)})
140
141        domain.default_order = 2
142        domain.set_timestepping_method('rk2')
143        yieldstep=1.0
144        finaltime=1.0
145
146        for t in domain.evolve(yieldstep=yieldstep, finaltime=finaltime):
147            pass
148       
149        assert numpy.allclose( 10.0*numpy.ones(10), domain.quantities['stage'].centroid_values )
150        assert numpy.allclose( numpy.zeros(10), domain.quantities['xmomentum'].centroid_values )
151
152    def test_evolve_second_order_space_time(self):
153        """
154        Compare still lake solution for various versions of shallow_water_domain
155        """
156
157        def slope_square(x):
158            return numpy.maximum(4.0-(x-5.0)*(x-5.0), 0.0)
159
160        domain = Domain(self.points2)
161        domain.set_quantity('stage',10.0)
162        domain.set_quantity('elevation',slope_square)
163        domain.set_boundary({'exterior' : Reflective_boundary(domain)})
164
165        domain.default_order = 2
166        domain.set_timestepping_method('rk3')
167
168        yieldstep=1.0
169        finaltime=1.0
170
171        for t in domain.evolve(yieldstep=yieldstep, finaltime=finaltime):
172            pass
173       
174        assert numpy.allclose( 10.0*numpy.ones(10), domain.quantities['stage'].centroid_values )
175        assert numpy.allclose( numpy.zeros(10), domain.quantities['xmomentum'].centroid_values )
176
177
178
179#==============================================================================
180
181if __name__ == "__main__":
182    suite = unittest.makeSuite(Test_Shallow_Water, 'test')
183    #suite = unittest.makeSuite(Test_Shallow_Water, 'test_evolve_first_order')
184
185
186    runner = unittest.TextTestRunner()
187    runner.run(suite)
Note: See TracBrowser for help on using the repository browser.