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

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

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