source: anuga_work/development/pipeflow/test_sww_domain.py @ 7788

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

cleaning up pipeflow directory

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