source: anuga_work/development/pipeflow/sww_domain.py @ 7781

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

cleaning up shallow water domain

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 Domain from module domain.py
6consisting of methods specific to the Shallow Water Wave Equation
7
8
9U_t + E_x = S
10
11where
12
13U = [w, uh]
14E = [uh, u^2h + gh^2/2]
15S represents source terms forcing the system
16(e.g. gravity, friction, wind stress, ...)
17
18and _t, _x, _y denote the derivative with respect to t, x and y respectiely.
19
20The quantities are
21
22symbol    variable name    explanation
23x         x                horizontal distance from origin [m]
24z         elevation        elevation of bed on which flow is modelled [m]
25h         height           water height above z [m]
26w         stage            absolute water level, w = z+h [m]
27u                          speed in the x direction [m/s]
28uh        xmomentum        momentum in the x direction [m^2/s]
29
30eta                        mannings friction coefficient [to appear]
31nu                         wind stress coefficient [to appear]
32
33The conserved quantities are w, uh
34
35For details see e.g.
36Christopher Zoppou and Stephen Roberts,
37Catastrophic Collapse of Water Supply Reservoirs in Urban Areas,
38Journal of Hydraulic Engineering, vol. 127, No. 7 July 1999
39
40
41John Jakeman, Ole Nielsen, Stephen Roberts, Duncan Gray, Christopher Zoppou
42Geoscience Australia, 2006
43"""
44
45import numpy
46
47from generic_domain import *
48from sww_boundary_conditions import *
49from sww_forcing_terms import *
50
51#Shallow water domain
52class Domain(Generic_Domain):
53
54    def __init__(self, coordinates, boundary = None, tagged_elements = None):
55
56        conserved_quantities = ['stage', 'xmomentum']
57        evolved_quantities   = ['stage', 'xmomentum', 'elevation', 'height', 'velocity']
58        other_quantities     = ['friction']
59        Generic_Domain.__init__(self,
60                                coordinates          = coordinates,
61                                boundary             = boundary,
62                                conserved_quantities = conserved_quantities,
63                                evolved_quantities   = evolved_quantities,
64                                other_quantities     = other_quantities,
65                                tagged_elements      = tagged_elements)
66       
67        from config import minimum_allowed_height, g, h0
68        self.minimum_allowed_height = minimum_allowed_height
69        self.g = g
70        self.h0 = h0
71
72        #forcing terms
73        self.forcing_terms.append(gravity)
74        #self.forcing_terms.append(manning_friction)
75       
76       
77        #Stored output
78        self.store = True
79        self.format = 'sww'
80        self.smooth = True
81
82        #Evolve parametrs
83        self.cfl = 1.0
84       
85        #Reduction operation for get_vertex_values
86        from util import mean
87        self.reduction = mean
88        #self.reduction = min  #Looks better near steep slopes
89
90        self.quantities_to_be_stored = ['stage','xmomentum']
91
92        self.__doc__ = 'shallow_water_domain'
93
94
95    def set_quantities_to_be_stored(self, q):
96        """Specify which quantities will be stored in the sww file.
97
98        q must be either:
99          - the name of a quantity
100          - a list of quantity names
101          - None
102
103        In the two first cases, the named quantities will be stored at each
104        yieldstep
105        (This is in addition to the quantities elevation and friction) 
106        If q is None, storage will be switched off altogether.
107        """
108
109
110        if q is None:
111            self.quantities_to_be_stored = []           
112            self.store = False
113            return
114
115        if isinstance(q, basestring):
116            q = [q] # Turn argument into a list
117
118        #Check correcness   
119        for quantity_name in q:
120            msg = 'Quantity %s is not a valid conserved quantity' %quantity_name
121            assert quantity_name in self.conserved_quantities, msg
122       
123        self.quantities_to_be_stored = q
124       
125
126
127    def check_integrity(self):
128        Generic_Domain.check_integrity(self)
129        #Check that we are solving the shallow water wave equation
130
131        msg = 'First conserved quantity must be "stage"'
132        assert self.conserved_quantities[0] == 'stage', msg
133        msg = 'Second conserved quantity must be "xmomentum"'
134        assert self.conserved_quantities[1] == 'xmomentum', msg
135
136
137
138    def compute_fluxes(self):
139        #Call correct module function
140        #(either from this module or C-extension)
141
142        import sys
143
144
145        timestep = numpy.float(sys.maxint)
146
147        stage    = self.quantities['stage']
148        xmom     = self.quantities['xmomentum']
149        bed      = self.quantities['elevation']
150        height   = self.quantities['height']
151        velocity = self.quantities['velocity']
152
153
154        from sww_comp_flux_ext import compute_fluxes_ext
155
156        self.flux_timestep = compute_fluxes_ext(timestep,self,stage,xmom,bed,height,velocity)
157
158
159    def distribute_to_vertices_and_edges(self):
160        #Call correct module function
161        #(either from this module or C-extension)
162        distribute_to_vertices_and_edges(self)
163       
164    def evolve(self, yieldstep = None, finaltime = None, duration = None,
165               skip_initial_step = False):
166        """Specialisation of basic evolve method from parent class
167        """
168
169        #Call basic machinery from parent class
170        for t in Generic_Domain.evolve(self, yieldstep, finaltime,duration,
171                                       skip_initial_step):
172
173            #Pass control on to outer loop for more specific actions
174            yield(t)
175
176    def initialise_storage(self):
177        """Create and initialise self.writer object for storing data.
178        Also, save x and bed elevation
179        """
180
181        import data_manager
182
183        #Initialise writer
184        self.writer = data_manager.get_dataobject(self, mode = 'w')
185
186        #Store vertices and connectivity
187        self.writer.store_connectivity()
188
189
190    def store_timestep(self, name):
191        """Store named quantity and time.
192
193        Precondition:
194           self.write has been initialised
195        """
196        self.writer.store_timestep(name)
197
198
199#=============== End of Shallow Water Domain ===============================
200
201
202# Module functions for gradient limiting (distribute_to_vertices_and_edges)
203
204def distribute_to_vertices_and_edges(domain):
205    """Distribution from centroids to vertices specific to the
206    shallow water wave
207    equation.
208
209    It will ensure that h (w-z) is always non-negative even in the
210    presence of steep bed-slopes by taking a weighted average between shallow
211    and deep cases.
212
213    In addition, all conserved quantities get distributed as per either a
214    constant (order==1) or a piecewise linear function (order==2).
215
216    FIXME: more explanation about removal of artificial variability etc
217
218    Precondition:
219      All quantities defined at centroids and bed elevation defined at
220      vertices.
221
222    Postcondition
223      Conserved quantities defined at vertices
224
225    """
226
227    #from config import optimised_gradient_limiter
228
229    #Remove very thin layers of water
230    #protect_against_infinitesimal_and_negative_heights(domain) 
231
232    import sys
233    from config import epsilon, h0
234       
235    N = domain.number_of_elements
236
237    #Shortcuts
238    Stage      = domain.quantities['stage']
239    Xmom       = domain.quantities['xmomentum']
240    Bed        = domain.quantities['elevation']
241    Height     = domain.quantities['height']
242    Velocity   = domain.quantities['velocity']
243
244    #Arrays   
245    w_C   = Stage.centroid_values   
246    uh_C  = Xmom.centroid_values   
247    z_C   = Bed.centroid_values
248    h_C   = Height.centroid_values
249    u_C   = Velocity.centroid_values
250       
251    #print id(h_C)
252    #FIXME replace with numpy.where
253    for i in range(N):
254        h_C[i] = w_C[i] - z_C[i]                                               
255        if h_C[i] <= 0.0:
256            #print 'h_C[%d]= %15.5e\n' % (i,h_C[i])
257            h_C[i] = 1.0e-15
258            w_C[i] = z_C[i]
259            uh_C[i] = 0.0
260               
261   
262##    for i in range(len(h_C)):
263##      if h_C[i] < epsilon:
264##          u_C[i] = 0.0  #Could have been negative
265##          h_C[i] = 0.0
266##      else:
267
268    u_C[:]  = uh_C/(h_C + h0/h_C)
269       
270    for name in [ 'velocity', 'stage' ]:
271        Q = domain.quantities[name]
272        if domain.order == 1:
273            Q.extrapolate_first_order()
274        elif domain.order == 2:
275            #print "add extrapolate second order to shallow water"
276            #if name != 'height':
277            Q.extrapolate_second_order()
278            #Q.limit()
279        else:
280            raise 'Unknown order'
281
282    stage_V = domain.quantities['stage'].vertex_values                 
283    bed_V   = domain.quantities['elevation'].vertex_values     
284    h_V     = domain.quantities['height'].vertex_values
285    u_V     = domain.quantities['velocity'].vertex_values
286    xmom_V  = domain.quantities['xmomentum'].vertex_values
287
288    h_V[:]    = stage_V - bed_V
289    for i in range(len(h_C)):
290        for j in range(2):
291            if h_V[i,j] < 0.0 :
292                #print 'h_V[%d,%d] = %f \n' % (i,j,h_V[i,j])                 
293                dh = h_V[i,j]
294                h_V[i,j] = 0.0
295                stage_V[i,j] = bed_V[i,j]
296                h_V[i,(j+1)%2] = h_V[i,(j+1)%2] + dh
297                stage_V[i,(j+1)%2] = stage_V[i,(j+1)%2] + dh
298               
299    xmom_V[:] = u_V * h_V
300   
301    return
302#
303
304
305
306
307
308
309
310#
311def protect_against_infinitesimal_and_negative_heights(domain):
312    """Protect against infinitesimal heights and associated high velocities
313    """
314
315    #Shortcuts
316    wc = domain.quantities['stage'].centroid_values
317    zc = domain.quantities['elevation'].centroid_values
318    xmomc = domain.quantities['xmomentum'].centroid_values
319    hc = wc - zc  #Water depths at centroids
320
321    zv = domain.quantities['elevation'].vertex_values
322    wv = domain.quantities['stage'].vertex_values
323    hv = wv-zv
324    xmomv = domain.quantities['xmomentum'].vertex_values
325    #remove the above two lines and corresponding code below
326
327    #Update
328    #FIXME replace with numpy.where
329    for k in range(domain.number_of_elements):
330
331        if hc[k] < domain.minimum_allowed_height:
332            #Control stage
333            if hc[k] < domain.epsilon:
334                wc[k] = zc[k] # Contain 'lost mass' error
335                wv[k,0] = zv[k,0]
336                wv[k,1] = zv[k,1]
337               
338            xmomc[k] = 0.0
339           
340        #N = domain.number_of_elements
341        #if (k == 0) | (k==N-1):
342        #    wc[k] = zc[k] # Contain 'lost mass' error
343        #    wv[k,0] = zv[k,0]
344        #    wv[k,1] = zv[k,1]
345   
346def h_limiter(domain):
347    """Limit slopes for each volume to eliminate artificial variance
348    introduced by e.g. second order extrapolator
349
350    limit on h = w-z
351
352    This limiter depends on two quantities (w,z) so it resides within
353    this module rather than within quantity.py
354    """
355
356    N = domain.number_of_elements
357    beta_h = domain.beta_h
358
359    #Shortcuts
360    wc = domain.quantities['stage'].centroid_values
361    zc = domain.quantities['elevation'].centroid_values
362    hc = wc - zc
363
364    wv = domain.quantities['stage'].vertex_values
365    zv = domain.quantities['elevation'].vertex_values
366    hv = wv-zv
367
368    hvbar = zeros(hv.shape, numpy.float) #h-limited values
369
370    #Find min and max of this and neighbour's centroid values
371    hmax = zeros(hc.shape, numpy.float)
372    hmin = zeros(hc.shape, numpy.float)
373
374    for k in range(N):
375        hmax[k] = hmin[k] = hc[k]
376        #for i in range(3):
377        for i in range(2):   
378            n = domain.neighbours[k,i]
379            if n >= 0:
380                hn = hc[n] #Neighbour's centroid value
381
382                hmin[k] = min(hmin[k], hn)
383                hmax[k] = max(hmax[k], hn)
384
385
386    #Diffences between centroids and maxima/minima
387    dhmax = hmax - hc
388    dhmin = hmin - hc
389
390    #Deltas between vertex and centroid values
391    dh = zeros(hv.shape, numpy.float)
392    #for i in range(3):
393    for i in range(2):
394        dh[:,i] = hv[:,i] - hc
395
396    #Phi limiter
397    for k in range(N):
398
399        #Find the gradient limiter (phi) across vertices
400        phi = 1.0
401        #for i in range(3):
402        for i in range(2):
403            r = 1.0
404            if (dh[k,i] > 0): r = dhmax[k]/dh[k,i]
405            if (dh[k,i] < 0): r = dhmin[k]/dh[k,i]
406
407            phi = min( min(r*beta_h, 1), phi )
408
409        #Then update using phi limiter
410        #for i in range(3):
411        for i in range(2):
412            hvbar[k,i] = hc[k] + phi*dh[k,i]
413
414    return hvbar
415
416def balance_deep_and_shallow(domain):
417    """Compute linear combination between stage as computed by
418    gradient-limiters limiting using w, and stage computed by
419    gradient-limiters limiting using h (h-limiter).
420    The former takes precedence when heights are large compared to the
421    bed slope while the latter takes precedence when heights are
422    relatively small.  Anything in between is computed as a balanced
423    linear combination in order to avoid numpyal disturbances which
424    would otherwise appear as a result of hard switching between
425    modes.
426
427    The h-limiter is always applied irrespective of the order.
428    """
429
430    #Shortcuts
431    wc = domain.quantities['stage'].centroid_values
432    zc = domain.quantities['elevation'].centroid_values
433    hc = wc - zc
434
435    wv = domain.quantities['stage'].vertex_values
436    zv = domain.quantities['elevation'].vertex_values
437    hv = wv-zv
438
439    #Limit h
440    hvbar = h_limiter(domain)
441
442    for k in range(domain.number_of_elements):
443        #Compute maximal variation in bed elevation
444        #  This quantitiy is
445        #    dz = max_i abs(z_i - z_c)
446        #  and it is independent of dimension
447        #  In the 1d case zc = (z0+z1)/2
448        #  In the 2d case zc = (z0+z1+z2)/3
449
450        dz = max(abs(zv[k,0]-zc[k]),
451                 abs(zv[k,1]-zc[k]))#,
452#                 abs(zv[k,2]-zc[k]))
453
454
455        hmin = min( hv[k,:] )
456
457        #Create alpha in [0,1], where alpha==0 means using the h-limited
458        #stage and alpha==1 means using the w-limited stage as
459        #computed by the gradient limiter (both 1st or 2nd order)
460
461        #If hmin > dz/2 then alpha = 1 and the bed will have no effect
462        #If hmin < 0 then alpha = 0 reverting to constant height above bed.
463
464        if dz > 0.0:
465            alpha = max( min( 2*hmin/dz, 1.0), 0.0 )
466        else:
467            #Flat bed
468            alpha = 1.0
469
470        alpha  = 0.0
471        #Let
472        #
473        #  wvi be the w-limited stage (wvi = zvi + hvi)
474        #  wvi- be the h-limited state (wvi- = zvi + hvi-)
475        #
476        #
477        #where i=0,1,2 denotes the vertex ids
478        #
479        #Weighted balance between w-limited and h-limited stage is
480        #
481        #  wvi := (1-alpha)*(zvi+hvi-) + alpha*(zvi+hvi)
482        #
483        #It follows that the updated wvi is
484        #  wvi := zvi + (1-alpha)*hvi- + alpha*hvi
485        #
486        # Momentum is balanced between constant and limited
487
488
489        #for i in range(3):
490        #    wv[k,i] = zv[k,i] + hvbar[k,i]
491
492        #return
493
494        if alpha < 1:
495
496            #for i in range(3):
497            for i in range(2):
498                wv[k,i] = zv[k,i] + (1.0-alpha)*hvbar[k,i] + alpha*hv[k,i]
499
500            #Momentums at centroids
501            xmomc = domain.quantities['xmomentum'].centroid_values
502#            ymomc = domain.quantities['ymomentum'].centroid_values
503
504            #Momentums at vertices
505            xmomv = domain.quantities['xmomentum'].vertex_values
506#           ymomv = domain.quantities['ymomentum'].vertex_values
507
508            # Update momentum as a linear combination of
509            # xmomc and ymomc (shallow) and momentum
510            # from extrapolator xmomv and ymomv (deep).
511            xmomv[k,:] = (1.0-alpha)*xmomc[k] + alpha*xmomv[k,:]
512#            ymomv[k,:] = (1-alpha)*ymomc[k] + alpha*ymomv[k,:]
513
514
515
Note: See TracBrowser for help on using the repository browser.