source: anuga_work/development/2010-projects/anuga_1d/sww/sww_python.py @ 7839

Last change on this file since 7839 was 7839, checked in by steve, 12 years ago

Changing name of 1d projects so that it will be easy to moveto the trunk

File size: 6.1 KB
Line 
1#! /usr/bin/python
2
3
4
5__author__="Stephen Roberts"
6__date__ ="$05/06/2010 4:54:58 PM$"
7
8
9def compute_fluxes(domain):
10    """Compute all fluxes and the timestep suitable for all volumes
11    in domain.
12
13    Compute total flux for each conserved quantity using "flux_function"
14
15    Fluxes across each edge are scaled by edgelengths and summed up
16    Resulting flux is then scaled by area and stored in
17    explicit_update for each of the three conserved quantities
18    stage, xmomentum and ymomentum
19
20    The maximal allowable speed computed by the flux_function for each volume
21    is converted to a timestep that must not be exceeded. The minimum of
22    those is computed as the next overall timestep.
23
24    Post conditions:
25      domain.explicit_update is reset to computed flux values
26      domain.timestep is set to the largest step satisfying all volumes.
27    """
28
29    import sys
30    import numpy
31
32    N = domain.number_of_elements
33
34    tmp0 = numpy.zeros(N,numpy.float)
35    tmp1 = numpy.zeros(N,numpy.float)
36
37    #Shortcuts
38    Stage = domain.quantities['stage']
39    Xmom = domain.quantities['xmomentum']
40    Bed = domain.quantities['elevation']
41
42    #Arrays
43    stage = Stage.vertex_values
44    xmom =  Xmom.vertex_values
45    bed =   Bed.vertex_values
46
47
48    stage_bdry = Stage.boundary_values
49    xmom_bdry =  Xmom.boundary_values
50
51    flux = numpy.zeros(2, numpy.float) #Work numpy.array for summing up fluxes
52    ql = numpy.zeros(2, numpy.float)
53    qr = numpy.zeros(2, numpy.float)
54
55    #Loop
56    timestep = float(sys.maxint)
57
58    for k in range(N):
59
60        flux[:,] = 0.  #Reset work numpy.array
61        #for i in range(3):
62        for i in range(2):
63            #Quantities inside volume facing neighbour i
64            #ql[0] = stage[k, i]
65            #ql[1] = xmom[k, i]
66            ql = [stage[k, i], xmom[k, i]]
67            zl = bed[k, i]
68
69            #Quantities at neighbour on nearest face
70            n = domain.neighbours[k,i]
71            if n < 0:
72                m = -n-1 #Convert negative flag to index
73                qr[0] = stage_bdry[m]
74                qr[1] = xmom_bdry[m]
75                zr = zl #Extend bed elevation to boundary
76            else:
77                #m = domain.neighbour_edges[k,i]
78                m = domain.neighbour_vertices[k,i]
79                #print i, ' ' , m
80                #qr = [stage[n, m], xmom[n, m], ymom[n, m]]
81                qr[0] = stage[n, m]
82                qr[1] = xmom[n, m]
83                zr = bed[n, m]
84
85
86            #Outward pointing normal vector
87            normal = domain.normals[k, i]
88
89            #Flux computation using provided function
90
91
92            edgeflux, max_speed = flux_function(normal, ql, qr, zl, zr)
93
94            #print 'edgeflux', edgeflux
95
96            # THIS IS THE LINE TO DEAL WITH LEFT AND RIGHT FLUXES
97            # flux = edgefluxleft - edgefluxright
98            flux -= edgeflux
99            #Update optimal_timestep
100            try:
101                #timestep = min(timestep, 0.5*domain.radii[k]/max_speed)
102                timestep = min(timestep, domain.CFL*0.5*domain.areas[k]/max_speed)
103            except ZeroDivisionError:
104                pass
105
106        #Normalise by area and store for when all conserved
107        #quantities get updated
108        flux /= domain.areas[k]
109
110        #Stage.explicit_update[k] = flux[0]
111        tmp0[k] = flux[0]
112        tmp1[k] = flux[1]
113
114
115    return tmp0, tmp1
116
117
118def flux_function(normal, ql, qr, zl, zr):
119    """Compute fluxes between volumes for the shallow water wave equation
120    cast in terms of w = h+z using the 'central scheme' as described in
121
122    Kurganov, Noelle, Petrova. 'Semidiscrete Central-Upwind Schemes For
123    Hyperbolic Conservation Laws and Hamilton-Jacobi Equations'.
124    Siam J. Sci. Comput. Vol. 23, No. 3, pp. 707-740.
125
126    The implemented formula is given in equation (3.15) on page 714
127
128    Conserved quantities w, uh, are stored as elements 0 and 1
129    in the numerical vectors ql an qr.
130
131    Bed elevations zl and zr.
132    """
133
134    from flow_1d.config import g, epsilon, h0
135    from math import sqrt
136    import numpy
137
138    #print 'ql',ql
139
140    #Align momentums with x-axis
141    #q_left  = rotate(ql, normal, direction = 1)
142    #q_right = rotate(qr, normal, direction = 1)
143    q_left = ql
144    q_left[1] = q_left[1]*normal
145    q_right = qr
146    q_right[1] = q_right[1]*normal
147
148    #z = (zl+zr)/2 #Take average of field values
149    z = 0.5*(zl+zr) #Take average of field values
150
151    w_left  = q_left[0]   #w=h+z
152    h_left  = w_left-z
153    uh_left = q_left[1]
154
155    if h_left < epsilon:
156        u_left = 0.0  #Could have been negative
157        h_left = 0.0
158    else:
159        u_left  = uh_left/(h_left +  h0/h_left)
160
161
162    uh_left = u_left*h_left
163
164    w_right  = q_right[0]  #w=h+z
165    h_right  = w_right-z
166    uh_right = q_right[1]
167
168
169    if h_right < epsilon:
170        u_right = 0.0  #Could have been negative
171        h_right = 0.0
172    else:
173        u_right  = uh_right/(h_right + h0/h_right)
174
175    uh_right = u_right*h_right
176
177    #vh_left  = q_left[2]
178    #vh_right = q_right[2]
179
180    #print h_right
181    #print u_right
182    #print h_left
183    #print u_right
184
185    soundspeed_left  = sqrt(g*h_left)
186    soundspeed_right = sqrt(g*h_right)
187
188    #Maximal wave speed
189    s_max = max(u_left+soundspeed_left, u_right+soundspeed_right, 0)
190
191    #Minimal wave speed
192    s_min = min(u_left-soundspeed_left, u_right-soundspeed_right, 0)
193
194    #Flux computation
195
196    #flux_left  = numpy.array([u_left*h_left,
197    #                    u_left*uh_left + 0.5*g*h_left**2])
198    #flux_right = numpy.array([u_right*h_right,
199    #                    u_right*uh_right + 0.5*g*h_right**2])
200    flux_left  = numpy.array([u_left*h_left,
201                        u_left*uh_left + 0.5*g*h_left*h_left])
202    flux_right = numpy.array([u_right*h_right,
203                        u_right*uh_right + 0.5*g*h_right*h_right])
204
205    denom = s_max-s_min
206    if denom == 0.0:
207        edgeflux = numpy.array([0.0, 0.0])
208        max_speed = 0.0
209    else:
210        edgeflux = (s_max*flux_left - s_min*flux_right)/denom
211        edgeflux += s_max*s_min*(q_right-q_left)/denom
212
213        edgeflux[1] = edgeflux[1]*normal
214
215        max_speed = max(abs(s_max), abs(s_min))
216
217    return edgeflux, max_speed
218
219
220if __name__ == "__main__":
221    print "Hello World";
Note: See TracBrowser for help on using the repository browser.