source: anuga_work/development/2010-projects/anuga_1d/channel/channel_domain.py @ 7842

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

Adding in a few new files

File size: 15.3 KB
Line 
1"""Class Domain -
21D interval domains for finite-volume computations of
3the shallow water wave equation.
4
5This module contains a specialisation of class Generic_domain from module
6generic_domain.py
7consisting of methods specific to Channel flow using the Shallow Water Wave Equation
8
9This particular modification of the Domain class implements the ability to
10vary the width of the 1D channel that the water flows in. As a result the
11conserved variables are different than previous implementations and so are the
12equations.
13
14U_t + E_x = S
15
16where
17
18U = [A, Q] = [b*h, u*b*h]
19E = [Q, Q^2/A + g*b*h^2/2]
20S represents source terms forcing the system
21(e.g. gravity, boundary_stree, friction, wind stress, ...)
22gravity = -g*b*h*z_x
23boundary_stress = 1/2*g*b_x*h^2
24
25and _t, _x, _y denote the derivative with respect to t, x and y respectiely.
26
27The quantities are
28
29symbol    variable name    explanation
30A         area             Wetted area = b*h
31Q         discharge        flux of water = u*b*h
32x         x                horizontal distance from origin [m]
33z         elevation        elevation of bed on which flow is modelled [m]
34h         height           water height above z [m]
35w         stage            absolute water level, w = z+h [m]
36u                          speed in the x direction [m/s]
37uh        xmomentum        momentum in the x direction [m^2/s]
38b         width            width of channel
39eta                        mannings friction coefficient [to appear]
40nu                         wind stress coefficient [to appear]
41
42The conserved quantities are A, Q
43--------------------------------------------------------------------------
44For details see e.g.
45Christopher Zoppou and Stephen Roberts,
46Catastrophic Collapse of Water Supply Reservoirs in Urban Areas,
47Journal of Hydraulic Engineering, vol. 127, No. 7 July 1999
48
49
50John Jakeman, Ole Nielsen, Stephen Roberts, Duncan Gray, Christopher Zoppou,
51Padarn Wilson, Geoscience Australia, 2008
52"""
53
54
55from anuga_1d.generic.generic_domain import *
56import numpy
57
58
59#Shallow water domain
60class Domain(Generic_domain):
61
62    def __init__(self, coordinates, boundary = None, tagged_elements = None):
63
64        conserved_quantities = ['area', 'discharge']
65        evolved_quantities = ['area', 'discharge', 'elevation', 'height', 'velocity','width','stage']
66        other_quantities = ['friction']
67        Generic_domain.__init__(self,
68                                coordinates = coordinates,
69                                boundary = boundary,
70                                conserved_quantities = conserved_quantities,
71                                evolved_quantities = evolved_quantities,
72                                other_quantities = other_quantities,
73                                tagged_elements = tagged_elements)
74       
75        from anuga_1d.config import minimum_allowed_height, g, h0
76        self.minimum_allowed_height = minimum_allowed_height
77        self.g = g
78        self.h0 = h0
79        self.setstageflag = False
80        self.discontinousb = False
81
82
83        # forcing terms gravity and boundary stress are included in the flux calculation
84        #self.forcing_terms.append(gravity)
85        #self.forcing_terms.append(boundary_stress)
86        #self.forcing_terms.append(manning_friction)
87       
88
89       
90        #Stored output
91        self.store = True
92        self.format = 'sww'
93        self.smooth = True
94
95       
96        #Reduction operation for get_vertex_values
97        from anuga_1d.generic.util import mean
98        self.reduction = mean
99        #self.reduction = min  #Looks better near steep slopes
100
101        self.set_quantities_to_be_stored(['area','discharge'])
102
103        self.__doc__ = 'channel_domain'
104
105        self.check_integrity()
106
107
108    def check_integrity(self):
109
110        #Check that we are solving the shallow water wave equation
111
112        msg = 'First conserved quantity must be "area"'
113        assert self.conserved_quantities[0] == 'area', msg
114        msg = 'Second conserved quantity must be "discharge"'
115        assert self.conserved_quantities[1] == 'discharge', msg
116
117        msg = 'First evolved quantity must be "area"'
118        assert self.evolved_quantities[0] == 'area', msg
119        msg = 'Second evolved quantity must be "discharge"'
120        assert self.evolved_quantities[1] == 'discharge', msg
121        msg = 'Third evolved quantity must be "elevation"'
122        assert self.evolved_quantities[2] == 'elevation', msg
123        msg = 'Fourth evolved quantity must be "height"'
124        assert self.evolved_quantities[3] == 'height', msg
125        msg = 'Fifth evolved quantity must be "velocity"'
126        assert self.evolved_quantities[4] == 'velocity', msg
127        msg = 'Sixth evolved quantity must be "width"'
128        assert self.evolved_quantities[5] == 'width', msg
129        msg = 'Seventh evolved quantity must be "stage"'
130        assert self.evolved_quantities[6] == 'stage', msg
131
132        Generic_domain.check_integrity(self)
133
134    def compute_fluxes(self):
135        #Call correct module function
136        #(either from this module or C-extension)
137        compute_fluxes_channel(self)
138
139    def distribute_to_vertices_and_edges(self):
140        #Call correct module function
141        #(either from this module or C-extension)
142        distribute_to_vertices_and_edges_limit_a_d(self)
143       
144
145#=============== End of Channel Domain ===============================
146
147#-----------------------------------
148# Compute flux definition with channel
149#-----------------------------------
150def compute_fluxes_channel(domain):
151    import sys
152    timestep = float(sys.maxint)
153
154    area       = domain.quantities['area']
155    discharge  = domain.quantities['discharge']
156    bed        = domain.quantities['elevation']
157    height     = domain.quantities['height']
158    velocity   = domain.quantities['velocity']
159    width      = domain.quantities['width']
160
161
162    from anuga_1d.channel.channel_domain_ext import compute_fluxes_channel_ext
163    domain.flux_timestep = compute_fluxes_channel_ext(timestep,domain,area,discharge,bed,height,velocity,width)
164
165#-----------------------------------------------------------------------
166# Distribute to verticies with stage reconstructed and then extrapolated
167#-----------------------------------------------------------------------
168def distribute_to_vertices_and_edges_limit_a_d(domain):
169   
170    #Remove very thin layers of water
171    #protect_against_infinitesimal_and_negative_heights(domain) 
172
173
174
175    import sys
176    from anuga_1d.config import epsilon, h0
177
178   
179    N = domain.number_of_elements
180
181    #Shortcuts
182    area      = domain.quantities['area']
183    discharge = domain.quantities['discharge']
184    bed       = domain.quantities['elevation']
185    height    = domain.quantities['height']
186    velocity  = domain.quantities['velocity']
187    width     = domain.quantities['width']
188    stage     = domain.quantities['stage']
189
190    #Arrays   
191    a_C   = area.centroid_values
192    d_C   = discharge.centroid_values
193    z_C   = bed.centroid_values
194    h_C   = height.centroid_values
195    u_C   = velocity.centroid_values
196    b_C   = width.centroid_values
197    w_C   = stage.centroid_values
198
199    if domain.setstageflag:
200        a_C[:,] = (w_C-z_C)*b_C
201        domain.setstageflag = False
202
203    if domain.discontinousb:
204        domain.quantities['width'].extrapolate_second_order()
205       
206
207    # FIXME (SGR): Replace with numpy.where to speed up code
208    h0 = 1.0e-12
209
210    h_C[:] = numpy.where( (a_C>h0) | (b_C>h0), a_C/(b_C+h0/b_C), 0.0 )
211    u_C[:] = numpy.where( (a_C>h0) | (b_C>h0), d_C/(a_C+h0/a_C), 0.0 )
212
213    a_C[:] = numpy.where( (a_C>h0) | (b_C>h0), a_C, 0.0 )
214    b_C[:] = numpy.where( (a_C>h0) | (b_C>h0), b_C, 0.0 )
215    d_C[:] = numpy.where( (a_C>h0) | (b_C>h0), d_C, 0.0 )
216
217    w_C[:] = h_C + z_C
218   
219#    for i in range(N):
220#
221#        if a_C[i] <= h0:
222#           a_C[i] = 0.0
223#           h_C[i] = 0.0
224#            d_C[i] = 0.0
225#            u_C[i] = 0.0
226#            w_C[i] = z_C[i]
227#        else:
228#
229#            if b_C[i]<=h0:
230#                a_C[i] = 0.0
231#                h_C[i] = 0.0
232#                d_C[i] = 0.0
233#                u_C[i] = 0.0
234#                w_C[i] = z_C[i]
235#
236#            else:
237#                h_C[i] = a_C[i]/(b_C[i]+h0/b_C[i])
238#                w_C[i] = h_C[i]+z_C[i]
239#                u_C[i] = d_C[i]/(a_C[i]+h0/a_C[i])
240               
241
242    bed.extrapolate_second_order()
243   
244    for name in ['velocity','stage']:
245        Q = domain.quantities[name]
246        if domain.order == 1:
247            Q.extrapolate_first_order()
248        elif domain.order == 2:
249            Q.extrapolate_second_order()
250        else:
251            raise 'Unknown order'
252
253
254    a_V  = area.vertex_values
255    w_V  = stage.vertex_values
256    z_V  = bed.vertex_values
257    h_V  = height.vertex_values
258    u_V  = velocity.vertex_values
259    d_V  = discharge.vertex_values
260    b_V  = width.vertex_values
261
262   
263    #FIXME (SGR): replace with numpy.where
264
265
266    h_V[:] = w_V-z_V
267
268    w_V[:] = numpy.where(h_V > h0, w_V, z_V)
269    h_V[:] = numpy.where(h_V > h0, h_V, 0.0)
270    a_V[:] = b_V*h_V
271    d_V[:] = u_V*a_V
272
273#    for i in range(len(h_C)):
274#        for j in range(2):
275#            h_V[i,j] = w_V[i,j]-z_V[i,j]
276#            if h_V[i,j]<h0:
277#                h_V[i,j]=0
278#                w_V[i,j]=z_V[i,j]
279#            a_V[i,j] = b_V[i,j]*h_V[i,j]
280#            d_V[i,j]=u_V[i,j]*a_V[i,j]
281   
282    return
283
284
285#--------------------------------------------------------
286#Boundaries - specific to the channel_domain
287#--------------------------------------------------------
288class Reflective_boundary(Boundary):
289    """Reflective boundary returns same conserved quantities as
290    those present in its neighbour volume but reflected.
291
292    This class is specific to the shallow water equation as it
293    works with the momentum quantities assumed to be the second
294    and third conserved quantities.
295    """
296
297    def __init__(self, domain = None):
298        Boundary.__init__(self)
299
300        if domain is None:
301            msg = 'Domain must be specified for reflective boundary'
302            raise msg
303
304        #Handy shorthands
305        self.normals  = domain.normals
306        self.area     = domain.quantities['area'].vertex_values
307        self.discharge     = domain.quantities['discharge'].vertex_values
308        self.bed      = domain.quantities['elevation'].vertex_values
309        self.height   = domain.quantities['height'].vertex_values
310        self.velocity = domain.quantities['velocity'].vertex_values
311        self.width    = domain.quantities['width'].vertex_values
312        self.stage    = domain.quantities['stage'].vertex_values
313
314        self.evolved_quantities = numpy.zeros(7, numpy.float)
315
316    def __repr__(self):
317        return 'Reflective_boundary'
318
319
320    def evaluate(self, vol_id, edge_id):
321        """Reflective boundaries reverses the outward momentum
322        of the volume they serve.
323        """
324
325        q = self.evolved_quantities
326        q[0] =  self.area[vol_id, edge_id]
327        q[1] = -self.discharge[vol_id, edge_id]
328        q[2] =  self.bed[vol_id, edge_id]
329        q[3] =  self.height[vol_id, edge_id]
330        q[4] = -self.velocity[vol_id, edge_id]
331        q[5] =  self.width[vol_id,edge_id]
332        q[6] =  self.stage[vol_id,edge_id]
333
334        return q
335
336class Dirichlet_boundary(Boundary):
337    """Dirichlet boundary returns constant values for the
338    conserved quantities
339if k>5 and k<15:
340            print discharge_ud[k],-g*zx*avg_h*avg_b
341        discharge_ud[k] +=-g*zx*avg_h*avg_b    """
342
343
344    def __init__(self, evolved_quantities=None):
345        Boundary.__init__(self)
346
347        if evolved_quantities is None:
348            msg = 'Must specify one value for each evolved quantity'
349            raise msg
350
351        assert len(evolved_quantities) == 7
352
353        self.evolved_quantities=numpy.array(evolved_quantities,numpy.float)
354
355    def __repr__(self):
356        return 'Dirichlet boundary (%s)' %self.evolved_quantities
357
358    def evaluate(self, vol_id=None, edge_id=None):
359        return self.evolved_quantities
360
361
362#----------------------------
363#Standard forcing terms:
364#---------------------------
365def gravity(domain):
366    """Apply gravitational pull in the presence of bed slope
367    """
368
369    from util import gradient
370    from Numeric import zeros, Float, array, sum
371
372
373
374    Area      = domain.quantities['area']
375    Discharge  = domain.quantities['discharge']
376    Elevation = domain.quantities['elevation']
377    Height    = domain.quantities['height']
378    Width     = domain.quantities['width']
379
380    discharge_ud  = Discharge.explicit_update
381
382
383       
384    h = Height.vertex_values
385    b = Width.vertex_values
386    a = Area.vertex_values
387    z = Elevation.vertex_values
388
389    x = domain.get_vertex_coordinates()
390    g = domain.g
391    for k in range(domain.number_of_elements):
392        avg_h = 0.5*(h[k,0] + h[k,1])
393        avg_b = 0.5*(b[k,0] + b[k,1])
394       
395        #Compute bed slope
396        x0, x1 = x[k,:]
397        z0, z1 = z[k,:]
398        zx = gradient(x0, x1, z0, z1)
399       
400        #Update momentum (explicit update is reset to source values)
401        discharge_ud[k]+= -g*zx*avg_h*avg_b
402     
403
404def boundary_stress(domain):
405
406    from util import gradient
407    from Numeric import zeros, Float, array, sum
408
409    Area     = domain.quantities['area']
410    Discharge  = domain.quantities['discharge']
411    Elevation = domain.quantities['elevation']
412    Height    = domain.quantities['height']
413    Width     = domain.quantities['width']
414
415    discharge_ud  = Discharge.explicit_update
416 
417    h = Height.vertex_values
418    b = Width.vertex_values
419    a = Area.vertex_values
420    z = Elevation.vertex_values
421
422    x = domain.get_vertex_coordinates()
423    g = domain.g
424
425    for k in range(domain.number_of_elements):
426        avg_h = 0.5*(h[k,0] + h[k,1])
427       
428
429        #Compute bed slope
430        x0, x1 = x[k,:]
431        b0, b1 = b[k,:]
432        bx = gradient(x0, x1, b0, b1)
433       
434        #Update momentum (explicit update is reset to source values)
435        discharge_ud[k] += 0.5*g*bx*avg_h*avg_h
436        #stage_ud[k] = 0.0
437 
438 
439def manning_friction(domain):
440    """Apply (Manning) friction to water momentum
441    """
442
443    from math import sqrt
444
445    w = domain.quantities['stage'].centroid_values
446    z = domain.quantities['elevation'].centroid_values
447    h = w-z
448
449    uh = domain.quantities['xmomentum'].centroid_values
450    eta = domain.quantities['friction'].centroid_values
451
452    xmom_update = domain.quantities['xmomentum'].semi_implicit_update
453
454    N = domain.number_of_elements
455    eps = domain.minimum_allowed_height
456    g = domain.g
457
458    for k in range(N):
459        if eta[k] >= eps:
460            if h[k] >= eps:
461                #S = -g * eta[k]**2 * sqrt((uh[k]**2 + vh[k]**2))
462                S = -g * eta[k]**2 * uh[k]
463                S /= h[k]**(7.0/3)
464
465                #Update momentum
466                xmom_update[k] += S*uh[k]
467                #ymom_update[k] += S*vh[k]
468
469def linear_friction(domain):
470    """Apply linear friction to water momentum
471
472    Assumes quantity: 'linear_friction' to be present
473    """
474
475    from math import sqrt
476
477    w = domain.quantities['stage'].centroid_values
478    z = domain.quantities['elevation'].centroid_values
479    h = w-z
480
481    uh = domain.quantities['xmomentum'].centroid_values
482    tau = domain.quantities['linear_friction'].centroid_values
483
484    xmom_update = domain.quantities['xmomentum'].semi_implicit_update
485
486    N = domain.number_of_elements
487    eps = domain.minimum_allowed_height
488
489    for k in range(N):
490        if tau[k] >= eps:
491            if h[k] >= eps:
492                S = -tau[k]/h[k]
493
494                #Update momentum
495                xmom_update[k] += S*uh[k]
496
497
498
499
500def linearb(domain):
501
502    bC = domain.quantities['width'].vertex_values
503   
504    for i in range(len(bC)-1):
505        temp= 0.5*(bC[i,1]+bC[i+1,0])
506        bC[i,1]=temp
507        bC[i+1,0]=temp
508
509
510
Note: See TracBrowser for help on using the repository browser.