source: anuga_work/development/anuga_1d/test_shallow_water_vel_domain.py @ 6042

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

Found bug associated with changing max_speed_array to an integer
File size: 16.8 KB
Line 
1#!/usr/bin/env python
2
3import unittest
4from math import sqrt, pi
5
6
7from shallow_water_vel_domain import *
8
9
10from Numeric import allclose, array, ones, Float, maximum, zeros
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        self.points2 = [0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0]
18       
19    def tearDown(self):
20        pass
21        #print "  Tearing down"
22
23
24    def test_creation(self):
25        domain = Domain(self.points)
26        assert allclose(domain.centroids, [0.5, 1.5, 2.5])
27
28    def test_reflective_boundary(self):
29        """
30        Test setting reflective boundary
31        """
32       
33        domain = Domain(self.points)
34        domain.set_quantity('stage',2.0)
35        domain.set_quantity('xmomentum',6.0)
36       
37        domain.set_boundary({'exterior' : Reflective_boundary(domain)})
38
39
40        domain.distribute_to_vertices_and_edges()
41        domain.update_boundary()
42
43##         print 'In test reflective'
44##         print domain.quantities['stage'].vertex_values
45##         print domain.quantities['xmomentum'].vertex_values
46##         print domain.quantities['elevation'].vertex_values
47##         print domain.quantities['height'].vertex_values
48##         print domain.quantities['velocity'].vertex_values
49       
50##         print domain.quantities['stage'].boundary_values
51##         print domain.quantities['xmomentum'].boundary_values
52##         print domain.quantities['elevation'].boundary_values
53##         print domain.quantities['height'].boundary_values
54##         print domain.quantities['velocity'].boundary_values
55
56        assert allclose( domain.quantities['stage'    ].boundary_values, [2.0, 2.0])
57        assert allclose( domain.quantities['xmomentum'].boundary_values, [-6.0, -6.0])
58        assert allclose( domain.quantities['elevation'].boundary_values, [0.0, 0.0])
59        assert allclose( domain.quantities['height'   ].boundary_values, [2.0, 2.0])
60        assert allclose( domain.quantities['velocity' ].boundary_values, [-3.0, -3.0])
61
62       
63
64    def test_dirichlet_boundary(self):
65        """
66        Test setting dirichlet boundary
67        """
68       
69        domain = Domain(self.points)
70        domain.set_quantity('stage',2.0)
71        domain.set_quantity('xmomentum',6.0)
72       
73        domain.set_boundary({'exterior' : Dirichlet_boundary([3.0, 8.0, 1.0, 2.0, 4.0])})
74
75
76        domain.distribute_to_vertices_and_edges()
77        domain.update_boundary()
78
79##         print 'In test dirichlet'
80##         print domain.quantities['stage'].vertex_values
81##         print domain.quantities['xmomentum'].vertex_values
82##         print domain.quantities['elevation'].vertex_values
83##         print domain.quantities['height'].vertex_values
84##         print domain.quantities['velocity'].vertex_values
85       
86##         print domain.quantities['stage'].boundary_values
87##         print domain.quantities['xmomentum'].boundary_values
88##         print domain.quantities['elevation'].boundary_values
89##         print domain.quantities['height'].boundary_values
90##         print domain.quantities['velocity'].boundary_values
91
92        assert allclose( domain.quantities['stage'    ].boundary_values, [3.0, 3.0])
93        assert allclose( domain.quantities['xmomentum'].boundary_values, [8.0, 8.0])
94        assert allclose( domain.quantities['elevation'].boundary_values, [1.0, 1.0])
95        assert allclose( domain.quantities['height'   ].boundary_values, [2.0, 2.0])
96        assert allclose( domain.quantities['velocity' ].boundary_values, [4.0, 4.0])
97
98
99    def test_compute_fluxes(self):
100        """
101        Compare shallow_water_domain flux calculation against a previous
102        Python implementation (defined in this file)
103        """
104        domain = Domain(self.points)
105        domain.set_quantity('stage',2.0)
106        domain.set_boundary({'exterior' : Dirichlet_boundary([0.0, 0.0, 0.0, 0.0, 0.0])})
107
108
109        domain.distribute_to_vertices_and_edges()
110        domain.update_boundary()
111
112        stage_ud, xmom_ud = local_compute_fluxes(domain)
113
114        domain.compute_fluxes()
115       
116        print domain.quantities['stage'].vertex_values
117        print domain.quantities['xmomentum'].vertex_values
118        print domain.quantities['elevation'].vertex_values
119        print domain.quantities['height'].vertex_values
120        print domain.quantities['velocity'].vertex_values
121
122        print domain.quantities['stage'].boundary_values
123        print domain.quantities['xmomentum'].boundary_values
124        print domain.quantities['elevation'].boundary_values
125        print domain.quantities['height'].boundary_values
126        print domain.quantities['velocity'].boundary_values
127       
128        print domain.quantities['stage'].explicit_update
129        print domain.quantities['xmomentum'].explicit_update
130        print stage_ud
131        print xmom_ud
132
133        assert allclose( domain.quantities['stage'].explicit_update, stage_ud )
134        assert allclose( domain.quantities['xmomentum'].explicit_update, xmom_ud )
135
136
137    def test_local_flux_function(self):
138        normal = 1.0
139        ql = array([1.0, 2.0],Float)
140        qr = array([1.0, 2.0],Float)
141        zl = 0.0
142        zr = 0.0
143
144        #This assumes h0 = 1.0e-3!!
145        edgeflux, maxspeed = local_flux_function(normal, ql,qr,zl,zr)
146        #print maxspeed
147        #print edgeflux
148       
149        assert allclose(array([2.0, 8.9],Float), edgeflux, rtol=1.0e-005)
150        assert allclose(5.1305, maxspeed, rtol=1.0e-005)
151
152        normal = -1.0
153        ql = array([1.0, 2.0],Float)
154        qr = array([1.0, 2.0],Float)
155        zl = 0.0
156        zr = 0.0
157
158        edgeflux, maxspeed = local_flux_function(normal, ql,qr,zl,zr)
159
160
161        #print maxspeed
162        #print edgeflux       
163       
164        assert allclose(array([-2.0, -8.9],Float), edgeflux, rtol=1.0e-005)
165        assert allclose(5.1305, maxspeed, rtol=1.0e-005)
166
167
168    def test_gravity(self):
169        """
170        Compare shallow_water_domain gravity calculation
171        """
172
173        def slope_one(x):
174            return x
175
176        domain = Domain(self.points)
177        domain.set_quantity('stage',4.0)
178        domain.set_quantity('elevation',slope_one)
179        domain.set_boundary({'exterior' : Reflective_boundary(domain)})
180
181        domain.distribute_to_vertices_and_edges()
182        domain.update_boundary()
183
184        gravity(domain)
185       
186        #print domain.quantities['stage'].vertex_values
187        #print domain.quantities['elevation'].vertex_values
188        #print domain.quantities['xmomentum'].explicit_update       
189
190        assert allclose( array([-34.3, -24.5, -14.7], Float), domain.quantities['xmomentum'].explicit_update )
191
192
193    def test_evolve_first_order(self):
194        """
195        Compare still lake solution for various versions of shallow_water_domain
196        """
197       
198        def slope_square(x):
199            return maximum(4.0-(x-5.0)*(x-5.0), 0.0)
200
201        domain = Domain(self.points2)
202        domain.set_quantity('stage',10.0)
203        domain.set_quantity('elevation',slope_square)
204        domain.set_boundary({'exterior' : Reflective_boundary(domain)})
205
206        domain.default_order = 1
207        domain.set_timestepping_method('euler')
208       
209        yieldstep=0.25
210        finaltime=1.0
211
212        for t in domain.evolve(yieldstep=yieldstep, finaltime=finaltime):
213            pass
214       
215        ##      print domain.quantities['stage'].vertex_values
216        ##      print domain.quantities['elevation'].vertex_values
217        ##        print domain.quantities['xmomentum'].vertex_values
218        ## 
219        ##
220        ##        print domain.quantities['stage'].centroid_values 
221        ##      print domain.quantities['elevation'].centroid_values   
222        ##        print domain.quantities['xmomentum'].centroid_values   
223
224        #assert allclose( 10.0*ones(10), domain.quantities['stage'].centroid_values )
225        #assert allclose( zeros(10), domain.quantities['xmomentum'].centroid_values )
226
227
228##     def test_evolve_euler_second_order_space(self):
229##         """
230##         Compare still lake solution for various versions of shallow_water_domain
231##         """
232
233##         def slope_square(x):
234##             return maximum(4.0-(x-5.0)*(x-5.0), 0.0)
235
236##         domain = Domain(self.points2)
237##         domain.set_quantity('stage',10.0)
238##         domain.set_quantity('elevation',slope_square)
239##         domain.set_boundary({'exterior' : Reflective_boundary(domain)})
240
241##         domain.default_order = 2
242##         domain.set_timestepping_method('rk2')
243##         yieldstep=1.0
244##         finaltime=1.0
245
246##         for t in domain.evolve(yieldstep=yieldstep, finaltime=finaltime):
247##             pass
248       
249##         assert allclose( 10.0*ones(10), domain.quantities['stage'].centroid_values )
250##         assert allclose( zeros(10), domain.quantities['xmomentum'].centroid_values )
251
252##     def test_evolve_second_order_space_time(self):
253##         """
254##         Compare still lake solution for various versions of shallow_water_domain
255##         """
256
257##         def slope_square(x):
258##             return maximum(4.0-(x-5.0)*(x-5.0), 0.0)
259
260##         domain = Domain(self.points2)
261##         domain.set_quantity('stage',10.0)
262##         domain.set_quantity('elevation',slope_square)
263##         domain.set_boundary({'exterior' : Reflective_boundary(domain)})
264
265##         domain.default_order = 2
266##         domain.set_timestepping_method('rk3')
267
268##         yieldstep=1.0
269##         finaltime=1.0
270
271##         for t in domain.evolve(yieldstep=yieldstep, finaltime=finaltime):
272##             pass
273       
274##         assert allclose( 10.0*ones(10), domain.quantities['stage'].centroid_values )
275##         assert allclose( zeros(10), domain.quantities['xmomentum'].centroid_values )
276
277
278
279#==============================================================================
280
281def local_compute_fluxes(domain):
282    """Compute all fluxes and the timestep suitable for all volumes
283    in domain.
284
285    Compute total flux for each conserved quantity using "flux_function"
286
287    Fluxes across each edge are scaled by edgelengths and summed up
288    Resulting flux is then scaled by area and stored in
289    explicit_update for each of the three conserved quantities
290    stage, xmomentum and ymomentum
291
292    The maximal allowable speed computed by the flux_function for each volume
293    is converted to a timestep that must not be exceeded. The minimum of
294    those is computed as the next overall timestep.
295
296    Post conditions:
297      domain.explicit_update is reset to computed flux values
298      domain.timestep is set to the largest step satisfying all volumes.
299    """
300
301    import sys
302    from Numeric import zeros, Float
303
304    N = domain.number_of_elements
305
306    tmp0 = zeros(N,Float)
307    tmp1 = zeros(N,Float)
308
309    #Shortcuts
310    Stage = domain.quantities['stage']
311    Xmom = domain.quantities['xmomentum']
312    #    Ymom = domain.quantities['ymomentum']
313    Bed = domain.quantities['elevation']
314
315    #Arrays
316    #stage = Stage.edge_values
317    #xmom =  Xmom.edge_values
318    #   ymom =  Ymom.edge_values
319    #bed =   Bed.edge_values
320   
321    stage = Stage.vertex_values
322    xmom =  Xmom.vertex_values
323    bed =   Bed.vertex_values
324   
325    #print 'stage edge values', stage
326    #print 'xmom edge values', xmom
327    #print 'bed values', bed
328
329    stage_bdry = Stage.boundary_values
330    xmom_bdry =  Xmom.boundary_values
331    #print 'stage_bdry',stage_bdry
332    #print 'xmom_bdry', xmom_bdry
333    #    ymom_bdry =  Ymom.boundary_values
334
335    #    flux = zeros(3, Float) #Work array for summing up fluxes
336    flux = zeros(2, Float) #Work array for summing up fluxes
337    ql = zeros(2, Float)
338    qr = zeros(2, Float)
339
340    #Loop
341    timestep = float(sys.maxint)
342    enter = True
343    for k in range(N):
344
345        flux[:] = 0.  #Reset work array
346        #for i in range(3):
347        for i in range(2):
348            #Quantities inside volume facing neighbour i
349            #ql[0] = stage[k, i]
350            #ql[1] = xmom[k, i]
351            ql = [stage[k, i], xmom[k, i]]
352            zl = bed[k, i]
353
354            #Quantities at neighbour on nearest face
355            n = domain.neighbours[k,i]
356            if n < 0:
357                m = -n-1 #Convert negative flag to index
358                qr[0] = stage_bdry[m]
359                qr[1] = xmom_bdry[m]
360                zr = zl #Extend bed elevation to boundary
361            else:
362                #m = domain.neighbour_edges[k,i]
363                m = domain.neighbour_vertices[k,i]
364                #print i, ' ' , m
365                #qr = [stage[n, m], xmom[n, m], ymom[n, m]]
366                qr[0] = stage[n, m]
367                qr[1] = xmom[n, m]
368                zr = bed[n, m]
369
370
371            #Outward pointing normal vector
372            normal = domain.normals[k, i]
373       
374            #Flux computation using provided function
375            #edgeflux, max_speed = flux_function(normal, ql, qr, zl, zr)
376            #print 'ql',ql
377            #print 'qr',qr
378           
379
380            edgeflux, max_speed = local_flux_function(normal, ql, qr, zl, zr)
381
382            #print 'edgeflux', edgeflux
383
384            # THIS IS THE LINE TO DEAL WITH LEFT AND RIGHT FLUXES
385            # flux = edgefluxleft - edgefluxright
386            flux -= edgeflux #* domain.edgelengths[k,i]
387            #Update optimal_timestep
388            try:
389                #timestep = min(timestep, 0.5*domain.radii[k]/max_speed)
390                timestep = min(timestep, domain.CFL*0.5*domain.areas[k]/max_speed)
391            except ZeroDivisionError:
392                pass
393
394        #Normalise by area and store for when all conserved
395        #quantities get updated
396        flux /= domain.areas[k]
397
398        #Stage.explicit_update[k] = flux[0]
399        tmp0[k] = flux[0]
400        tmp1[k] = flux[1]
401
402
403    return tmp0, tmp1
404
405
406def local_flux_function(normal, ql, qr, zl, zr):
407    """Compute fluxes between volumes for the shallow water wave equation
408    cast in terms of w = h+z using the 'central scheme' as described in
409
410    Kurganov, Noelle, Petrova. 'Semidiscrete Central-Upwind Schemes For
411    Hyperbolic Conservation Laws and Hamilton-Jacobi Equations'.
412    Siam J. Sci. Comput. Vol. 23, No. 3, pp. 707-740.
413
414    The implemented formula is given in equation (3.15) on page 714
415
416    Conserved quantities w, uh, are stored as elements 0 and 1
417    in the numerical vectors ql an qr.
418
419    Bed elevations zl and zr.
420    """
421
422    from config import g, epsilon, h0
423    from math import sqrt
424    from Numeric import array
425
426    #print 'ql',ql
427
428    #Align momentums with x-axis
429    #q_left  = rotate(ql, normal, direction = 1)
430    #q_right = rotate(qr, normal, direction = 1)
431    q_left = ql
432    q_left[1] = q_left[1]*normal
433    q_right = qr
434    q_right[1] = q_right[1]*normal
435
436    #z = (zl+zr)/2 #Take average of field values
437    z = 0.5*(zl+zr) #Take average of field values
438
439    w_left  = q_left[0]   #w=h+z
440    h_left  = w_left-z
441    uh_left = q_left[1]
442
443    if h_left < epsilon:
444        u_left = 0.0  #Could have been negative
445        h_left = 0.0
446    else:
447        u_left  = uh_left/(h_left +  h0/h_left)
448
449
450    uh_left = u_left*h_left
451
452    w_right  = q_right[0]  #w=h+z
453    h_right  = w_right-z
454    uh_right = q_right[1]
455
456
457    if h_right < epsilon:
458        u_right = 0.0  #Could have been negative
459        h_right = 0.0
460    else:
461        u_right  = uh_right/(h_right + h0/h_right)
462
463    uh_right = u_right*h_right
464   
465    #vh_left  = q_left[2]
466    #vh_right = q_right[2]
467
468    #print h_right
469    #print u_right
470    #print h_left
471    #print u_right
472
473    soundspeed_left  = sqrt(g*h_left)
474    soundspeed_right = sqrt(g*h_right)
475
476    #Maximal wave speed
477    s_max = max(u_left+soundspeed_left, u_right+soundspeed_right, 0)
478   
479    #Minimal wave speed
480    s_min = min(u_left-soundspeed_left, u_right-soundspeed_right, 0)
481   
482    #Flux computation
483
484    #flux_left  = array([u_left*h_left,
485    #                    u_left*uh_left + 0.5*g*h_left**2])
486    #flux_right = array([u_right*h_right,
487    #                    u_right*uh_right + 0.5*g*h_right**2])
488    flux_left  = array([u_left*h_left,
489                        u_left*uh_left + 0.5*g*h_left*h_left])
490    flux_right = array([u_right*h_right,
491                        u_right*uh_right + 0.5*g*h_right*h_right])
492
493    denom = s_max-s_min
494    if denom == 0.0:
495        edgeflux = array([0.0, 0.0])
496        max_speed = 0.0
497    else:
498        edgeflux = (s_max*flux_left - s_min*flux_right)/denom
499        edgeflux += s_max*s_min*(q_right-q_left)/denom
500       
501        edgeflux[1] = edgeflux[1]*normal
502
503        max_speed = max(abs(s_max), abs(s_min))
504
505    return edgeflux, max_speed
506
507
508#-------------------------------------------------------------
509if __name__ == "__main__":
510    suite = unittest.makeSuite(Test_Shallow_Water, 'test')   
511    #suite = unittest.makeSuite(Test_Quantity, 'test_set_values_from_file_using_polygon')
512
513    #suite = unittest.makeSuite(Test_Quantity, 'test_set_vertex_values_using_general_interface_with_subset')
514    #print "restricted test"
515    #suite = unittest.makeSuite(Test_Quantity,'verbose_test_set_values_from_UTM_pts')
516    runner = unittest.TextTestRunner()
517    runner.run(suite)
Note: See TracBrowser for help on using the repository browser.