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

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

Added in minimum height

File size: 9.7 KB
Line 
1#!/usr/bin/env python
2
3import unittest
4from math import sqrt, pi
5
6
7from shallow_water_domain import *
8from Numeric import allclose, array, ones, Float
9
10
11class Test_Shallow_Water(unittest.TestCase):
12    def setUp(self):
13        self.points = [0.0, 1.0, 2.0, 3.0]
14        self.vertex_values = [[1.0,2.0],[4.0,5.0],[-1.0,2.0]]
15       
16    def tearDown(self):
17        pass
18        #print "  Tearing down"
19
20
21    def test_creation(self):
22        domain = Domain(self.points)
23        assert allclose(domain.centroids, [0.5, 1.5, 2.5])
24
25    def test_compute_fluxes(self):
26        """
27        Compare shallow_water_domain flux calculation against a previous
28        Python implementation (defined in this file)
29        """
30        domain = Domain(self.points)
31        domain.set_quantity('stage',2.0)
32        domain.set_boundary({'exterior' : Reflective_boundary(domain)})
33
34        stage_ud, xmom_ud = compute_fluxes_python(domain)
35       
36        domain.compute_fluxes()
37
38        #print doamin.quantities['xmomentum'].explicit_update
39        #print compute_fluxes_python(domain)
40
41        assert allclose( domain.quantities['stage'].explicit_update, stage_ud )
42        assert allclose( domain.quantities['xmomentum'].explicit_update, xmom_ud )
43
44
45    def test_local_flux_function(self):
46        normal = 1.0
47        ql = array([1.0, 2.0],Float)
48        qr = array([1.0, 2.0],Float)
49        zl = 0.0
50        zr = 0.0
51
52        #This assumes h0 = 1.0e-3!!
53        edgeflux, maxspeed = flux_function(normal, ql,qr,zl,zr)
54        #print maxspeed
55        #print edgeflux
56       
57        assert allclose(array([1.998002, 8.89201198],Float), edgeflux)
58        assert allclose(5.1284971665, maxspeed)
59
60        normal = -1.0
61        ql = array([1.0, 2.0],Float)
62        qr = array([1.0, 2.0],Float)
63        zl = 0.0
64        zr = 0.0
65
66        edgeflux, maxspeed = flux_function(normal, ql,qr,zl,zr)
67
68
69        #print maxspeed
70        #print edgeflux       
71       
72        assert allclose(array([-1.998002, -8.89201198],Float), edgeflux)
73        assert allclose(5.1284971665, maxspeed)
74
75    def test_domain_flux_function(self):
76        normal = 1.0
77        ql = array([1.0, 2.0],Float)
78        qr = array([1.0, 2.0],Float)
79        zl = 0.0
80        zr = 0.0
81
82        edgeflux, maxspeed = flux_function(normal, ql,qr,zl,zr)
83
84        #print edgeflux
85
86        from shallow_water_domain import flux_function as domain_flux_function
87
88        domainedgeflux, domainmaxspeed = domain_flux_function(normal, ql,qr,zl,zr)
89
90        #print domainedgeflux
91
92        assert allclose(domainedgeflux, edgeflux)
93        assert allclose(domainmaxspeed, maxspeed)
94
95
96
97#==============================================================================
98
99def compute_fluxes_python(domain):
100    """Compute all fluxes and the timestep suitable for all volumes
101    in domain.
102
103    Compute total flux for each conserved quantity using "flux_function"
104
105    Fluxes across each edge are scaled by edgelengths and summed up
106    Resulting flux is then scaled by area and stored in
107    explicit_update for each of the three conserved quantities
108    stage, xmomentum and ymomentum
109
110    The maximal allowable speed computed by the flux_function for each volume
111    is converted to a timestep that must not be exceeded. The minimum of
112    those is computed as the next overall timestep.
113
114    Post conditions:
115      domain.explicit_update is reset to computed flux values
116      domain.timestep is set to the largest step satisfying all volumes.
117    """
118
119    import sys
120    from Numeric import zeros, Float
121
122    N = domain.number_of_elements
123
124    tmp0 = zeros((N,),Float)
125    tmp1 = zeros((N,),Float)
126
127    #Shortcuts
128    Stage = domain.quantities['stage']
129    Xmom = domain.quantities['xmomentum']
130#    Ymom = domain.quantities['ymomentum']
131    Bed = domain.quantities['elevation']
132
133    #Arrays
134    #stage = Stage.edge_values
135    #xmom =  Xmom.edge_values
136 #   ymom =  Ymom.edge_values
137    #bed =   Bed.edge_values
138   
139    stage = Stage.vertex_values
140    xmom =  Xmom.vertex_values
141    bed =   Bed.vertex_values
142   
143    #print 'stage edge values', stage
144    #print 'xmom edge values', xmom
145    #print 'bed values', bed
146
147    stage_bdry = Stage.boundary_values
148    xmom_bdry =  Xmom.boundary_values
149    #print 'stage_bdry',stage_bdry
150    #print 'xmom_bdry', xmom_bdry
151#    ymom_bdry =  Ymom.boundary_values
152
153#    flux = zeros(3, Float) #Work array for summing up fluxes
154    flux = zeros(2, Float) #Work array for summing up fluxes
155    ql = zeros(2, Float)
156    qr = zeros(2, Float)
157
158    #Loop
159    timestep = float(sys.maxint)
160    enter = True
161    for k in range(N):
162
163        flux[:] = 0.  #Reset work array
164        #for i in range(3):
165        for i in range(2):
166            #Quantities inside volume facing neighbour i
167            #ql[0] = stage[k, i]
168            #ql[1] = xmom[k, i]
169            ql = [stage[k, i], xmom[k, i]]
170            zl = bed[k, i]
171
172            #Quantities at neighbour on nearest face
173            n = domain.neighbours[k,i]
174            if n < 0:
175                m = -n-1 #Convert negative flag to index
176                qr[0] = stage_bdry[m]
177                qr[1] = xmom_bdry[m]
178                zr = zl #Extend bed elevation to boundary
179            else:
180                #m = domain.neighbour_edges[k,i]
181                m = domain.neighbour_vertices[k,i]
182                #print i, ' ' , m
183                #qr = [stage[n, m], xmom[n, m], ymom[n, m]]
184                qr[0] = stage[n, m]
185                qr[1] = xmom[n, m]
186                zr = bed[n, m]
187
188
189            #Outward pointing normal vector
190            normal = domain.normals[k, i]
191       
192            #Flux computation using provided function
193            #edgeflux, max_speed = flux_function(normal, ql, qr, zl, zr)
194            #print 'ql',ql
195            #print 'qr',qr
196           
197
198            edgeflux, max_speed = flux_function(normal, ql, qr, zl, zr)
199
200            #print 'edgeflux', edgeflux
201
202            # THIS IS THE LINE TO DEAL WITH LEFT AND RIGHT FLUXES
203            # flux = edgefluxleft - edgefluxright
204            flux -= edgeflux #* domain.edgelengths[k,i]
205            #Update optimal_timestep
206            try:
207                #timestep = min(timestep, 0.5*domain.radii[k]/max_speed)
208                timestep = min(timestep, domain.cfl*0.5*domain.areas[k]/max_speed)
209            except ZeroDivisionError:
210                pass
211
212        #Normalise by area and store for when all conserved
213        #quantities get updated
214        flux /= domain.areas[k]
215
216        #Stage.explicit_update[k] = flux[0]
217        tmp0[k] = flux[0]
218        tmp1[k] = flux[1]
219
220
221    return tmp0, tmp1
222
223
224def flux_function(normal, ql, qr, zl, zr):
225    """Compute fluxes between volumes for the shallow water wave equation
226    cast in terms of w = h+z using the 'central scheme' as described in
227
228    Kurganov, Noelle, Petrova. 'Semidiscrete Central-Upwind Schemes For
229    Hyperbolic Conservation Laws and Hamilton-Jacobi Equations'.
230    Siam J. Sci. Comput. Vol. 23, No. 3, pp. 707-740.
231
232    The implemented formula is given in equation (3.15) on page 714
233
234    Conserved quantities w, uh, are stored as elements 0 and 1
235    in the numerical vectors ql an qr.
236
237    Bed elevations zl and zr.
238    """
239
240    from config import g, epsilon, h0
241    from math import sqrt
242    from Numeric import array
243
244    #print 'ql',ql
245
246    #Align momentums with x-axis
247    #q_left  = rotate(ql, normal, direction = 1)
248    #q_right = rotate(qr, normal, direction = 1)
249    q_left = ql
250    q_left[1] = q_left[1]*normal
251    q_right = qr
252    q_right[1] = q_right[1]*normal
253
254    #z = (zl+zr)/2 #Take average of field values
255    z = 0.5*(zl+zr) #Take average of field values
256
257    w_left  = q_left[0]   #w=h+z
258    h_left  = w_left-z
259    uh_left = q_left[1]
260
261    if h_left < epsilon:
262        u_left = 0.0  #Could have been negative
263        h_left = 0.0
264    else:
265        u_left  = uh_left/(h_left +  h0/h_left)
266
267
268    uh_left = u_left*h_left
269
270    w_right  = q_right[0]  #w=h+z
271    h_right  = w_right-z
272    uh_right = q_right[1]
273
274
275    if h_right < epsilon:
276        u_right = 0.0  #Could have been negative
277        h_right = 0.0
278    else:
279        u_right  = uh_right/(h_right + h0/h_right)
280
281    uh_right = u_right*h_right
282   
283    #vh_left  = q_left[2]
284    #vh_right = q_right[2]
285
286    #print h_right
287    #print u_right
288    #print h_left
289    #print u_right
290
291    soundspeed_left  = sqrt(g*h_left)
292    soundspeed_right = sqrt(g*h_right)
293
294    #Maximal wave speed
295    s_max = max(u_left+soundspeed_left, u_right+soundspeed_right, 0)
296   
297    #Minimal wave speed
298    s_min = min(u_left-soundspeed_left, u_right-soundspeed_right, 0)
299   
300    #Flux computation
301
302    #flux_left  = array([u_left*h_left,
303    #                    u_left*uh_left + 0.5*g*h_left**2])
304    #flux_right = array([u_right*h_right,
305    #                    u_right*uh_right + 0.5*g*h_right**2])
306    flux_left  = array([u_left*h_left,
307                        u_left*uh_left + 0.5*g*h_left*h_left])
308    flux_right = array([u_right*h_right,
309                        u_right*uh_right + 0.5*g*h_right*h_right])
310
311    denom = s_max-s_min
312    if denom == 0.0:
313        edgeflux = array([0.0, 0.0])
314        max_speed = 0.0
315    else:
316        edgeflux = (s_max*flux_left - s_min*flux_right)/denom
317        edgeflux += s_max*s_min*(q_right-q_left)/denom
318       
319        edgeflux[1] = edgeflux[1]*normal
320
321        max_speed = max(abs(s_max), abs(s_min))
322
323    return edgeflux, max_speed
324
325
326#-------------------------------------------------------------
327if __name__ == "__main__":
328    suite = unittest.makeSuite(Test_Shallow_Water, 'test')   
329    #suite = unittest.makeSuite(Test_Quantity, 'test_set_values_from_file_using_polygon')
330
331    #suite = unittest.makeSuite(Test_Quantity, 'test_set_vertex_values_using_general_interface_with_subset')
332    #print "restricted test"
333    #suite = unittest.makeSuite(Test_Quantity,'verbose_test_set_values_from_UTM_pts')
334    runner = unittest.TextTestRunner()
335    runner.run(suite)
Note: See TracBrowser for help on using the repository browser.