source: trunk/anuga_work/development/2010-projects/anuga_1d/sww/test_sww_vel_domain.py @ 7884

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

Moving 2010 project

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