source: inundation/ga/storm_surge/pyvolution-1d/test_quantity.py @ 279

Last change on this file since 279 was 279, checked in by steve, 20 years ago

Working on Quantities
File size: 15.2 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
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
80##     def test_boundary_allocation(self):
81##         quantity = Conserved_quantity(self.mesh4,
82##                             [[1,2,3], [5,5,5], [0,0,9], [-6, 3, 3]])
83
84##         assert quantity.boundary_values.shape[0] == len(self.mesh4.boundary)
85
86
87    def test_set_values(self):
88        quantity = Quantity(self.domain5)
89
90
91        quantity.set_values([[1,2], [5,5], [0,0], [-6, 3], [-2,4]],
92                            location = 'vertices')
93        assert allclose(quantity.vertex_values,
94                        [[1,2], [5,5], [0,0], [-6, 3], [-2,4]])
95        assert allclose(quantity.centroid_values, [1.5, 5., 0., -1.5, 1.0]) #Centroid
96
97        #Test default
98        quantity.set_values([[1,2], [5,5], [0,0], [-6, 3], [-2,4]])
99        assert allclose(quantity.vertex_values,
100                        [[1,2], [5,5], [0,0], [-6, 3], [-2,4]])       
101        assert allclose(quantity.centroid_values, [1.5, 5., 0., -1.5, 1.0]) #Centroid
102
103        #Test centroids
104        quantity.set_values([1,2,3,4,5], location = 'centroids')
105        assert allclose(quantity.centroid_values, [1., 2., 3., 4., 5.]) #Centroid
106
107        #Test exceptions
108        try:
109            quantity.set_values([[1,2], [5,5], [0,0], [-6, 3], [-2,4]],
110                                location = 'bas kamel tuba')
111        except:
112            pass
113
114
115        try:
116            quantity.set_values([[1,2], [0,0]])
117        except AssertionError:
118            pass
119        except:
120            raise 'should have raised AssertionError'
121
122
123
124    def test_set_values_const(self):
125        quantity = Quantity(self.domain5)
126
127        quantity.set_values(1.0, location = 'vertices')
128        assert allclose(quantity.vertex_values,
129                        [[1,1], [1,1], [1,1], [1,1], [1,1]])
130        assert allclose(quantity.centroid_values, [1, 1, 1, 1, 1]) #Centroid
131
132
133        quantity.set_values(2.0, location = 'centroids')
134        assert allclose(quantity.centroid_values, [2, 2, 2, 2, 2])
135
136
137    def test_set_values_func(self):
138        quantity = Quantity(self.domain5)
139
140        def f(x):
141            return x**2
142
143        quantity.set_values(f, location = 'vertices')
144        assert allclose(quantity.vertex_values,
145                        [[0,1], [1,4], [4,6.25], [6.25,9.61], [9.61,16]])       
146        assert allclose(quantity.centroid_values,
147                        [0.5, 2.5, 5.125, 7.93, 12.805 ])
148       
149        quantity.set_values(f, location = 'centroids')
150        assert allclose(quantity.centroid_values,
151                        [0.25, 1.5**2, 2.25**2, 2.8**2, 3.55**2])       
152
153
154##     def test_gradient(self):
155##         quantity = Conserved_quantity(self.mesh4)
156
157##         #Set up for a gradient of (3,0) at mid triangle
158##         quantity.set_values([2.0, 4.0, 8.0, 2.0],
159##                             location = 'centroids')
160       
161##         #Gradients
162##         a, b = quantity.compute_gradients()
163
164
165##         #gradient bewteen t0 and t1 is undefined as det==0
166##         assert a[0] == 0.0
167##         assert b[0] == 0.0
168##         #The others are OK
169##         for i in range(1,4):
170##             assert a[i] == 3.0
171##             assert b[i] == 0.0
172
173
174##         quantity.extrapolate_second_order()
175       
176##         assert allclose(quantity.vertex_values, [[2., 2.,  2.],
177##                                                  [0., 6.,  6.],
178##                                                  [6., 6., 12.],
179##                                                  [0., 0.,  6.]])
180
181
182
183##     def test_second_order_extrapolation2(self):
184##         quantity = Conserved_quantity(self.mesh4)       
185
186##         #Set up for a gradient of (3,1), f(x) = 3x+y
187##         quantity.set_values([2.0+2.0/3, 4.0+4.0/3, 8.0+2.0/3, 2.0+8.0/3],
188##                             location = 'centroids')
189       
190##         #Gradients
191##         a, b = quantity.compute_gradients()
192
193##         #gradient bewteen t0 and t1 is undefined as det==0
194##         assert a[0] == 0.0
195##         assert b[0] == 0.0
196##         #The others are OK
197##      163.968025   for i in range(1,4):
198##             assert allclose(a[i], 3.0)
199##             assert allclose(b[i], 1.0)
200
201
202##         quantity.extrapolate_second_order()
203
204##         assert allclose(quantity.vertex_values[1,0], 2.0)
205##         assert allclose(quantity.vertex_values[1,1], 6.0)       
206##         assert allclose(quantity.vertex_values[1,2], 8.0)
207
208
209
210## #     def test_limiter(self):
211
212## #         initialise_consecutive_datastructure(points=6+4, elements=4)         
213       
214## #         a = Point (0.0, 0.0)
215## #         b = Point (0.0, 2.0)
216## #         c = Point (2.0, 0.0)
217## #         d = Point (0.0, 4.0)
218## #         e = Point (2.0, 2.0)
219## #         f = Point (4.0, 0.0)
220
221## #         #Set up for a gradient of (3,1), f(x) = 3x+y
222## #         v1 = Volume(b,a,c,array([0.0,0,0]))       
223## #         v2 = Volume(b,c,e,array([1.0,0,0]))
224## #         v3 = Volume(e,c,f,array([10.0,0,0]))
225## #         v4 = Volume(d,b,e,array([0.0,0,0]))
226
227## #         #Setup neighbour structure
228## #         domain = Domain([v1,v2,v3,v4])
229## #         domain.precompute()
230       
231## #         #Lets's check first order first, hey
232## #    domain.order = 1
233## #    domain.limiter = None
234## #         distribute_to_vertices_and_edges(domain)
235## #         assert allclose(v2.conserved_quantities_vertex0,
236## #                         v2.conserved_quantities_centroid)
237## #         assert allclose(v2.conserved_quantities_vertex1,
238## #                         v2.conserved_quantities_centroid)
239## #         assert allclose(v2.conserved_quantities_vertex2,
240## #                         v2.conserved_quantities_centroid)
241
242
243## #         #Gradient of fitted pwl surface
244## #         a, b = compute_gradient(v2.id)
245       
246
247## #         assert abs(a[0] - 5.0) < epsilon
248## #         assert abs(b[0]) < epsilon       
249## #         #assert qminr[0] == 0.0
250## #         #assert qmaxr[0] == 10.0
251       
252## #         #And now for the second order stuff       
253## #         # - the full second order extrapolation
254## #    domain.order = 2       
255## #         distribute_to_vertices_and_edges(domain)
256
257
258## #         qmin = qmax = v2.conserved_quantities_centroid
259
260## #         qmin = minimum(qmin, v1.conserved_quantities_centroid)
261## #         qmax = maximum(qmax, v1.conserved_quantities_centroid)
262
263## #         qmin = minimum(qmin, v3.conserved_quantities_centroid)
264## #         qmax = maximum(qmax, v3.conserved_quantities_centroid)
265
266## #         qmin = minimum(qmin, v4.conserved_quantities_centroid)
267## #         qmax = maximum(qmax, v4.conserved_quantities_centroid)                   
268## #         #assert qminr == qmin
269## #         #assert qmaxr == qmax
270
271## #         assert v2.conserved_quantities_vertex0 <= qmax
272## #         assert v2.conserved_quantities_vertex0 >= qmin
273## #         assert v2.conserved_quantities_vertex1 <= qmax
274## #         assert v2.conserved_quantities_vertex1 >= qmin
275## #         assert v2.conserved_quantities_vertex2 <= qmax
276## #         assert v2.conserved_quantities_vertex2 >= qmin                   
277163.968025
278
279## #         #Check that volume has been preserved   
280
281## #         q = v2.conserved_quantities_centroid[0]
282## #         w = (v2.conserved_quantities_vertex0[0] +
283## #              v2.conserved_quantities_vertex1[0] +
284## #              v2.conserved_quantities_vertex2[0])/3
285
286## #         assert allclose(q, w)
287
288
289       
290
291
292##     def test_first_order_extrapolator(self):
293##         quantity = Conserved_quantity(self.mesh4)               
294
295##         #Test centroids
296##         quantity.set_values([1.,2.,3.,4.], location = 'centroids')
297##         assert allclose(quantity.centroid_values, [1, 2, 3, 4]) #Centroid
298
299##         #Extrapolate
300##         quantity.extrapolate_first_order()
301
302##         #Check vertices but not edge values
303##         assert allclose(quantity.vertex_values,
304##                         [[1,1,1], [2,2,2], [3,3,3], [4, 4, 4]])
305
306
307##     def test_second_order_extrapolator(self):
308##         quantity = Conserved_quantity(self.mesh4)                       
309
310##         #Set up for a gradient of (3,0) at mid triangle
311##         quantity.set_values([2.0, 4.0, 8.0, 2.0],
312##                             location = 'centroids')
313       
314
315
316##         quantity.extrapolate_second_order()
317##         quantity.limit()
318
319
320##         #Assert that central triangle is limited by neighbours
321##         assert quantity.vertex_values[1,0] >= quantity.vertex_values[0,0]
322##         assert quantity.vertex_values[1,0] >= quantity.vertex_values[3,1]
323       
324##         assert quantity.vertex_values[1,1] <= quantity.vertex_values[2,1]
325##         assert quantity.vertex_values[1,1] >= quantity.vertex_values[0,2]
326       
327##         assert quantity.vertex_values[1,2] <= quantity.vertex_values[2,0]
328##         assert quantity.vertex_values[1,2] >= quantity.vertex_values[3,1]
329
330       
331##         #Assert that quantities are conserved
332##         from Numeric import sum
333##         for k in range(quantity.centroid_values.shape[0]):
334##             assert allclose (quantity.centroid_values[k],
335##                              sum(quantity.vertex_values[k,:])/3)
336       
337       
338
339       
340
341##     def test_limiter(self):
342##         quantity = Conserved_quantity(self.mesh4)
343
344##         #Create a deliberate overshoot (e.g. from gradient computation)
345##         quantity.set_values([[3,0,3], [2,2,6], [5,3,8], [8,3,5]])
346
347
348##         #Limit
349##         quantity.limit()
350
351##         #Assert that central triangle is limited by neighbours
352##         assert quantity.vertex_values[1,0] >= quantity.vertex_values[0,0]
353##         assert quantity.vertex_values[1,0] <= quantity.vertex_values[3,1]
354       
355##         assert quantity.vertex_values[1,1] <= quantity.vertex_values[2,1]
356##         assert quantity.vertex_values[1,1] >= quantity.vertex_values[0,2]
357       
358##         assert quantity.vertex_values[1,2] <= quantity.vertex_values[2,0]
359##         assert quantity.vertex_values[1,2] <= quantity.vertex_values[3,1]
360
361
362       
363##         #Assert that quantities are conserved
364##         from Numeric import sum
365##         for k in range(quantity.centroid_values.shape[0]):
366##             assert allclose (quantity.centroid_values[k],
367##                              sum(quantity.vertex_values[k,:])/3)
368       
369
370
371##     def test_distribute_first_order(self):
372##         quantity = Conserved_quantity(self.mesh4)       
373
374##         #Test centroids
375##         quantity.set_values([1.,2.,3.,4.], location = 'centroids')
376##         assert allclose(quantity.centroid_values, [1, 2, 3, 4]) #Centroid
377
378
379##         #Extrapolate
380##         quantity.extrapolate_first_order()
381
382##         #Interpolate
383##         quantity.interpolate_from_vertices_to_edges()       
384
385##         assert allclose(quantity.vertex_values,
386##                         [[1,1,1], [2,2,2], [3,3,3], [4, 4, 4]])
387##         assert allclose(quantity.edge_values, [[1,1,1], [2,2,2],
388##                                                [3,3,3], [4, 4, 4]])
389
390
391
392##     def test_update_explicit(self):
393##         quantity = Conserved_quantity(self.mesh4)
394
395##         #Test centroids
396##         quantity.set_values([1.,2.,3.,4.], location = 'centroids')
397##         assert allclose(quantity.centroid_values, [1, 2, 3, 4]) #Centroid
398
399##         #Set explicit_update
400##         quantity.explicit_update = array( [1.,1.,1.,1.] )
401
402##         #Update with given timestep
403##         quantity.update(0.1)
404
405##         x = array([1, 2, 3, 4]) + array( [.1,.1,.1,.1] )
406##         assert allclose( quantity.centroid_values, x)
407
408##     def test_update_semi_implicit(self):
409##         quantity = Conserved_quantity(self.mesh4)
410
411##         #Test centroids
412##         quantity.set_values([1.,2.,3.,4.], location = 'centroids')
413##         assert allclose(quantity.centroid_values, [1, 2, 3, 4]) #Centroid
414
415##         #Set semi implicit update
416##         quantity.semi_implicit_update = array( [1.,1.,1.,1.] )
417
418##         #Update with given timestep
419##         quantity.update(0.1)
420
421##         x = array([1, 2, 3, 4])/array( [.9,.9,.9,.9] )
422##         assert allclose( quantity.centroid_values, x)
423
424##     def test_both_updates(self):
425##         quantity = Conserved_quantity(self.mesh4)
426
427##         #Test centroids
428##         quantity.set_values([1.,2.,3.,4.], location = 'centroids')
429##         assert allclose(quantity.centroid_values, [1, 2, 3, 4]) #Centroid
430
431##         #Set explicit_update
432##         quantity.explicit_update = array( [4.,3.,2.,1.] )
433       
434##         #Set semi implicit update
435##         quantity.semi_implicit_update = array( [1.,1.,1.,1.] )
436
437##         #Update with given timestep
438##         quantity.update(0.1)
439
440##         x = array([1, 2, 3, 4]) + array( [.4,.3,.2,.1] )
441##         x /= array( [.9,.9,.9,.9] )
442##         assert allclose( quantity.centroid_values, x)       
443
444       
445       
446#-------------------------------------------------------------
447if __name__ == "__main__":
448    suite = unittest.makeSuite(TestCase,'test')
449    runner = unittest.TextTestRunner()
450    runner.run(suite)
451
452
Note: See TracBrowser for help on using the repository browser.