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

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

Adding old files as well as workinng files from anuga_1d

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        print domain.quantities['stage'].explicit_update
44        print domain.quantities['xmomentum'].explicit_update
45
46        print stage_ud
47        print xmom_ud
48
49        assert numpy.allclose( domain.quantities['stage'].explicit_update, stage_ud )
50        assert numpy.allclose( domain.quantities['xmomentum'].explicit_update, xmom_ud )
51
52
53    def test_sww_python_flux_function(self):
54        normal = 1.0
55        ql = numpy.array([1.0, 2.0],numpy.float)
56        qr = numpy.array([1.0, 2.0],numpy.float)
57        zl = 0.0
58        zr = 0.0
59
60        #This assumes h0 = 1.0e-3!!
61        edgeflux, maxspeed = sww_python.flux_function(normal, ql,qr,zl,zr)
62        #print maxspeed
63        #print edgeflux
64       
65        assert numpy.allclose([2.0, 8.9], edgeflux, rtol=1.0e-005)
66        assert numpy.allclose(5.1305, maxspeed, rtol=1.0e-005)
67
68        normal = -1.0
69        ql = numpy.array([1.0, 2.0],numpy.float)
70        qr = numpy.array([1.0, 2.0],numpy.float)
71        zl = 0.0
72        zr = 0.0
73
74        edgeflux, maxspeed = sww_python.flux_function(normal, ql,qr,zl,zr)
75
76
77        #print maxspeed
78        #print edgeflux       
79       
80        assert numpy.allclose([-2.0, -8.9], edgeflux, rtol=1.0e-005)
81        assert numpy.allclose(5.1305, maxspeed, rtol=1.0e-005)
82
83
84
85    def test_gravity(self):
86        """
87        Compare shallow_water_domain gravity calculation
88        """
89
90        def slope_one(x):
91            return x
92
93        domain = Domain(self.points)
94        domain.set_quantity('stage',4.0)
95        domain.set_quantity('elevation',slope_one)
96        domain.set_boundary({'exterior' : Reflective_boundary(domain)})
97
98        gravity(domain)
99       
100        #print domain.quantities['stage'].vertex_values
101        #print domain.quantities['elevation'].vertex_values
102        #print domain.quantities['xmomentum'].explicit_update       
103
104        assert numpy.allclose( [-34.3, -24.5, -14.7], domain.quantities['xmomentum'].explicit_update )
105
106
107    def test_evolve_first_order(self):
108        """
109        Compare still lake solution for various versions of shallow_water_domain
110        """
111       
112        def slope_square(x):
113            return numpy.maximum(4.0-(x-5.0)*(x-5.0), 0.0)
114
115        domain = Domain(self.points2)
116        domain.set_quantity('stage',10.0)
117        domain.set_quantity('elevation',slope_square)
118        domain.set_boundary({'exterior' : Reflective_boundary(domain)})
119
120        domain.default_order = 1
121        domain.set_timestepping_method('euler')
122       
123        yieldstep=1.0
124        finaltime=1.0
125
126        for t in domain.evolve(yieldstep=yieldstep, finaltime=finaltime):
127            #pass
128       
129            print domain.quantities['stage'].vertex_values
130            print domain.quantities['elevation'].vertex_values
131            print domain.quantities['xmomentum'].vertex_values
132         
133       
134            print domain.quantities['stage'].centroid_values 
135            print domain.quantities['elevation'].centroid_values
136            print domain.quantities['xmomentum'].centroid_values   
137
138        assert numpy.allclose( 10.0*numpy.ones(10), domain.quantities['stage'].centroid_values )
139        assert numpy.allclose( numpy.zeros(10), domain.quantities['xmomentum'].centroid_values )
140
141
142    def test_evolve_euler_second_order_space(self):
143        """
144        Compare still lake solution for various versions of shallow_water_domain
145        """
146
147        def slope_square(x):
148            return numpy.maximum(4.0-(x-5.0)*(x-5.0), 0.0)
149
150        domain = Domain(self.points2)
151        domain.set_quantity('stage',10.0)
152        domain.set_quantity('elevation',slope_square)
153        domain.set_boundary({'exterior' : Reflective_boundary(domain)})
154
155        domain.default_order = 2
156        domain.set_timestepping_method('rk2')
157        yieldstep=1.0
158        finaltime=1.0
159
160        for t in domain.evolve(yieldstep=yieldstep, finaltime=finaltime):
161            pass
162       
163        assert numpy.allclose( 10.0*ones(10), domain.quantities['stage'].centroid_values )
164        assert numpy.allclose( zeros(10), domain.quantities['xmomentum'].centroid_values )
165
166    def test_evolve_second_order_space_time(self):
167        """
168        Compare still lake solution for various versions of shallow_water_domain
169        """
170
171        def slope_square(x):
172            return numpy.maximum(4.0-(x-5.0)*(x-5.0), 0.0)
173
174        domain = Domain(self.points2)
175        domain.set_quantity('stage',10.0)
176        domain.set_quantity('elevation',slope_square)
177        domain.set_boundary({'exterior' : Reflective_boundary(domain)})
178
179        domain.default_order = 2
180        domain.set_timestepping_method('rk3')
181
182        yieldstep=1.0
183        finaltime=1.0
184
185        for t in domain.evolve(yieldstep=yieldstep, finaltime=finaltime):
186            pass
187       
188        assert numpy.allclose( 10.0*ones(10), domain.quantities['stage'].centroid_values )
189        assert numpy.allclose( zeros(10), domain.quantities['xmomentum'].centroid_values )
190
191
192
193#==============================================================================
194
195if __name__ == "__main__":
196    suite = unittest.makeSuite(Test_Shallow_Water, 'test')
197    #suite = unittest.makeSuite(Test_Shallow_Water, 'test_evolve_first_order')
198
199
200    runner = unittest.TextTestRunner()
201    runner.run(suite)
Note: See TracBrowser for help on using the repository browser.