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

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

Moving 2010 project

File size: 14.4 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.base.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.base.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        width.extrapolate_second_order()
205       
206
207    h0 = 1.0e-12
208
209    h_C[:] = numpy.where( (a_C>h0) | (b_C>h0), a_C/(b_C+h0/b_C), 0.0 )
210    u_C[:] = numpy.where( (a_C>h0) | (b_C>h0), d_C/(a_C+h0/a_C), 0.0 )
211
212    a_C[:] = numpy.where( (a_C>h0) | (b_C>h0), a_C, 0.0 )
213    b_C[:] = numpy.where( (a_C>h0) | (b_C>h0), b_C, 0.0 )
214    d_C[:] = numpy.where( (a_C>h0) | (b_C>h0), d_C, 0.0 )
215
216    w_C[:] = h_C + z_C
217   
218               
219
220    bed.extrapolate_second_order()
221   
222    for name in ['velocity','stage']:
223        Q = domain.quantities[name]
224        if domain.order == 1:
225            Q.extrapolate_first_order()
226        elif domain.order == 2:
227            Q.extrapolate_second_order()
228        else:
229            raise 'Unknown order'
230
231
232    a_V  = area.vertex_values
233    w_V  = stage.vertex_values
234    z_V  = bed.vertex_values
235    h_V  = height.vertex_values
236    u_V  = velocity.vertex_values
237    d_V  = discharge.vertex_values
238    b_V  = width.vertex_values
239
240
241    h_V[:] = w_V-z_V
242
243    w_V[:] = numpy.where(h_V > h0, w_V, z_V)
244    h_V[:] = numpy.where(h_V > h0, h_V, 0.0)
245    a_V[:] = b_V*h_V
246    d_V[:] = u_V*a_V
247
248   
249    return
250
251
252#--------------------------------------------------------
253#Boundaries - specific to the channel_domain
254#--------------------------------------------------------
255class Reflective_boundary(Boundary):
256    """Reflective boundary returns same conserved quantities as
257    those present in its neighbour volume but reflected.
258
259    This class is specific to the shallow water equation as it
260    works with the momentum quantities assumed to be the second
261    and third conserved quantities.
262    """
263
264    def __init__(self, domain = None):
265        Boundary.__init__(self)
266
267        if domain is None:
268            msg = 'Domain must be specified for reflective boundary'
269            raise msg
270
271        #Handy shorthands
272        self.normals  = domain.normals
273        self.area     = domain.quantities['area'].vertex_values
274        self.discharge     = domain.quantities['discharge'].vertex_values
275        self.bed      = domain.quantities['elevation'].vertex_values
276        self.height   = domain.quantities['height'].vertex_values
277        self.velocity = domain.quantities['velocity'].vertex_values
278        self.width    = domain.quantities['width'].vertex_values
279        self.stage    = domain.quantities['stage'].vertex_values
280
281        self.evolved_quantities = numpy.zeros(7, numpy.float)
282
283    def __repr__(self):
284        return 'Reflective_boundary'
285
286
287    def evaluate(self, vol_id, edge_id):
288        """Reflective boundaries reverses the outward momentum
289        of the volume they serve.
290        """
291
292        q = self.evolved_quantities
293        q[0] =  self.area[vol_id, edge_id]
294        q[1] = -self.discharge[vol_id, edge_id]
295        q[2] =  self.bed[vol_id, edge_id]
296        q[3] =  self.height[vol_id, edge_id]
297        q[4] = -self.velocity[vol_id, edge_id]
298        q[5] =  self.width[vol_id,edge_id]
299        q[6] =  self.stage[vol_id,edge_id]
300
301        return q
302
303class Dirichlet_boundary(Boundary):
304    """Dirichlet boundary returns constant values for the
305    conserved quantities
306if k>5 and k<15:
307            print discharge_ud[k],-g*zx*avg_h*avg_b
308        discharge_ud[k] +=-g*zx*avg_h*avg_b    """
309
310
311    def __init__(self, evolved_quantities=None):
312        Boundary.__init__(self)
313
314        if evolved_quantities is None:
315            msg = 'Must specify one value for each evolved quantity'
316            raise msg
317
318        assert len(evolved_quantities) == 7
319
320        self.evolved_quantities=numpy.array(evolved_quantities,numpy.float)
321
322    def __repr__(self):
323        return 'Dirichlet boundary (%s)' %self.evolved_quantities
324
325    def evaluate(self, vol_id=None, edge_id=None):
326        return self.evolved_quantities
327
328
329#----------------------------
330#Standard forcing terms:
331#---------------------------
332def gravity(domain):
333    """Apply gravitational pull in the presence of bed slope
334    """
335
336    from util import gradient
337    from Numeric import zeros, Float, array, sum
338
339
340
341    Area      = domain.quantities['area']
342    Discharge  = domain.quantities['discharge']
343    Elevation = domain.quantities['elevation']
344    Height    = domain.quantities['height']
345    Width     = domain.quantities['width']
346
347    discharge_ud  = Discharge.explicit_update
348
349
350       
351    h = Height.vertex_values
352    b = Width.vertex_values
353    a = Area.vertex_values
354    z = Elevation.vertex_values
355
356    x = domain.get_vertex_coordinates()
357    g = domain.g
358    for k in range(domain.number_of_elements):
359        avg_h = 0.5*(h[k,0] + h[k,1])
360        avg_b = 0.5*(b[k,0] + b[k,1])
361       
362        #Compute bed slope
363        x0, x1 = x[k,:]
364        z0, z1 = z[k,:]
365        zx = gradient(x0, x1, z0, z1)
366       
367        #Update momentum (explicit update is reset to source values)
368        discharge_ud[k]+= -g*zx*avg_h*avg_b
369     
370
371def boundary_stress(domain):
372
373    from util import gradient
374    from Numeric import zeros, Float, array, sum
375
376    Area     = domain.quantities['area']
377    Discharge  = domain.quantities['discharge']
378    Elevation = domain.quantities['elevation']
379    Height    = domain.quantities['height']
380    Width     = domain.quantities['width']
381
382    discharge_ud  = Discharge.explicit_update
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
392    for k in range(domain.number_of_elements):
393        avg_h = 0.5*(h[k,0] + h[k,1])
394       
395
396        #Compute bed slope
397        x0, x1 = x[k,:]
398        b0, b1 = b[k,:]
399        bx = gradient(x0, x1, b0, b1)
400       
401        #Update momentum (explicit update is reset to source values)
402        discharge_ud[k] += 0.5*g*bx*avg_h*avg_h
403        #stage_ud[k] = 0.0
404 
405 
406def manning_friction(domain):
407    """Apply (Manning) friction to water momentum
408    """
409
410    from math import sqrt
411
412    w = domain.quantities['stage'].centroid_values
413    z = domain.quantities['elevation'].centroid_values
414    h = w-z
415
416    uh = domain.quantities['xmomentum'].centroid_values
417    eta = domain.quantities['friction'].centroid_values
418
419    xmom_update = domain.quantities['xmomentum'].semi_implicit_update
420
421    N = domain.number_of_elements
422    eps = domain.minimum_allowed_height
423    g = domain.g
424
425    for k in range(N):
426        if eta[k] >= eps:
427            if h[k] >= eps:
428                #S = -g * eta[k]**2 * sqrt((uh[k]**2 + vh[k]**2))
429                S = -g * eta[k]**2 * uh[k]
430                S /= h[k]**(7.0/3)
431
432                #Update momentum
433                xmom_update[k] += S*uh[k]
434                #ymom_update[k] += S*vh[k]
435
436def linear_friction(domain):
437    """Apply linear friction to water momentum
438
439    Assumes quantity: 'linear_friction' to be present
440    """
441
442    from math import sqrt
443
444    w = domain.quantities['stage'].centroid_values
445    z = domain.quantities['elevation'].centroid_values
446    h = w-z
447
448    uh = domain.quantities['xmomentum'].centroid_values
449    tau = domain.quantities['linear_friction'].centroid_values
450
451    xmom_update = domain.quantities['xmomentum'].semi_implicit_update
452
453    N = domain.number_of_elements
454    eps = domain.minimum_allowed_height
455
456    for k in range(N):
457        if tau[k] >= eps:
458            if h[k] >= eps:
459                S = -tau[k]/h[k]
460
461                #Update momentum
462                xmom_update[k] += S*uh[k]
463
464
465
466
467def linearb(domain):
468
469    bC = domain.quantities['width'].vertex_values
470   
471    for i in range(len(bC)-1):
472        temp= 0.5*(bC[i,1]+bC[i+1,0])
473        bC[i,1]=temp
474        bC[i+1,0]=temp
475
476
477
Note: See TracBrowser for help on using the repository browser.