source: inundation/ga/storm_surge/pyvolution/test_quantity.py @ 255

Last change on this file since 255 was 242, checked in by ole, 21 years ago

Subclassed class Quantity into Conserved_quantity

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