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

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

Changing name of 1d projects so that it will be easy to moveto the trunk

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.