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

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

Moving 2010 project

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 anuga_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.