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

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

anuga_1d works with numeric on 64bit machine

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