source: trunk/anuga_work/shallow_water_balanced_steve/test_swb_domain_ext.py @ 8288

Last change on this file since 8288 was 8175, checked in by steve, 13 years ago

Working on balanced scheme
File size: 17.1 KB
Line 
1#!/usr/bin/env python
2
3import unittest, os
4import os.path
5from math import pi, sqrt
6
7from swb_domain import *
8
9from anuga.config import g
10
11import numpy as num
12
13# Get gateway to C implementation of flux function for direct testing
14from swb_domain_ext import flux_function_c
15from swb_domain_ext import gravity_c
16from swb_domain_ext import compute_fluxes_c
17
18
19class Test_swb_domain_ext(unittest.TestCase):
20    def setUp(self):
21        pass
22
23    def tearDown(self):
24        pass
25
26    def test_compute_fluxes0(self):
27        # Do a full triangle and check that fluxes cancel out for
28        # the constant stage case
29
30        a = [0.0, 0.0]
31        b = [0.0, 2.0]
32        c = [2.0, 0.0]
33        d = [0.0, 4.0]
34        e = [2.0, 2.0]
35        f = [4.0, 0.0]
36
37        points = [a, b, c, d, e, f]
38        #             bac,     bce,     ecf,     dbe
39        vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
40
41        domain = Domain(points, vertices)
42        domain.set_quantity('stage', [[2,2,2], [2,2,2], [2,2,2], [2,2,2]])
43        domain.check_integrity()
44        domain.distribute_to_vertices_and_edges()
45
46        assert num.allclose(domain.neighbours,
47                            [[-1,1,-2], [2,3,0], [-3,-4,1],[1,-5,-6]])
48        assert num.allclose(domain.neighbour_edges,
49                            [[-1,2,-1], [2,0,1], [-1,-1,0],[1,-1,-1]])
50
51        edgeflux = num.zeros(3, num.float)
52        edgeflux0 = num.zeros(3, num.float)
53        edgeflux1 = num.zeros(3, num.float)
54        edgeflux2 = num.zeros(3, num.float)
55
56        # Flux across right edge of volume 1
57        normal = domain.get_normal(1, 0)
58        ql = domain.get_evolved_quantities(vol_id=1, edge=0)
59        qr = domain.get_evolved_quantities(vol_id=2, edge=2)
60
61        #print ql
62        #print qr
63        local_max_speed = flux_function_c(normal, ql, qr, edgeflux0, g)
64
65        # Check that flux seen from other triangles is inverse
66        (ql, qr) = (qr, ql)
67        normal = domain.get_normal(2, 2)
68        local_max_speed = flux_function_c(normal, ql, qr, edgeflux, g)
69
70        assert num.allclose(edgeflux0 + edgeflux, 0.)
71
72        # Flux across upper edge of volume 1
73        normal = domain.get_normal(1, 1)
74        ql = domain.get_evolved_quantities(vol_id=1, edge=1)
75        qr = domain.get_evolved_quantities(vol_id=3, edge=0)
76        #print ql
77        #print qr
78        local_max_speed = flux_function_c(normal, ql, qr, edgeflux1, g)
79
80        # Check that flux seen from other triangles is inverse
81        (ql, qr) = (qr, ql)
82        normal = domain.get_normal(3, 0)
83        local_max_speed = flux_function_c(normal, ql, qr, edgeflux, g)
84
85        assert num.allclose(edgeflux1 + edgeflux, 0.)
86
87        # Flux across lower left hypotenuse of volume 1
88        normal = domain.get_normal(1, 2)
89        ql = domain.get_evolved_quantities(vol_id=1, edge=2)
90        qr = domain.get_evolved_quantities(vol_id=0, edge=1)
91        #print ql
92        #print qr
93        local_max_speed = flux_function_c(normal, ql, qr, edgeflux2, g)
94
95        # Check that flux seen from other triangles is inverse
96        (ql, qr) = (qr, ql)
97        normal = domain.get_normal(0, 1)
98        max_speed = flux_function_c(normal, ql, qr, edgeflux, g)
99        assert num.allclose(edgeflux2 + edgeflux, 0.)
100
101        # Scale by edgelengths, add up anc check that total flux is zero
102        e0 = domain.edgelengths[1, 0]
103        e1 = domain.edgelengths[1, 1]
104        e2 = domain.edgelengths[1, 2]
105
106        assert num.allclose(e0*edgeflux0 + e1*edgeflux1 + e2*edgeflux2, 0.)
107
108        # Now check that compute_flux yields zeros as well
109        w  = domain.quantities['stage']
110        uh = domain.quantities['xmomentum']
111        vh = domain.quantities['ymomentum']
112        h  = domain.quantities['height']
113        z  = domain.quantities['elevation']
114        u  = domain.quantities['xvelocity']
115        v  = domain.quantities['yvelocity']
116        timestep = compute_fluxes_c(2.0,domain,w,uh,vh,h,z,u,v)
117
118        assert num.allclose(timestep,  0.1064794275)
119
120
121        for name in ['stage', 'xmomentum', 'ymomentum']:
122            assert num.allclose(domain.quantities[name].explicit_update[1], 0.0)
123
124    def test_compute_fluxes1(self):
125        #Use values from previous version
126        a = [0.0, 0.0]
127        b = [0.0, 2.0]
128        c = [2.0, 0.0]
129        d = [0.0, 4.0]
130        e = [2.0, 2.0]
131        f = [4.0, 0.0]
132
133        points = [a, b, c, d, e, f]
134        #             bac,     bce,     ecf,     dbe
135        vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
136
137        domain = Domain(points, vertices)
138        val0 = 2. + 2.0/3
139        val1 = 4. + 4.0/3
140        val2 = 8. + 2.0/3
141        val3 = 2. + 8.0/3
142
143        domain.set_quantity('stage', [[val0, val0, val0], [val1, val1, val1],
144                                      [val2, val2, val2], [val3, val3, val3]])
145        domain.check_integrity()
146        domain.set_default_order(1)
147        domain.distribute_to_vertices_and_edges()
148
149
150        w  = domain.quantities['stage']
151        uh = domain.quantities['xmomentum']
152        vh = domain.quantities['ymomentum']
153        h  = domain.quantities['height']
154        z  = domain.quantities['elevation']
155        u  = domain.quantities['xvelocity']
156        v  = domain.quantities['yvelocity']
157
158
159        edgeflux = num.zeros(3, num.float)
160        edgeflux0 = num.zeros(3, num.float)
161        edgeflux1 = num.zeros(3, num.float)
162        edgeflux2 = num.zeros(3, num.float)
163
164
165        # Flux across right edge of volume 1
166        normal = domain.get_normal(1, 0)    # Get normal 0 of triangle 1
167        assert num.allclose(normal, [1, 0])
168
169 
170        ql = domain.get_evolved_quantities(vol_id=1, edge=0)
171        #print ql, val1
172        assert num.allclose(ql, [val1, 0, 0, val1, 0, 0, 0])
173
174        qr = domain.get_evolved_quantities(vol_id=2, edge=2)
175        assert num.allclose(qr, [val2, 0, 0, val2, 0, 0, 0])
176
177        local_max_speed = flux_function_c(normal, ql, qr, edgeflux0, g)
178
179        #print edgeflux0
180        #print local_max_speed
181       
182        # Flux across edge in the east direction (as per normal vector)
183        assert num.allclose(edgeflux0, [-15.3598804, 253.71111111, 0.])
184        assert num.allclose(local_max_speed, 9.21592824046)
185
186        #Flux across edge in the west direction (opposite sign for xmomentum)
187        normal_opposite = domain.get_normal(2, 2)   # Get normal 2 of triangle 2
188        assert num.allclose(normal_opposite, [-1, 0])
189
190        max_speed = flux_function_c(normal_opposite, ql, qr, edgeflux, g)
191        assert num.allclose(edgeflux, [-15.3598804, -253.71111111, 0.])
192
193        #Flux across upper edge of volume 1
194        normal = domain.get_normal(1, 1)
195        ql = domain.get_evolved_quantities(vol_id=1, edge=1)
196        qr = domain.get_evolved_quantities(vol_id=3, edge=0)
197        max_speed = flux_function_c(normal, ql, qr, edgeflux1, g)
198
199        assert num.allclose(edgeflux1, [2.4098563, 0., 123.04444444])
200        assert num.allclose(max_speed, 7.22956891292)
201
202        #Flux across lower left hypotenuse of volume 1
203        normal = domain.get_normal(1, 2)
204        ql = domain.get_evolved_quantities(vol_id=1, edge=2)
205        qr = domain.get_evolved_quantities(vol_id=0, edge=1)
206        max_speed = flux_function_c(normal, ql, qr, edgeflux2, g)
207
208        assert num.allclose(edgeflux2, [9.63942522, -61.59685738, -61.59685738])
209        assert num.allclose(max_speed, 7.22956891292)
210
211        #Scale, add up and check that compute_fluxes is correct for vol 1
212        e0 = domain.edgelengths[1, 0]
213        e1 = domain.edgelengths[1, 1]
214        e2 = domain.edgelengths[1, 2]
215
216        total_flux = -(e0*edgeflux0 +
217                       e1*edgeflux1 +
218                       e2*edgeflux2) / domain.areas[1]
219
220        assert num.allclose(total_flux, [-0.68218178, -166.6, -35.93333333])
221
222
223        timestep = compute_fluxes_c(2.0,domain,w,uh,vh,h,z,u,v)
224        assert num.allclose(timestep, 0.0511510624314)
225
226        for i, name in enumerate(['stage', 'xmomentum', 'ymomentum']):
227            assert num.allclose(total_flux[i],
228                                domain.quantities[name].explicit_update[1])
229
230
231        assert num.allclose(domain.quantities['stage'].explicit_update,
232                            [0., -0.68218178, -111.77316251, -35.68522449])
233        assert num.allclose(domain.quantities['xmomentum'].explicit_update,
234                            [-69.68888889, -166.6, 69.68888889, 0])
235        assert num.allclose(domain.quantities['ymomentum'].explicit_update,
236                            [-69.68888889, -35.93333333, 0., 69.68888889])
237
238    def test_compute_fluxes2(self):
239        #Random values, incl momentum
240        a = [0.0, 0.0]
241        b = [0.0, 2.0]
242        c = [2.0, 0.0]
243        d = [0.0, 4.0]
244        e = [2.0, 2.0]
245        f = [4.0, 0.0]
246
247        points = [a, b, c, d, e, f]
248        #             bac,     bce,     ecf,     dbe
249        vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
250
251        domain = Domain(points, vertices)
252        domain.set_default_order(1)
253
254        w  = domain.quantities['stage']
255        uh = domain.quantities['xmomentum']
256        vh = domain.quantities['ymomentum']
257        h  = domain.quantities['height']
258        z  = domain.quantities['elevation']
259        u  = domain.quantities['xvelocity']
260        v  = domain.quantities['yvelocity']
261
262        zl = 10    # Assume flat bed
263        val0 = zl + 2. + 2.0/3
264        val1 = zl + 4. + 4.0/3
265        val2 = zl + 8. + 2.0/3
266        val3 = zl + 2. + 8.0/3
267
268
269        edgeflux = num.zeros(3, num.float)
270        edgeflux0 = num.zeros(3, num.float)
271        edgeflux1 = num.zeros(3, num.float)
272        edgeflux2 = num.zeros(3, num.float)
273
274        domain.set_quantity('elevation', zl*num.ones((4, 3), num.int)) #array default#
275
276        domain.set_quantity('stage', [[val0, val0-1, val0-2],
277                                      [val1, val1+1, val1],
278                                      [val2, val2-2, val2],
279                                      [val3-0.5, val3, val3]])
280
281        domain.set_quantity('xmomentum',
282                            [[1,2,3], [3,4,5], [1,-1,0], [0,-2,2]])
283
284        domain.set_quantity('ymomentum',
285                            [[1,-1,0], [0,-3,2], [0,1,0], [-1,2,2]])
286
287        domain.check_integrity()
288        domain.distribute_to_vertices_and_edges()
289
290       
291        # Flux across right edge of volume 1
292        normal = domain.get_normal(1, 0)
293        ql = domain.get_evolved_quantities(vol_id=1, edge=0)
294        qr = domain.get_evolved_quantities(vol_id=2, edge=2)
295        max_speed = flux_function_c(normal, ql, qr, edgeflux0, g)
296
297        # Flux across upper edge of volume 1
298        normal = domain.get_normal(1, 1)
299        ql = domain.get_evolved_quantities(vol_id=1, edge=1)
300        qr = domain.get_evolved_quantities(vol_id=3, edge=0)
301        max_speed = flux_function_c(normal, ql, qr, edgeflux1, g)
302
303        # Flux across lower left hypotenuse of volume 1
304        normal = domain.get_normal(1, 2)
305        ql = domain.get_evolved_quantities(vol_id=1, edge=2)
306        qr = domain.get_evolved_quantities(vol_id=0, edge=1)
307        max_speed = flux_function_c(normal, ql, qr, edgeflux2, g)
308
309        # Scale, add up and check that compute_fluxes is correct for vol 1
310        e0 = domain.edgelengths[1, 0]
311        e1 = domain.edgelengths[1, 1]
312        e2 = domain.edgelengths[1, 2]
313
314        total_flux = -(e0*edgeflux0 +
315                       e1*edgeflux1 +
316                       e2*edgeflux2) / domain.areas[1]
317
318        timestep = compute_fluxes_c(2.0,domain,w,uh,vh,h,z,u,v)
319        assert num.allclose(timestep, 0.0529903533177)
320
321        for i, name in enumerate(['stage', 'xmomentum', 'ymomentum']):
322            assert num.allclose(total_flux[i],
323                                domain.quantities[name].explicit_update[1])
324
325    # FIXME (Ole): Need test like this for fluxes in very shallow water.
326    def test_compute_fluxes3(self):
327        #Random values, incl momentum
328        a = [0.0, 0.0]
329        b = [0.0, 2.0]
330        c = [2.0, 0.0]
331        d = [0.0, 4.0]
332        e = [2.0, 2.0]
333        f = [4.0, 0.0]
334
335        points = [a, b, c, d, e, f]
336        #             bac,     bce,     ecf,     dbe
337        vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
338
339        domain = Domain(points, vertices)
340        domain.set_default_order(1)
341
342
343        w  = domain.quantities['stage']
344        uh = domain.quantities['xmomentum']
345        vh = domain.quantities['ymomentum']
346        h  = domain.quantities['height']
347        z  = domain.quantities['elevation']
348        u  = domain.quantities['xvelocity']
349        v  = domain.quantities['yvelocity']
350       
351        val0 = 2.+2.0/3
352        val1 = 4.+4.0/3
353        val2 = 8.+2.0/3
354        val3 = 2.+8.0/3
355
356        zl = -3.75    # Assume constant bed (must be less than stage)
357        domain.set_quantity('elevation', zl*num.ones((4, 3), num.int)) #array default#
358
359        edgeflux = num.zeros(3, num.float)
360        edgeflux0 = num.zeros(3, num.float)
361        edgeflux1 = num.zeros(3, num.float)
362        edgeflux2 = num.zeros(3, num.float)
363
364
365        domain.set_quantity('stage', [[val0, val0-1, val0-2],
366                                      [val1, val1+1, val1],
367                                      [val2, val2-2, val2],
368                                      [val3-0.5, val3, val3]])
369
370        domain.set_quantity('xmomentum',
371                            [[1,2,3], [3,4,5], [1,-1,0], [0,-2,2]])
372
373        domain.set_quantity('ymomentum',
374                            [[1,-1,0], [0,-3,2], [0,1,0], [-1,2,2]])
375
376        domain.check_integrity()
377        domain.distribute_to_vertices_and_edges()       
378
379        # Flux across right edge of volume 1
380        normal = domain.get_normal(1, 0)
381        ql = domain.get_evolved_quantities(vol_id=1, edge=0)
382        qr = domain.get_evolved_quantities(vol_id=2, edge=2)
383        max_speed = flux_function_c(normal, ql, qr, edgeflux0, g)
384
385        #print max_speed, edgeflux0
386        assert num.allclose(edgeflux0, [-10.51926294, 577.81462335, -3.64772873], rtol=1.0e-2)
387
388        # Flux across upper edge of volume 1
389        normal = domain.get_normal(1, 1)
390        ql = domain.get_evolved_quantities(vol_id=1, edge=1)
391        qr = domain.get_evolved_quantities(vol_id=3, edge=0)
392        max_speed = flux_function_c(normal, ql, qr, edgeflux1, g)
393
394        #print max_speed, edgeflux1
395        assert num.allclose(edgeflux1, [ 5.93945967, 19.14204647,  377.48851296], rtol=1.0e-2)
396       
397        # Flux across lower left hypotenuse of volume 1
398        normal = domain.get_normal(1, 2)
399        ql = domain.get_evolved_quantities(vol_id=1, edge=2)
400        qr = domain.get_evolved_quantities(vol_id=0, edge=1)
401        max_speed = flux_function_c(normal, ql, qr, edgeflux2, g)
402
403        #print max_speed, edgeflux2
404        assert num.allclose(edgeflux2, [17.21047971, -192.90578267, -203.05771337], rtol=1.0e-2)
405
406       
407        # Scale, add up and check that compute_fluxes is correct for vol 1
408        e0 = domain.edgelengths[1, 0]
409        e1 = domain.edgelengths[1, 1]
410        e2 = domain.edgelengths[1, 2]
411
412        total_flux = -(e0*edgeflux0 +
413                       e1*edgeflux1 +
414                       e2*edgeflux2) / domain.areas[1]
415
416        timestep = compute_fluxes_c(100.0,domain,w,uh,vh,h,z,u,v)
417
418        assert num.allclose(timestep, 0.0438142267244, rtol=1.0e-2)
419
420        for i, name in enumerate(['stage', 'xmomentum', 'ymomentum']):
421            assert num.allclose(total_flux[i],
422                                domain.quantities[name].explicit_update[1])
423
424
425
426    def test_compute_fluxes4(self):
427        #Random values, incl momentum
428        a = [0.0, 0.0]
429        b = [0.0, 2.0]
430        c = [2.0, 0.0]
431        d = [0.0, 4.0]
432        e = [2.0, 2.0]
433        f = [4.0, 0.0]
434
435        points = [a, b, c, d, e, f]
436        #             bac,     bce,     ecf,     dbe
437        vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
438
439        domain = Domain(points, vertices)
440        domain.set_default_order(1)
441
442
443        w  = domain.quantities['stage']
444        uh = domain.quantities['xmomentum']
445        vh = domain.quantities['ymomentum']
446        h  = domain.quantities['height']
447        z  = domain.quantities['elevation']
448        u  = domain.quantities['xvelocity']
449        v  = domain.quantities['yvelocity']
450       
451        val0 = 2.+2.0/3
452        val1 = 4.+4.0/3
453        val2 = 8.+2.0/3
454        val3 = 2.+8.0/3
455
456        zl = -3.75    # Assume constant bed (must be less than stage)
457        domain.set_quantity('elevation', [0.0, -2.0, 0.0 ,0.0], location='centroids') 
458
459        edgeflux = num.zeros(3, num.float)
460        edgeflux0 = num.zeros(3, num.float)
461        edgeflux1 = num.zeros(3, num.float)
462        edgeflux2 = num.zeros(3, num.float)
463
464
465        domain.set_quantity('stage', -1.0)
466
467        domain.set_quantity('xmomentum',0.0)
468
469        domain.set_quantity('ymomentum',0.0)
470
471        domain.check_integrity()
472        domain.distribute_to_vertices_and_edges()       
473
474
475        #print w.edge_values
476        #print u.edge_values
477        #print v.edge_values
478        #print h.edge_values
479        #print z.edge_values
480
481
482        domain.compute_fluxes()
483
484        #for i, name in enumerate(['stage', 'xmomentum', 'ymomentum']):
485            #print domain.quantities[name].explicit_update[1]
486            #assert num.allclose(total_flux[i],
487            #                    domain.quantities[name].explicit_update[1])
488
489
490#################################################################################
491
492if __name__ == "__main__":
493    suite = unittest.makeSuite(Test_swb_domain_ext, 'test')
494    runner = unittest.TextTestRunner(verbosity=1)
495    runner.run(suite)
Note: See TracBrowser for help on using the repository browser.