source: anuga_work/development/2010-projects/anuga_1d/sww/sww_domain.py @ 7868

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

Added in some more c based limiters

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