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

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