source: anuga_work/development/anuga_1d/test_shallow_water.py @ 5565

Last change on this file since 5565 was 5563, checked in by steve, 16 years ago

Working version of comp_flux_ext

File size: 9.4 KB
Line 
1#!/usr/bin/env python
2
3import unittest
4from math import sqrt, pi
5
6
7from shallow_water_domain import *
8from Numeric import allclose, array, ones, Float
9
10
11class Test_Shallow_Water(unittest.TestCase):
12    def setUp(self):
13        self.points = [0.0, 1.0, 2.0, 3.0]
14        self.vertex_values = [[1.0,2.0],[4.0,5.0],[-1.0,2.0]]
15       
16    def tearDown(self):
17        pass
18        #print "  Tearing down"
19
20
21    def test_creation(self):
22        domain = Domain(self.points)
23        assert allclose(domain.centroids, [0.5, 1.5, 2.5])
24
25    def test_compute_fluxes(self):
26        """
27        Compare shallow_water_domain flux calculation against a previous
28        Python implementation (defined in this file)
29        """
30        domain = Domain(self.points)
31        domain.set_quantity('stage',2.0)
32        domain.set_boundary({'exterior' : Reflective_boundary(domain)})
33
34        stage_ud, xmom_ud = compute_fluxes_python(domain)
35       
36        domain.compute_fluxes()
37
38        #print doamin.quantities['xmomentum'].explicit_update
39        #print compute_fluxes_python(domain)
40
41        assert allclose( domain.quantities['stage'].explicit_update, stage_ud )
42        assert allclose( domain.quantities['xmomentum'].explicit_update, xmom_ud )
43
44
45    def test_local_flux_function(self):
46        normal = 1.0
47        ql = array([1.0, 2.0],Float)
48        qr = array([1.0, 2.0],Float)
49        zl = 0.0
50        zr = 0.0
51
52        edgeflux, maxspeed = flux_function(normal, ql,qr,zl,zr)
53
54        assert allclose(array([2.0, 8.9],Float), edgeflux)
55        assert allclose(5.1304951685, maxspeed)
56
57        normal = -1.0
58        ql = array([1.0, 2.0],Float)
59        qr = array([1.0, 2.0],Float)
60        zl = 0.0
61        zr = 0.0
62
63        edgeflux, maxspeed = flux_function(normal, ql,qr,zl,zr)
64
65        assert allclose(array([-2.0, -8.9],Float), edgeflux)
66        assert allclose(5.1304951685, maxspeed)
67
68    def test_domain_flux_function(self):
69        normal = 1.0
70        ql = array([1.0, 2.0],Float)
71        qr = array([1.0, 2.0],Float)
72        zl = 0.0
73        zr = 0.0
74
75        edgeflux, maxspeed = flux_function(normal, ql,qr,zl,zr)
76
77        #print edgeflux
78
79        from shallow_water_domain import flux_function as domain_flux_function
80
81        domainedgeflux, domainmaxspeed = domain_flux_function(normal, ql,qr,zl,zr)
82
83        #print domainedgeflux
84
85        assert allclose(domainedgeflux, edgeflux)
86        assert allclose(domainmaxspeed, maxspeed)
87
88
89
90#==============================================================================
91
92def compute_fluxes_python(domain):
93    """Compute all fluxes and the timestep suitable for all volumes
94    in domain.
95
96    Compute total flux for each conserved quantity using "flux_function"
97
98    Fluxes across each edge are scaled by edgelengths and summed up
99    Resulting flux is then scaled by area and stored in
100    explicit_update for each of the three conserved quantities
101    stage, xmomentum and ymomentum
102
103    The maximal allowable speed computed by the flux_function for each volume
104    is converted to a timestep that must not be exceeded. The minimum of
105    those is computed as the next overall timestep.
106
107    Post conditions:
108      domain.explicit_update is reset to computed flux values
109      domain.timestep is set to the largest step satisfying all volumes.
110    """
111
112    import sys
113    from Numeric import zeros, Float
114
115    N = domain.number_of_elements
116
117    tmp0 = zeros((N,),Float)
118    tmp1 = zeros((N,),Float)
119
120    #Shortcuts
121    Stage = domain.quantities['stage']
122    Xmom = domain.quantities['xmomentum']
123#    Ymom = domain.quantities['ymomentum']
124    Bed = domain.quantities['elevation']
125
126    #Arrays
127    #stage = Stage.edge_values
128    #xmom =  Xmom.edge_values
129 #   ymom =  Ymom.edge_values
130    #bed =   Bed.edge_values
131   
132    stage = Stage.vertex_values
133    xmom =  Xmom.vertex_values
134    bed =   Bed.vertex_values
135   
136    #print 'stage edge values', stage
137    #print 'xmom edge values', xmom
138    #print 'bed values', bed
139
140    stage_bdry = Stage.boundary_values
141    xmom_bdry =  Xmom.boundary_values
142    #print 'stage_bdry',stage_bdry
143    #print 'xmom_bdry', xmom_bdry
144#    ymom_bdry =  Ymom.boundary_values
145
146#    flux = zeros(3, Float) #Work array for summing up fluxes
147    flux = zeros(2, Float) #Work array for summing up fluxes
148    ql = zeros(2, Float)
149    qr = zeros(2, Float)
150
151    #Loop
152    timestep = float(sys.maxint)
153    enter = True
154    for k in range(N):
155
156        flux[:] = 0.  #Reset work array
157        #for i in range(3):
158        for i in range(2):
159            #Quantities inside volume facing neighbour i
160            #ql[0] = stage[k, i]
161            #ql[1] = xmom[k, i]
162            ql = [stage[k, i], xmom[k, i]]
163            zl = bed[k, i]
164
165            #Quantities at neighbour on nearest face
166            n = domain.neighbours[k,i]
167            if n < 0:
168                m = -n-1 #Convert negative flag to index
169                qr[0] = stage_bdry[m]
170                qr[1] = xmom_bdry[m]
171                zr = zl #Extend bed elevation to boundary
172            else:
173                #m = domain.neighbour_edges[k,i]
174                m = domain.neighbour_vertices[k,i]
175                #print i, ' ' , m
176                #qr = [stage[n, m], xmom[n, m], ymom[n, m]]
177                qr[0] = stage[n, m]
178                qr[1] = xmom[n, m]
179                zr = bed[n, m]
180
181
182            #Outward pointing normal vector
183            normal = domain.normals[k, i]
184       
185            #Flux computation using provided function
186            #edgeflux, max_speed = flux_function(normal, ql, qr, zl, zr)
187            #print 'ql',ql
188            #print 'qr',qr
189           
190
191            edgeflux, max_speed = flux_function(normal, ql, qr, zl, zr)
192
193            #print 'edgeflux', edgeflux
194
195            # THIS IS THE LINE TO DEAL WITH LEFT AND RIGHT FLUXES
196            # flux = edgefluxleft - edgefluxright
197            flux -= edgeflux #* domain.edgelengths[k,i]
198            #Update optimal_timestep
199            try:
200                #timestep = min(timestep, 0.5*domain.radii[k]/max_speed)
201                timestep = min(timestep, domain.cfl*0.5*domain.areas[k]/max_speed)
202            except ZeroDivisionError:
203                pass
204
205        #Normalise by area and store for when all conserved
206        #quantities get updated
207        flux /= domain.areas[k]
208
209        #Stage.explicit_update[k] = flux[0]
210        tmp0[k] = flux[0]
211        tmp1[k] = flux[1]
212
213
214    return tmp0, tmp1
215
216
217def flux_function(normal, ql, qr, zl, zr):
218    """Compute fluxes between volumes for the shallow water wave equation
219    cast in terms of w = h+z using the 'central scheme' as described in
220
221    Kurganov, Noelle, Petrova. 'Semidiscrete Central-Upwind Schemes For
222    Hyperbolic Conservation Laws and Hamilton-Jacobi Equations'.
223    Siam J. Sci. Comput. Vol. 23, No. 3, pp. 707-740.
224
225    The implemented formula is given in equation (3.15) on page 714
226
227    Conserved quantities w, uh, are stored as elements 0 and 1
228    in the numerical vectors ql an qr.
229
230    Bed elevations zl and zr.
231    """
232
233    from config import g, epsilon
234    from math import sqrt
235    from Numeric import array
236
237    #print 'ql',ql
238
239    #Align momentums with x-axis
240    #q_left  = rotate(ql, normal, direction = 1)
241    #q_right = rotate(qr, normal, direction = 1)
242    q_left = ql
243    q_left[1] = q_left[1]*normal
244    q_right = qr
245    q_right[1] = q_right[1]*normal
246
247    #z = (zl+zr)/2 #Take average of field values
248    z = 0.5*(zl+zr) #Take average of field values
249
250    w_left  = q_left[0]   #w=h+z
251    h_left  = w_left-z
252    uh_left = q_left[1]
253
254    if h_left < epsilon:
255        u_left = 0.0  #Could have been negative
256        h_left = 0.0
257    else:
258        u_left  = uh_left/h_left
259
260
261    w_right  = q_right[0]  #w=h+z
262    h_right  = w_right-z
263    uh_right = q_right[1]
264
265
266    if h_right < epsilon:
267        u_right = 0.0  #Could have been negative
268        h_right = 0.0
269    else:
270        u_right  = uh_right/h_right
271
272    #vh_left  = q_left[2]
273    #vh_right = q_right[2]
274
275    #print h_right
276    #print u_right
277    #print h_left
278    #print u_right
279
280    soundspeed_left  = sqrt(g*h_left)
281    soundspeed_right = sqrt(g*h_right)
282
283    #Maximal wave speed
284    s_max = max(u_left+soundspeed_left, u_right+soundspeed_right, 0)
285   
286    #Minimal wave speed
287    s_min = min(u_left-soundspeed_left, u_right-soundspeed_right, 0)
288   
289    #Flux computation
290
291    #flux_left  = array([u_left*h_left,
292    #                    u_left*uh_left + 0.5*g*h_left**2])
293    #flux_right = array([u_right*h_right,
294    #                    u_right*uh_right + 0.5*g*h_right**2])
295    flux_left  = array([u_left*h_left,
296                        u_left*uh_left + 0.5*g*h_left*h_left])
297    flux_right = array([u_right*h_right,
298                        u_right*uh_right + 0.5*g*h_right*h_right])
299
300    denom = s_max-s_min
301    if denom == 0.0:
302        edgeflux = array([0.0, 0.0])
303        max_speed = 0.0
304    else:
305        edgeflux = (s_max*flux_left - s_min*flux_right)/denom
306        edgeflux += s_max*s_min*(q_right-q_left)/denom
307       
308        edgeflux[1] = edgeflux[1]*normal
309
310        max_speed = max(abs(s_max), abs(s_min))
311
312    return edgeflux, max_speed
313
314
315#-------------------------------------------------------------
316if __name__ == "__main__":
317    suite = unittest.makeSuite(Test_Shallow_Water, 'test')   
318    #suite = unittest.makeSuite(Test_Quantity, 'test_set_values_from_file_using_polygon')
319
320    #suite = unittest.makeSuite(Test_Quantity, 'test_set_vertex_values_using_general_interface_with_subset')
321    #print "restricted test"
322    #suite = unittest.makeSuite(Test_Quantity,'verbose_test_set_values_from_UTM_pts')
323    runner = unittest.TextTestRunner()
324    runner.run(suite)
Note: See TracBrowser for help on using the repository browser.