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

Last change on this file since 5727 was 5727, checked in by steve, 15 years ago

Setting up some more unit tests for shallow_water_domain

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