source: anuga_core/source/anuga/shallow_water_balanced/test_swb_compute_flux.py @ 7559

Last change on this file since 7559 was 7559, checked in by steve, 14 years ago

Broke up test_swb_domain so that we can concentrate on new functionality
File size: 16.4 KB
Line 
1#!/usr/bin/env python
2
3import unittest, os
4import os.path
5from math import pi, sqrt
6import tempfile
7
8from anuga.config import g, epsilon
9from anuga.config import netcdf_mode_r, netcdf_mode_w, netcdf_mode_a
10from anuga.utilities.numerical_tools import mean
11from anuga.utilities.polygon import is_inside_polygon
12from anuga.coordinate_transforms.geo_reference import Geo_reference
13from anuga.abstract_2d_finite_volumes.quantity import Quantity
14from anuga.geospatial_data.geospatial_data import Geospatial_data
15from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular_cross
16
17from anuga.utilities.system_tools import get_pathname_from_package
18from swb_domain import *
19
20import numpy as num
21
22# Get gateway to C implementation of flux function for direct testing
23from shallow_water_ext import flux_function_central as flux_function
24
25
26
27
28class Test_swb_clean(unittest.TestCase):
29    def setUp(self):
30        pass
31
32    def tearDown(self):
33        pass
34
35    def test_compute_fluxes0(self):
36        # Do a full triangle and check that fluxes cancel out for
37        # the constant stage case
38
39        a = [0.0, 0.0]
40        b = [0.0, 2.0]
41        c = [2.0, 0.0]
42        d = [0.0, 4.0]
43        e = [2.0, 2.0]
44        f = [4.0, 0.0]
45
46        points = [a, b, c, d, e, f]
47        #             bac,     bce,     ecf,     dbe
48        vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
49
50        domain = Domain(points, vertices)
51        domain.set_quantity('stage', [[2,2,2], [2,2,2], [2,2,2], [2,2,2]])
52        domain.check_integrity()
53
54        assert num.allclose(domain.neighbours,
55                            [[-1,1,-1], [2,3,0], [-1,-1,1],[1,-1,-1]])
56        assert num.allclose(domain.neighbour_edges,
57                            [[-1,2,-1], [2,0,1], [-1,-1,0],[1,-1,-1]])
58
59        zl = zr = 0.     # Assume flat bed
60
61        edgeflux = num.zeros(3, num.float)
62        edgeflux0 = num.zeros(3, num.float)
63        edgeflux1 = num.zeros(3, num.float)
64        edgeflux2 = num.zeros(3, num.float)
65        H0 = 0.0
66
67        # Flux across right edge of volume 1
68        normal = domain.get_normal(1, 0)
69        ql = domain.get_conserved_quantities(vol_id=1, edge=0)
70        qr = domain.get_conserved_quantities(vol_id=2, edge=2)
71        max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux0,
72                                  epsilon, g, H0)
73
74        # Check that flux seen from other triangles is inverse
75        (ql, qr) = (qr, ql)
76        normal = domain.get_normal(2, 2)
77        max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux,
78                                  epsilon, g, H0)
79
80        assert num.allclose(edgeflux0 + edgeflux, 0.)
81
82        # Flux across upper edge of volume 1
83        normal = domain.get_normal(1, 1)
84        ql = domain.get_conserved_quantities(vol_id=1, edge=1)
85        qr = domain.get_conserved_quantities(vol_id=3, edge=0)
86        max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux1,
87                                  epsilon, g, H0)
88
89        # Check that flux seen from other triangles is inverse
90        (ql, qr) = (qr, ql)
91        normal = domain.get_normal(3, 0)
92        max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux,
93                                  epsilon, g, H0)
94
95        assert num.allclose(edgeflux1 + edgeflux, 0.)
96
97        # Flux across lower left hypotenuse of volume 1
98        normal = domain.get_normal(1, 2)
99        ql = domain.get_conserved_quantities(vol_id=1, edge=2)
100        qr = domain.get_conserved_quantities(vol_id=0, edge=1)
101        max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux2,
102                                  epsilon, g, H0)
103
104        # Check that flux seen from other triangles is inverse
105        (ql, qr) = (qr, ql)
106        normal = domain.get_normal(0, 1)
107        max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux,
108                                  epsilon, g, H0)
109        assert num.allclose(edgeflux2 + edgeflux, 0.)
110
111        # Scale by edgelengths, add up anc check that total flux is zero
112        e0 = domain.edgelengths[1, 0]
113        e1 = domain.edgelengths[1, 1]
114        e2 = domain.edgelengths[1, 2]
115
116        assert num.allclose(e0*edgeflux0 + e1*edgeflux1 + e2*edgeflux2, 0.)
117
118        # Now check that compute_flux yields zeros as well
119        domain.compute_fluxes()
120
121        for name in ['stage', 'xmomentum', 'ymomentum']:
122            assert num.allclose(domain.quantities[name].explicit_update[1], 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
147        zl = zr = 0.    # Assume flat bed
148
149        edgeflux = num.zeros(3, num.float)
150        edgeflux0 = num.zeros(3, num.float)
151        edgeflux1 = num.zeros(3, num.float)
152        edgeflux2 = num.zeros(3, num.float)
153        H0 = 0.0
154
155        # Flux across right edge of volume 1
156        normal = domain.get_normal(1, 0)    # Get normal 0 of triangle 1
157        assert num.allclose(normal, [1, 0])
158
159        ql = domain.get_conserved_quantities(vol_id=1, edge=0)
160        assert num.allclose(ql, [val1, 0, 0])
161
162        qr = domain.get_conserved_quantities(vol_id=2, edge=2)
163        assert num.allclose(qr, [val2, 0, 0])
164
165        max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux0,
166                                  epsilon, g, H0)
167
168        # Flux across edge in the east direction (as per normal vector)
169        assert num.allclose(edgeflux0, [-15.3598804, 253.71111111, 0.])
170        assert num.allclose(max_speed, 9.21592824046)
171
172        #Flux across edge in the west direction (opposite sign for xmomentum)
173        normal_opposite = domain.get_normal(2, 2)   # Get normal 2 of triangle 2
174        assert num.allclose(normal_opposite, [-1, 0])
175
176        max_speed = flux_function(normal_opposite, ql, qr, zl, zr, edgeflux,
177                                  epsilon, g, H0)
178        assert num.allclose(edgeflux, [-15.3598804, -253.71111111, 0.])
179
180        #Flux across upper edge of volume 1
181        normal = domain.get_normal(1, 1)
182        ql = domain.get_conserved_quantities(vol_id=1, edge=1)
183        qr = domain.get_conserved_quantities(vol_id=3, edge=0)
184        max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux1,
185                                  epsilon, g, H0)
186
187        assert num.allclose(edgeflux1, [2.4098563, 0., 123.04444444])
188        assert num.allclose(max_speed, 7.22956891292)
189
190        #Flux across lower left hypotenuse of volume 1
191        normal = domain.get_normal(1, 2)
192        ql = domain.get_conserved_quantities(vol_id=1, edge=2)
193        qr = domain.get_conserved_quantities(vol_id=0, edge=1)
194        max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux2,
195                                  epsilon, g, H0)
196
197        assert num.allclose(edgeflux2, [9.63942522, -61.59685738, -61.59685738])
198        assert num.allclose(max_speed, 7.22956891292)
199
200        #Scale, add up and check that compute_fluxes is correct for vol 1
201        e0 = domain.edgelengths[1, 0]
202        e1 = domain.edgelengths[1, 1]
203        e2 = domain.edgelengths[1, 2]
204
205        total_flux = -(e0*edgeflux0 +
206                       e1*edgeflux1 +
207                       e2*edgeflux2) / domain.areas[1]
208
209        assert num.allclose(total_flux, [-0.68218178, -166.6, -35.93333333])
210
211        domain.compute_fluxes()
212
213        for i, name in enumerate(['stage', 'xmomentum', 'ymomentum']):
214            assert num.allclose(total_flux[i],
215                                domain.quantities[name].explicit_update[1])
216
217        assert num.allclose(domain.quantities['stage'].explicit_update,
218                            [0., -0.68218178, -111.77316251, -35.68522449])
219        assert num.allclose(domain.quantities['xmomentum'].explicit_update,
220                            [-69.68888889, -166.6, 69.68888889, 0])
221        assert num.allclose(domain.quantities['ymomentum'].explicit_update,
222                            [-69.68888889, -35.93333333, 0., 69.68888889])
223
224    def test_compute_fluxes2(self):
225        #Random values, incl momentum
226        a = [0.0, 0.0]
227        b = [0.0, 2.0]
228        c = [2.0, 0.0]
229        d = [0.0, 4.0]
230        e = [2.0, 2.0]
231        f = [4.0, 0.0]
232
233        points = [a, b, c, d, e, f]
234        #             bac,     bce,     ecf,     dbe
235        vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
236
237        domain = Domain(points, vertices)
238        val0 = 2. + 2.0/3
239        val1 = 4. + 4.0/3
240        val2 = 8. + 2.0/3
241        val3 = 2. + 8.0/3
242
243        zl = zr = 0    # Assume flat zero bed
244        edgeflux = num.zeros(3, num.float)
245        edgeflux0 = num.zeros(3, num.float)
246        edgeflux1 = num.zeros(3, num.float)
247        edgeflux2 = num.zeros(3, num.float)
248        H0 = 0.0
249
250        domain.set_quantity('elevation', zl*num.ones((4, 3), num.int)) #array default#
251
252        domain.set_quantity('stage', [[val0, val0-1, val0-2],
253                                      [val1, val1+1, val1],
254                                      [val2, val2-2, val2],
255                                      [val3-0.5, val3, val3]])
256
257        domain.set_quantity('xmomentum',
258                            [[1,2,3], [3,4,5], [1,-1,0], [0,-2,2]])
259
260        domain.set_quantity('ymomentum',
261                            [[1,-1,0], [0,-3,2], [0,1,0], [-1,2,2]])
262
263        domain.check_integrity()
264
265        # Flux across right edge of volume 1
266        normal = domain.get_normal(1, 0)
267        ql = domain.get_conserved_quantities(vol_id=1, edge=0)
268        qr = domain.get_conserved_quantities(vol_id=2, edge=2)
269        max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux0,
270                                  epsilon, g, H0)
271
272        # Flux across upper edge of volume 1
273        normal = domain.get_normal(1, 1)
274        ql = domain.get_conserved_quantities(vol_id=1, edge=1)
275        qr = domain.get_conserved_quantities(vol_id=3, edge=0)
276        max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux1,
277                                  epsilon, g, H0)
278
279        # Flux across lower left hypotenuse of volume 1
280        normal = domain.get_normal(1, 2)
281        ql = domain.get_conserved_quantities(vol_id=1, edge=2)
282        qr = domain.get_conserved_quantities(vol_id=0, edge=1)
283        max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux2,
284                                  epsilon, g, H0)
285
286        # Scale, add up and check that compute_fluxes is correct for vol 1
287        e0 = domain.edgelengths[1, 0]
288        e1 = domain.edgelengths[1, 1]
289        e2 = domain.edgelengths[1, 2]
290
291        total_flux = -(e0*edgeflux0 +
292                       e1*edgeflux1 +
293                       e2*edgeflux2) / domain.areas[1]
294
295        domain.compute_fluxes()
296
297        for i, name in enumerate(['stage', 'xmomentum', 'ymomentum']):
298            assert num.allclose(total_flux[i],
299                                domain.quantities[name].explicit_update[1])
300
301    # FIXME (Ole): Need test like this for fluxes in very shallow water.
302    def test_compute_fluxes3(self):
303        #Random values, incl momentum
304        a = [0.0, 0.0]
305        b = [0.0, 2.0]
306        c = [2.0, 0.0]
307        d = [0.0, 4.0]
308        e = [2.0, 2.0]
309        f = [4.0, 0.0]
310
311        points = [a, b, c, d, e, f]
312        #             bac,     bce,     ecf,     dbe
313        vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
314
315        domain = Domain(points, vertices)
316
317        val0 = 2.+2.0/3
318        val1 = 4.+4.0/3
319        val2 = 8.+2.0/3
320        val3 = 2.+8.0/3
321
322        zl = zr = -3.75    # Assume constant bed (must be less than stage)
323        domain.set_quantity('elevation', zl*num.ones((4, 3), num.int)) #array default#
324
325        edgeflux = num.zeros(3, num.float)
326        edgeflux0 = num.zeros(3, num.float)
327        edgeflux1 = num.zeros(3, num.float)
328        edgeflux2 = num.zeros(3, num.float)
329        H0 = 0.0
330
331        domain.set_quantity('stage', [[val0, val0-1, val0-2],
332                                      [val1, val1+1, val1],
333                                      [val2, val2-2, val2],
334                                      [val3-0.5, val3, val3]])
335
336        domain.set_quantity('xmomentum',
337                            [[1,2,3], [3,4,5], [1,-1,0], [0,-2,2]])
338
339        domain.set_quantity('ymomentum',
340                            [[1,-1,0], [0,-3,2], [0,1,0], [-1,2,2]])
341
342        domain.check_integrity()
343
344        # Flux across right edge of volume 1
345        normal = domain.get_normal(1, 0)
346        ql = domain.get_conserved_quantities(vol_id=1, edge=0)
347        qr = domain.get_conserved_quantities(vol_id=2, edge=2)
348        max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux0,
349                                  epsilon, g, H0)
350
351        # Flux across upper edge of volume 1
352        normal = domain.get_normal(1, 1)
353        ql = domain.get_conserved_quantities(vol_id=1, edge=1)
354        qr = domain.get_conserved_quantities(vol_id=3, edge=0)
355        max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux1,
356                                  epsilon, g, H0)
357
358        # Flux across lower left hypotenuse of volume 1
359        normal = domain.get_normal(1, 2)
360        ql = domain.get_conserved_quantities(vol_id=1, edge=2)
361        qr = domain.get_conserved_quantities(vol_id=0, edge=1)
362        max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux2,
363                                  epsilon, g, H0)
364
365        # Scale, add up and check that compute_fluxes is correct for vol 1
366        e0 = domain.edgelengths[1, 0]
367        e1 = domain.edgelengths[1, 1]
368        e2 = domain.edgelengths[1, 2]
369
370        total_flux = -(e0*edgeflux0 +
371                       e1*edgeflux1 +
372                       e2*edgeflux2) / domain.areas[1]
373
374        domain.compute_fluxes()
375
376        for i, name in enumerate(['stage', 'xmomentum', 'ymomentum']):
377            assert num.allclose(total_flux[i],
378                                domain.quantities[name].explicit_update[1])
379
380    def test_flux_optimisation(self):
381        """test_flux_optimisation
382
383        Test that fluxes are correctly computed using
384        dry and still cell exclusions
385        """
386
387        from anuga.config import g
388        import copy
389
390        a = [0.0, 0.0]
391        b = [0.0, 2.0]
392        c = [2.0, 0.0]
393        d = [0.0, 4.0]
394        e = [2.0, 2.0]
395        f = [4.0, 0.0]
396
397        points = [a, b, c, d, e, f]
398        #             bac,     bce,     ecf,     dbe
399        vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
400
401        domain = Domain(points, vertices)
402
403        #Set up for a gradient of (3,0) at mid triangle (bce)
404        def slope(x, y):
405            return 3*x
406
407        h = 0.1
408        def stage(x, y):
409            return slope(x, y) + h
410
411        domain.set_quantity('elevation', slope)
412        domain.set_quantity('stage', stage)
413
414        # Allow slope limiters to work
415        # (FIXME (Ole): Shouldn't this be automatic in ANUGA?)
416        domain.distribute_to_vertices_and_edges()
417
418        initial_stage = copy.copy(domain.quantities['stage'].vertex_values)
419
420        domain.set_boundary({'exterior': Reflective_boundary(domain)})
421
422        #  Check that update arrays are initialised to zero
423        assert num.allclose(domain.get_quantity('stage').explicit_update, 0)
424        assert num.allclose(domain.get_quantity('xmomentum').explicit_update, 0)
425        assert num.allclose(domain.get_quantity('ymomentum').explicit_update, 0)
426
427        # Get true values
428        domain.optimise_dry_cells = False
429        domain.compute_fluxes()
430        stage_ref = copy.copy(domain.get_quantity('stage').explicit_update)
431        xmom_ref = copy.copy(domain.get_quantity('xmomentum').explicit_update)
432        ymom_ref = copy.copy(domain.get_quantity('ymomentum').explicit_update)
433
434        # Try with flux optimisation
435        domain.optimise_dry_cells = True
436        domain.compute_fluxes()
437
438        assert num.allclose(stage_ref,
439                            domain.get_quantity('stage').explicit_update)
440        assert num.allclose(xmom_ref,
441                            domain.get_quantity('xmomentum').explicit_update)
442        assert num.allclose(ymom_ref,
443                            domain.get_quantity('ymomentum').explicit_update)
444
445
446#################################################################################
447
448if __name__ == "__main__":
449    suite = unittest.makeSuite(Test_swb_clean, 'test')
450    runner = unittest.TextTestRunner(verbosity=1)
451    runner.run(suite)
Note: See TracBrowser for help on using the repository browser.