source: anuga_work/development/flow_1d/sww_flow/sww_domain.py @ 7831

Last change on this file since 7831 was 7831, checked in by steve, 14 years ago
File size: 15.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 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 flow_1d.generic_1d.generic_domain import Generic_domain
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 flow_1d.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 flow_1d.utilities.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 = 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 flow_1d.sww_flow.sww_comp_flux_ext import compute_fluxes_ext_short
155
156        #self.flux_timestep = compute_fluxes_ext(timestep,self,stage,xmom,bed,height,velocity)
157
158   
159        self.flux_timestep = compute_fluxes_ext_short(timestep,self,stage,xmom,bed)
160
161
162
163    def distribute_to_vertices_and_edges(self):
164        #Call correct module function
165        #(either from this module or C-extension)
166        distribute_to_vertices_and_edges(self)
167       
168    def evolve(self, yieldstep = None, finaltime = None, duration = None,
169               skip_initial_step = False):
170        """Specialisation of basic evolve method from parent class
171        """
172
173        #Call basic machinery from parent class
174        for t in Generic_domain.evolve(self, yieldstep, finaltime,duration,
175                                       skip_initial_step):
176
177            #Pass control on to outer loop for more specific actions
178            yield(t)
179
180    def initialise_storage(self):
181        """Create and initialise self.writer object for storing data.
182        Also, save x and bed elevation
183        """
184
185        import data_manager
186
187        #Initialise writer
188        self.writer = data_manager.get_dataobject(self, mode = 'w')
189
190        #Store vertices and connectivity
191        self.writer.store_connectivity()
192
193
194    def store_timestep(self, name):
195        """Store named quantity and time.
196
197        Precondition:
198           self.write has been initialised
199        """
200        self.writer.store_timestep(name)
201
202
203#=============== End of Shallow Water Domain ===============================
204
205
206# Module functions for gradient limiting (distribute_to_vertices_and_edges)
207
208def distribute_to_vertices_and_edges(domain):
209    """Distribution from centroids to vertices specific to the
210    shallow water wave
211    equation.
212
213    It will ensure that h (w-z) is always non-negative even in the
214    presence of steep bed-slopes by taking a weighted average between shallow
215    and deep cases.
216
217    In addition, all conserved quantities get distributed as per either a
218    constant (order==1) or a piecewise linear function (order==2).
219
220    FIXME: more explanation about removal of artificial variability etc
221
222    Precondition:
223      All quantities defined at centroids and bed elevation defined at
224      vertices.
225
226    Postcondition
227      Conserved quantities defined at vertices
228
229    """
230
231    #from config import optimised_gradient_limiter
232
233    #Remove very thin layers of water
234    #protect_against_infinitesimal_and_negative_heights(domain) 
235
236    import sys
237    from flow_1d.config import epsilon, h0
238       
239    N = domain.number_of_elements
240
241    #Shortcuts
242    Stage      = domain.quantities['stage']
243    Xmom       = domain.quantities['xmomentum']
244    Bed        = domain.quantities['elevation']
245    Height     = domain.quantities['height']
246    Velocity   = domain.quantities['velocity']
247
248    #Arrays   
249    w_C   = Stage.centroid_values   
250    uh_C  = Xmom.centroid_values   
251    z_C   = Bed.centroid_values
252    h_C   = Height.centroid_values
253    u_C   = Velocity.centroid_values
254       
255    #print id(h_C)
256    #FIXME replace with numpy.where
257    for i in range(N):
258        h_C[i] = w_C[i] - z_C[i]                                               
259        if h_C[i] <= 0.0:
260            #print 'h_C[%d]= %15.5e\n' % (i,h_C[i])
261            h_C[i] = 1.0e-15
262            w_C[i] = z_C[i]
263            uh_C[i] = 0.0
264               
265   
266##    for i in range(len(h_C)):
267##      if h_C[i] < epsilon:
268##          u_C[i] = 0.0  #Could have been negative
269##          h_C[i] = 0.0
270##      else:
271
272    u_C[:,]  = uh_C/(h_C + h0/h_C)
273       
274    for name in [ 'velocity', 'stage' ]:
275        Q = domain.quantities[name]
276        if domain.order == 1:
277            Q.extrapolate_first_order()
278        elif domain.order == 2:
279            #print "add extrapolate second order to shallow water"
280            #if name != 'height':
281            Q.extrapolate_second_order()
282            #Q.limit()
283        else:
284            raise 'Unknown order'
285
286    stage_V = domain.quantities['stage'].vertex_values                 
287    bed_V   = domain.quantities['elevation'].vertex_values     
288    h_V     = domain.quantities['height'].vertex_values
289    u_V     = domain.quantities['velocity'].vertex_values
290    xmom_V  = domain.quantities['xmomentum'].vertex_values
291
292    h_V[:,:]    = stage_V - bed_V
293    for i in range(len(h_C)):
294        for j in range(2):
295            if h_V[i,j] < 0.0 :
296                #print 'h_V[%d,%d] = %f \n' % (i,j,h_V[i,j])                 
297                dh = h_V[i,j]
298                h_V[i,j] = 0.0
299                stage_V[i,j] = bed_V[i,j]
300                h_V[i,(j+1)%2] = h_V[i,(j+1)%2] + dh
301                stage_V[i,(j+1)%2] = stage_V[i,(j+1)%2] + dh
302               
303    xmom_V[:,:] = u_V * h_V
304   
305    return
306#
307
308
309
310
311
312
313
314#
315def protect_against_infinitesimal_and_negative_heights(domain):
316    """Protect against infinitesimal heights and associated high velocities
317    """
318
319    #Shortcuts
320    wc    = domain.quantities['stage'].centroid_values
321    zc    = domain.quantities['elevation'].centroid_values
322    xmomc = domain.quantities['xmomentum'].centroid_values
323    hc = wc - zc  #Water depths at centroids
324
325    zv = domain.quantities['elevation'].vertex_values
326    wv = domain.quantities['stage'].vertex_values
327    hv = wv-zv
328    xmomv = domain.quantities['xmomentum'].vertex_values
329    #remove the above two lines and corresponding code below
330
331    #Update
332    #FIXME replace with numpy.where
333    for k in range(domain.number_of_elements):
334
335        if hc[k] < domain.minimum_allowed_height:
336            #Control stage
337            if hc[k] < domain.epsilon:
338                wc[k] = zc[k] # Contain 'lost mass' error
339                wv[k,0] = zv[k,0]
340                wv[k,1] = zv[k,1]
341               
342            xmomc[k] = 0.0
343           
344        #N = domain.number_of_elements
345        #if (k == 0) | (k==N-1):
346        #    wc[k] = zc[k] # Contain 'lost mass' error
347        #    wv[k,0] = zv[k,0]
348        #    wv[k,1] = zv[k,1]
349   
350def h_limiter(domain):
351    """Limit slopes for each volume to eliminate artificial variance
352    introduced by e.g. second order extrapolator
353
354    limit on h = w-z
355
356    This limiter depends on two quantities (w,z) so it resides within
357    this module rather than within quantity.py
358    """
359
360    N = domain.number_of_elements
361    beta_h = domain.beta_h
362
363    #Shortcuts
364    wc = domain.quantities['stage'].centroid_values
365    zc = domain.quantities['elevation'].centroid_values
366    hc = wc - zc
367
368    wv = domain.quantities['stage'].vertex_values
369    zv = domain.quantities['elevation'].vertex_values
370    hv = wv-zv
371
372    hvbar = zeros(hv.shape, numpy.float) #h-limited values
373
374    #Find min and max of this and neighbour's centroid values
375    hmax = zeros(hc.shape, numpy.float)
376    hmin = zeros(hc.shape, numpy.float)
377
378    for k in range(N):
379        hmax[k] = hmin[k] = hc[k]
380        #for i in range(3):
381        for i in range(2):   
382            n = domain.neighbours[k,i]
383            if n >= 0:
384                hn = hc[n] #Neighbour's centroid value
385
386                hmin[k] = min(hmin[k], hn)
387                hmax[k] = max(hmax[k], hn)
388
389
390    #Diffences between centroids and maxima/minima
391    dhmax = hmax - hc
392    dhmin = hmin - hc
393
394    #Deltas between vertex and centroid values
395    dh = zeros(hv.shape, numpy.float)
396    #for i in range(3):
397    for i in range(2):
398        dh[:,i] = hv[:,i] - hc
399
400    #Phi limiter
401    for k in range(N):
402
403        #Find the gradient limiter (phi) across vertices
404        phi = 1.0
405        #for i in range(3):
406        for i in range(2):
407            r = 1.0
408            if (dh[k,i] > 0): r = dhmax[k]/dh[k,i]
409            if (dh[k,i] < 0): r = dhmin[k]/dh[k,i]
410
411            phi = min( min(r*beta_h, 1), phi )
412
413        #Then update using phi limiter
414        #for i in range(3):
415        for i in range(2):
416            hvbar[k,i] = hc[k] + phi*dh[k,i]
417
418    return hvbar
419
420def balance_deep_and_shallow(domain):
421    """Compute linear combination between stage as computed by
422    gradient-limiters limiting using w, and stage computed by
423    gradient-limiters limiting using h (h-limiter).
424    The former takes precedence when heights are large compared to the
425    bed slope while the latter takes precedence when heights are
426    relatively small.  Anything in between is computed as a balanced
427    linear combination in order to avoid numpyal disturbances which
428    would otherwise appear as a result of hard switching between
429    modes.
430
431    The h-limiter is always applied irrespective of the order.
432    """
433
434    #Shortcuts
435    wc = domain.quantities['stage'].centroid_values
436    zc = domain.quantities['elevation'].centroid_values
437    hc = wc - zc
438
439    wv = domain.quantities['stage'].vertex_values
440    zv = domain.quantities['elevation'].vertex_values
441    hv = wv-zv
442
443    #Limit h
444    hvbar = h_limiter(domain)
445
446    for k in range(domain.number_of_elements):
447        #Compute maximal variation in bed elevation
448        #  This quantitiy is
449        #    dz = max_i abs(z_i - z_c)
450        #  and it is independent of dimension
451        #  In the 1d case zc = (z0+z1)/2
452        #  In the 2d case zc = (z0+z1+z2)/3
453
454        dz = max(abs(zv[k,0]-zc[k]),
455                 abs(zv[k,1]-zc[k]))#,
456#                 abs(zv[k,2]-zc[k]))
457
458
459        hmin = min( hv[k,:] )
460
461        #Create alpha in [0,1], where alpha==0 means using the h-limited
462        #stage and alpha==1 means using the w-limited stage as
463        #computed by the gradient limiter (both 1st or 2nd order)
464
465        #If hmin > dz/2 then alpha = 1 and the bed will have no effect
466        #If hmin < 0 then alpha = 0 reverting to constant height above bed.
467
468        if dz > 0.0:
469            alpha = max( min( 2*hmin/dz, 1.0), 0.0 )
470        else:
471            #Flat bed
472            alpha = 1.0
473
474        alpha  = 0.0
475        #Let
476        #
477        #  wvi be the w-limited stage (wvi = zvi + hvi)
478        #  wvi- be the h-limited state (wvi- = zvi + hvi-)
479        #
480        #
481        #where i=0,1,2 denotes the vertex ids
482        #
483        #Weighted balance between w-limited and h-limited stage is
484        #
485        #  wvi := (1-alpha)*(zvi+hvi-) + alpha*(zvi+hvi)
486        #
487        #It follows that the updated wvi is
488        #  wvi := zvi + (1-alpha)*hvi- + alpha*hvi
489        #
490        # Momentum is balanced between constant and limited
491
492
493        #for i in range(3):
494        #    wv[k,i] = zv[k,i] + hvbar[k,i]
495
496        #return
497
498        if alpha < 1:
499
500            #for i in range(3):
501            for i in range(2):
502                wv[k,i] = zv[k,i] + (1.0-alpha)*hvbar[k,i] + alpha*hv[k,i]
503
504            #Momentums at centroids
505            xmomc = domain.quantities['xmomentum'].centroid_values
506#            ymomc = domain.quantities['ymomentum'].centroid_values
507
508            #Momentums at vertices
509            xmomv = domain.quantities['xmomentum'].vertex_values
510#           ymomv = domain.quantities['ymomentum'].vertex_values
511
512            # Update momentum as a linear combination of
513            # xmomc and ymomc (shallow) and momentum
514            # from extrapolator xmomv and ymomv (deep).
515            xmomv[k,:] = (1.0-alpha)*xmomc[k] + alpha*xmomv[k,:]
516#            ymomv[k,:] = (1-alpha)*ymomc[k] + alpha*ymomv[k,:]
517
518
519
Note: See TracBrowser for help on using the repository browser.