source: branches/anuga_1_1/anuga_work/development/anuga_1d/test_shallow_water_domain_suggestion1.py @ 7963

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

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