source: development/pyvolution-1d/test_quantity.py @ 2716

Last change on this file since 2716 was 2716, checked in by jakeman, 18 years ago

Adding new files

File size: 12.9 KB
Line 
1#!/usr/bin/env python
2
3import unittest
4from math import sqrt, pi
5
6
7from quantity import *
8#from config import epsilon
9from Numeric import allclose, array, zeros, Float
10
11       
12class TestCase(unittest.TestCase):
13    def setUp(self):
14        from domain import Domain
15       
16        a = 0.0
17        b = 1.0
18        c = 2.0
19        d = 2.5
20        e = 3.1
21        f = 4.0
22
23        self.points          = [a, b, c, d, e, f]
24        self.centroids       = [(a+b)/2,(b+c)/2,(c+d)/2,(d+e)/2,(e+f)/2]
25        self.vertex_values   = [[1.0,2.0],[2.0,3.0],[3.0,4.0],[4.5,5],[5.5,5.6]]
26        self.centroid_values = [[1.5,      2.5,      3.5,      4.75,   5.55]]
27
28        self.domain1 = Domain(self.points[0:2])
29        self.domain5 = Domain(self.points)
30       
31    def tearDown(self):
32        pass
33        #print "  Tearing down"
34
35
36    def test_creation(self):
37       
38        quantity = Quantity(self.domain5, self.vertex_values)
39        assert allclose(quantity.centroid_values, self.centroid_values)
40
41        try:
42            quantity = Quantity()
43        except:
44            pass
45        else:
46            raise 'Should have raised empty quantity exception'
47
48
49        try:
50            quantity = Quantity([1,2])
51        except AssertionError:
52            pass
53        except:
54            raise 'Should have raised "missing domain object" error'       
55       
56
57    def test_creation_zeros(self):
58       
59        quantity = Quantity(self.domain1)
60        assert allclose(quantity.vertex_values, [[0.,0.]])
61
62
63        quantity = Quantity(self.domain5)
64        assert allclose(quantity.vertex_values, [[0.,0.], [0.,0.],
65                                                 [0.,0.], [0.,0.],
66                                                 [0.,0.]])       
67
68       
69    def test_interpolation(self):
70        quantity = Quantity(self.domain1, [[1,2]])
71        assert allclose(quantity.centroid_values, 1.5) #Centroid
72
73
74    def test_interpolation2(self):
75        quantity = Quantity(self.domain5,
76                            [[1,2], [5,5], [0,9], [-6, 3], [3,4]])
77        assert allclose(quantity.centroid_values, [1.5, 5., 4.5, -1.5, 3.5 ]) #Centroid
78
79    def test_set_values(self):
80        quantity = Quantity(self.domain5)
81
82
83        quantity.set_values([[1,2], [5,5], [0,0], [-6, 3], [-2,4]],
84                            location = 'vertices')
85        assert allclose(quantity.vertex_values,
86                        [[1,2], [5,5], [0,0], [-6, 3], [-2,4]])
87        assert allclose(quantity.centroid_values, [1.5, 5., 0., -1.5, 1.0]) #Centroid
88
89        #Test default
90        quantity.set_values([[1,2], [5,5], [0,0], [-6, 3], [-2,4]])
91        assert allclose(quantity.vertex_values,
92                        [[1,2], [5,5], [0,0], [-6, 3], [-2,4]])       
93        assert allclose(quantity.centroid_values, [1.5, 5., 0., -1.5, 1.0]) #Centroid
94
95        #Test centroids
96        quantity.set_values([1,2,3,4,5], location = 'centroids')
97        assert allclose(quantity.centroid_values, [1., 2., 3., 4., 5.]) #Centroid
98
99        #Test exceptions
100        try:
101            quantity.set_values([[1,2], [5,5], [0,0], [-6, 3], [-2,4]],
102                                location = 'bas kamel tuba')
103        except:
104            pass
105
106
107        try:
108            quantity.set_values([[1,2], [0,0]])
109        except AssertionError:
110            pass
111        except:
112            raise 'should have raised AssertionError'
113
114
115    def test_set_values_const(self):
116        quantity = Quantity(self.domain5)
117
118        quantity.set_values(1.0, location = 'vertices')
119        assert allclose(quantity.vertex_values,
120                        [[1,1], [1,1], [1,1], [1,1], [1,1]])
121        assert allclose(quantity.centroid_values, [1, 1, 1, 1, 1]) #Centroid
122
123        quantity.set_values(2.0, location = 'centroids')
124        assert allclose(quantity.centroid_values, [2, 2, 2, 2, 2])
125
126    def test_set_values_func(self):
127        quantity = Quantity(self.domain5)
128
129        def f(x):
130            return x**2
131
132        quantity.set_values(f, location = 'vertices')
133        assert allclose(quantity.vertex_values,
134                        [[0,1], [1,4], [4,6.25], [6.25,9.61], [9.61,16]])       
135        assert allclose(quantity.centroid_values,
136                        [0.5, 2.5, 5.125, 7.93, 12.805 ])
137       
138        quantity.set_values(f, location = 'centroids')
139        assert allclose(quantity.centroid_values,
140                        [0.25, 1.5**2, 2.25**2, 2.8**2, 3.55**2])       
141
142
143    def test_gradient(self):
144       
145        quantity = Conserved_quantity(self.domain5)
146
147        #Set up for a gradient of (3,0) at mid triangle
148        quantity.set_values([2.0, 4.0, 8.0, 2.0, 5.0],
149                            location = 'centroids')
150       
151        #Gradients
152        quantity.compute_gradients()
153        N = quantity.gradients.shape[0]
154       
155        G = quantity.gradients
156        G0 = zeros(N, Float)
157        Q = quantity.centroid_values
158        X = quantity.domain.centroids
159
160        def grad0(x0,x1,x2,q0,q1,q2):
161            return ((q0-q1)/(x0-x1)*(x2-x1) - (q2-q1)/(x2-x1)*(x0-x1))/(x2-x0)
162
163        def grad1(x0,x1,q0,q1):
164            return  (q1 - q0)/(x1 - x0)
165
166        #Check Gradients
167        G0[0] = grad1(X[0],X[1],Q[0],Q[1])
168        for i in range(1,4):
169            G0[i] = grad0(X[i-1],X[i],X[i+1],Q[i-1],Q[i],Q[i+1])
170        G0[4] = grad1(X[3],X[4],Q[3],Q[4])
171
172        assert allclose(G,G0)
173
174        quantity.extrapolate_first_order()
175       
176        assert allclose(quantity.vertex_values, [[2., 2.],
177                                                 [4., 4.],
178                                                 [8., 8.],
179                                                 [2., 2.],
180                                                 [5., 5.]])
181
182        quantity.extrapolate_second_order()
183
184        V = quantity.domain.vertices
185        for k in range(quantity.domain.number_of_elements):
186            G0[k] = G0[k]*(V[k,1]-V[k,0])*0.5
187
188        assert allclose(quantity.vertex_values, [[2.-G0[0], 2.+G0[0]],
189                                                 [4.-G0[1], 4.+G0[1]],
190                                                 [8.-G0[2], 8.+G0[2]],
191                                                 [2.-G0[3], 2.+G0[3]],
192                                                 [5.-G0[4], 5.+G0[4]]])
193
194
195
196    def test_second_order_extrapolation2(self):
197        quantity = Conserved_quantity(self.domain5)       
198
199        #Set up for a gradient of (2), f(x) = 2x
200        quantity.set_values([1., 3., 4.5, 5.6, 7.1],
201                            location = 'centroids')
202       
203        #Gradients
204        quantity.compute_gradients()
205        G = quantity.gradients
206       
207        allclose(G, 2.0)
208
209        quantity.extrapolate_second_order()
210
211        V = quantity.domain.vertices
212        Q = quantity.vertex_values
213       
214        assert allclose(Q[:,0], 2.0*V[:,0])
215        assert allclose(Q[:,1], 2.0*V[:,1])       
216
217
218
219
220
221       
222
223
224##     def test_first_order_extrapolator(self):
225##         quantity = Conserved_quantity(self.mesh4)               
226
227##         #Test centroids
228##         quantity.set_values([1.,2.,3.,4.], location = 'centroids')
229##         assert allclose(quantity.centroid_values, [1, 2, 3, 4]) #Centroid
230
231##         #Extrapolate
232##         quantity.extrapolate_first_order()
233
234##         #Check vertices but not edge values
235##         assert allclose(quantity.vertex_values,
236##                         [[1,1,1], [2,2,2], [3,3,3], [4, 4, 4]])
237
238
239    def test_second_order_limit_extrapolator(self):
240        quantity = Conserved_quantity(self.domain5)                       
241
242        #Set up for a gradient of (3,0) at mid triangle
243        quantity.set_values([2.0, 4.0, 8.0, 2.0, 0.0],
244                            location = 'centroids')
245       
246        quantity.extrapolate_second_order()
247        quantity.limit()
248
249        #Assert that central triangle is limited by neighbours
250        assert quantity.vertex_values[0,0] >= 2.0
251        assert quantity.vertex_values[0,0] <= 4.0
252        assert quantity.vertex_values[0,1] >= 2.0
253        assert quantity.vertex_values[0,1] <= 4.0
254       
255        assert quantity.vertex_values[1,0] >= 2.0
256        assert quantity.vertex_values[1,0] <= 8.0
257        assert quantity.vertex_values[1,1] >= 2.0
258        assert quantity.vertex_values[1,1] <= 8.0
259
260        assert quantity.vertex_values[2,0] >= 2.0
261        assert quantity.vertex_values[2,0] <= 8.0
262        assert quantity.vertex_values[2,1] >= 2.0
263        assert quantity.vertex_values[2,1] <= 8.0
264
265        assert quantity.vertex_values[3,0] >= 0.0
266        assert quantity.vertex_values[3,0] <= 8.0
267        assert quantity.vertex_values[3,1] >= 0.0
268        assert quantity.vertex_values[3,1] <= 8.0
269
270        assert quantity.vertex_values[4,0] >= 0.0
271        assert quantity.vertex_values[4,0] <= 2.0
272        assert quantity.vertex_values[4,1] >= 0.0
273        assert quantity.vertex_values[4,1] <= 2.0
274
275        #Assert that quantities are conserved
276        from Numeric import sum
277        for k in range(quantity.centroid_values.shape[0]):
278            assert allclose (quantity.centroid_values[k],
279                             sum(quantity.vertex_values[k,:])/2)
280       
281       
282       
283
284    def test_limiter(self):
285        quantity = Conserved_quantity(self.domain5)
286
287        #Create a deliberate overshoot (e.g. from gradient computation)
288        quantity.set_values([[3,4], [5,5], [0,0], [-6, 3], [-2,4]],
289                            location = 'vertices')       
290
291        #Limit
292        quantity.limit()
293
294        #Assert that central triangle is limited by neighbours
295        assert quantity.vertex_values[0,1] >= quantity.centroid_values[0]
296        assert quantity.vertex_values[1,0] <= quantity.centroid_values[1]
297       
298        assert quantity.vertex_values[1,1] >= quantity.centroid_values[1]
299
300       
301        #Assert that quantities are conserved
302        from Numeric import sum
303        for k in range(quantity.centroid_values.shape[0]):
304            assert allclose (quantity.centroid_values[k],
305                             sum(quantity.vertex_values[k,:])/2)
306       
307
308
309##     def test_distribute_first_order(self):
310##         quantity = Conserved_quantity(self.mesh4)       
311
312##         #Test centroids
313##         quantity.set_values([1.,2.,3.,4.], location = 'centroids')
314##         assert allclose(quantity.centroid_values, [1, 2, 3, 4]) #Centroid
315
316
317##         #Extrapolate
318##         quantity.extrapolate_first_order()
319
320##         #Interpolate
321##         quantity.interpolate_from_vertices_to_edges()       
322
323##         assert allclose(quantity.vertex_values,
324##                         [[1,1,1], [2,2,2], [3,3,3], [4, 4, 4]])
325##         assert allclose(quantity.edge_values, [[1,1,1], [2,2,2],
326##                                                [3,3,3], [4, 4, 4]])
327
328
329
330    def test_update_explicit(self):
331        quantity = Conserved_quantity(self.domain5)
332
333        #Test centroids
334        quantity.set_values([1.,2.,3.,4.,5.], location = 'centroids')
335        assert allclose(quantity.centroid_values, [1, 2, 3, 4, 5]) #Centroid
336
337        #Set explicit_update
338        quantity.explicit_update = array( [1.,1.,1.,1.,1.] )
339
340        #Update with given timestep
341        quantity.update(0.1)
342
343        x = array([1, 2, 3, 4, 5]) + array( [.1,.1,.1,.1,.1] )
344        assert allclose( quantity.centroid_values, x)
345
346    def test_update_semi_implicit(self):
347        quantity = Conserved_quantity(self.domain5)
348
349        #Test centroids
350        quantity.set_values([1.,2.,3.,4.,5.], location = 'centroids')
351        assert allclose(quantity.centroid_values, [1, 2, 3, 4, 5]) #Centroid
352
353        #Set semi implicit update
354        quantity.semi_implicit_update = array( [1.,1.,1.,1.,1.] )
355
356        #Update with given timestep
357        quantity.update(0.1)
358
359        x = array([1, 2, 3, 4, 5])/array( [.9,.9,.9,.9,.9] )
360        assert allclose( quantity.centroid_values, x)
361
362    def test_both_updates(self):
363        quantity = Conserved_quantity(self.domain5)
364
365        #Test centroids
366        quantity.set_values([1.,2.,3.,4.,5.], location = 'centroids')
367        assert allclose(quantity.centroid_values, [1, 2, 3, 4, 5]) #Centroid
368
369        #Set explicit_update
370        quantity.explicit_update = array( [4.,3.,2.,1.,0.0] )
371       
372        #Set semi implicit update
373        quantity.semi_implicit_update = array( [1.,1.,1.,1.,1.] )
374
375        #Update with given timestep
376        quantity.update(0.1)
377
378        x = array([1, 2, 3, 4, 5]) + array( [.4,.3,.2,.1,0.0] )
379        x /= array( [.9,.9,.9,.9,.9] )
380        assert allclose( quantity.centroid_values, x)       
381
382       
383       
384#-------------------------------------------------------------
385if __name__ == "__main__":
386    suite = unittest.makeSuite(TestCase,'test')
387    runner = unittest.TextTestRunner(verbosity=2)
388    runner.run(suite)
389
390
Note: See TracBrowser for help on using the repository browser.