source: trunk/anuga_work/development/sudi/sw_1d/shock_detector_shv/sww_domain_shm.py @ 7906

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

Adding sudi's code

File size: 18.6 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
40John Jakeman, Ole Nielsen, Stephen Roberts, Duncan Gray, Christopher Zoppou
41Geoscience Australia, 2006
42
43Sudi Mungkasi, ANU, 2010
44"""
45
46
47from domain import *
48Generic_Domain = Domain #Rename
49
50#Shallow water domain
51class Domain(Generic_Domain):
52
53    def __init__(self, coordinates, boundary = None, tagged_elements = None):
54        conserved_quantities = ['stage', 'xmomentum'] #['height', 'xmomentum']
55        evolved_quantities = ['stage', 'xmomentum', 'elevation', 'height', 'velocity']
56        other_quantities = ['friction']
57        Generic_Domain.__init__(self, coordinates, boundary,
58                                conserved_quantities, evolved_quantities, other_quantities,
59                                tagged_elements)
60       
61        from config import minimum_allowed_height, g, h0
62        self.minimum_allowed_height = minimum_allowed_height
63        self.g = g
64        self.h0 = h0
65
66        self.forcing_terms.append(gravity_F2)
67        #self.forcing_terms.append(manning_friction)
68
69        #Realtime visualisation
70        self.visualiser = None
71        self.visualise  = False
72        self.visualise_color_stage = False
73        self.visualise_stage_range = 1.0
74        self.visualise_timer = True
75        self.visualise_range_z = None
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__ = 'sww_domain_shm'
93        self.check_integrity()
94
95
96    def set_quantities_to_be_stored(self, q):
97        """Specify which quantities will be stored in the sww file.
98
99        q must be either:
100          - the name of a quantity
101          - a list of quantity names
102          - None
103
104        In the two first cases, the named quantities will be stored at each
105        yieldstep
106        (This is in addition to the quantities elevation and friction) 
107        If q is None, storage will be switched off altogether.
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    def initialise_visualiser(self,scale_z=1.0,rect=None):
127        #Realtime visualisation
128        if self.visualiser is None:
129            from realtime_visualisation_new import Visualiser
130            self.visualiser = Visualiser(self,scale_z,rect)
131            self.visualiser.setup['elevation']=True
132            self.visualiser.updating['stage']=True
133        self.visualise = True
134        if self.visualise_color_stage == True:
135            self.visualiser.coloring['stage'] = True
136            self.visualiser.qcolor['stage'] = (0.0, 0.0, 0.8)
137        print 'initialise visualiser'
138        print self.visualiser.setup
139        print self.visualiser.updating
140
141    def check_integrity(self):
142        Generic_Domain.check_integrity(self)
143        #Check that we are solving the shallow water wave equation
144
145        msg = 'First conserved quantity must be "stage"'
146        assert self.conserved_quantities[0] == 'stage', msg
147        msg = 'Second conserved quantity must be "xmomentum"'
148        assert self.conserved_quantities[1] == 'xmomentum', msg
149               
150        msg = 'First evolved quantity must be "stage"'
151        assert self.evolved_quantities[0] == 'stage', msg
152        msg = 'Second evolved quantity must be "xmomentum"'
153        assert self.evolved_quantities[1] == 'xmomentum', msg
154        msg = 'Third evolved quantity must be "elevation"'
155        assert self.evolved_quantities[2] == 'elevation', msg
156        msg = 'Fourth evolved quantity must be "height"'
157        assert self.evolved_quantities[3] == 'height', msg
158        msg = 'Fifth evolved quantity must be "velocity"'
159        assert self.evolved_quantities[4] == 'velocity', msg
160
161    def extrapolate_second_order_sw(self):
162        #Call correct module function
163        #(either from this module or C-extension)
164        extrapolate_second_order_sw(self)
165
166    def compute_fluxes(self):
167        #Call correct module function(either from this module or C-extension)
168        compute_fluxes_C_wellbalanced(self)
169        #compute_fluxes_C_nonwellbalanced2(self)
170
171    def compute_timestep(self):
172        #Call correct module function
173        compute_timestep(self)
174
175    def distribute_to_vertices_and_edges(self):
176        #Call correct module function(either from this module or C-extension)
177        #distribute_to_vertices_and_edges_shv(self)
178        distribute_to_vertices_and_edges_shm(self)
179       
180    def evolve(self, yieldstep = None, finaltime = None, duration = None, skip_initial_step = False):
181        #Call basic machinery from parent class
182        for t in Generic_Domain.evolve(self, yieldstep, finaltime, duration, skip_initial_step):
183            yield(t)
184
185    def initialise_storage(self):
186        """Create and initialise self.writer object for storing data.
187        Also, save x and bed elevation
188        """
189
190        import data_manager
191
192        #Initialise writer
193        self.writer = data_manager.get_dataobject(self, mode = 'w')
194
195        #Store vertices and connectivity
196        self.writer.store_connectivity()
197
198
199    def store_timestep(self, name):
200        """Store named quantity and time.
201
202        Precondition:
203           self.write has been initialised
204        """
205        self.writer.store_timestep(name)
206
207
208#=============== End of Shallow Water Domain ===============================
209
210# Compute flux definition
211# ###################################
212def compute_fluxes_C_wellbalanced(domain):
213    import sys
214    from Numeric import Float
215    from numpy import zeros   
216   
217    N = domain.number_of_elements
218    timestep = float(sys.maxint)
219    epsilon = domain.epsilon
220    g = domain.g
221    neighbours = domain.neighbours
222    neighbour_vertices = domain.neighbour_vertices
223    normals = domain.normals
224    areas = domain.areas
225   
226    Stage    = domain.quantities['stage']
227    Xmom     = domain.quantities['xmomentum']
228    Bed      = domain.quantities['elevation']
229    Height   = domain.quantities['height']
230    Velocity = domain.quantities['velocity']
231
232
233    stage_boundary_values = Stage.boundary_values
234    xmom_boundary_values  = Xmom.boundary_values
235    bed_boundary_values   = Bed.boundary_values
236    height_boundary_values= Height.boundary_values
237    vel_boundary_values   = Velocity.boundary_values
238   
239    stage_explicit_update = Stage.explicit_update
240    xmom_explicit_update  = Xmom.explicit_update
241    bed_explicit_values   = Bed.explicit_update
242    height_explicit_values= Height.explicit_update
243    vel_explicit_values   = Velocity.explicit_update
244   
245    max_speed_array = domain.max_speed_array
246   
247    domain.distribute_to_vertices_and_edges()
248    domain.update_boundary()
249   
250    stage_V = Stage.vertex_values
251    xmom_V  = Xmom.vertex_values
252    bed_V   = Bed.vertex_values
253    height_V= Height.vertex_values
254    vel_V   = Velocity.vertex_values
255
256    number_of_elements = len(stage_V)
257
258    from comp_flux_ext_wellbalanced import compute_fluxes_ext_wellbalanced
259    domain.flux_timestep = compute_fluxes_ext_wellbalanced(timestep,
260                                  epsilon,
261                                  g,
262                                                           
263                                  neighbours,
264                                  neighbour_vertices,
265                                  normals,
266                                  areas,
267                                                           
268                                  stage_V,
269                                  xmom_V,
270                                  bed_V,
271                                  height_V,
272                                  vel_V,
273                                                           
274                                  stage_boundary_values,
275                                  xmom_boundary_values,
276                                  bed_boundary_values,
277                                  height_boundary_values,
278                                  vel_boundary_values,
279                                                           
280                                  stage_explicit_update,
281                                  xmom_explicit_update,
282                                  bed_explicit_values,
283                                  height_explicit_values,
284                                  vel_explicit_values,
285                                                           
286                                  number_of_elements,
287                                  max_speed_array)
288
289
290# ###################################
291# Module functions for gradient limiting (distribute_to_vertices_and_edges)
292def distribute_to_vertices_and_edges_shv(domain):
293    """Distribution from centroids to vertices specific to the
294    shallow water wave
295    equation.
296
297    It will ensure that h (w-z) is always non-negative even in the
298    presence of steep bed-slopes by taking a weighted average between shallow
299    and deep cases.
300
301    In addition, all conserved quantities get distributed as per either a
302    constant (order==1) or a piecewise linear function (order==2).
303
304    FIXME: more explanation about removal of artificial variability etc
305
306    Precondition:
307      All quantities defined at centroids and bed elevation defined at
308      vertices.
309
310    Postcondition
311      Conserved quantities defined at vertices
312
313    """
314
315    #from config import optimised_gradient_limiter
316
317    #Remove very thin layers of water
318    #protect_against_infinitesimal_and_negative_heights(domain) 
319
320    import sys
321    from Numeric import Float
322    from numpy import zeros   
323    from config import epsilon, h0
324
325    N = domain.number_of_elements
326
327    #Shortcuts
328    Stage = domain.quantities['stage']
329    Xmom = domain.quantities['xmomentum']
330    Bed = domain.quantities['elevation']
331    Height = domain.quantities['height']
332    Velocity = domain.quantities['velocity']
333
334    #Arrays   
335    w_C   = Stage.centroid_values   
336    uh_C  = Xmom.centroid_values   
337    z_C   = Bed.centroid_values
338    h_C   = Height.centroid_values
339    u_C   = Velocity.centroid_values
340
341    for i in range(N):
342        h_C[i] = w_C[i] - z_C[i]                                               
343        if h_C[i] <= epsilon:
344            uh_C[i] = 0.0
345            u_C[i] = 0.0
346            #w_C[i] = z_C[i]
347        else:
348            u_C[i]  = uh_C[i]/(h_C[i] + h0/h_C[i])
349
350    for name in ['stage', 'height', 'velocity']:
351        Q = domain.quantities[name]
352        if domain.order == 1:
353            Q.extrapolate_first_order()
354        elif domain.order == 2:
355            Q.extrapolate_second_order()
356        else:
357            raise 'Unknown order'
358
359    stage_V = domain.quantities['stage'].vertex_values
360    bed_V   = domain.quantities['elevation'].vertex_values
361    h_V     = domain.quantities['height'].vertex_values
362    u_V     = domain.quantities['velocity'].vertex_values               
363    xmom_V  = domain.quantities['xmomentum'].vertex_values
364
365    bed_V[:] = stage_V - h_V
366    xmom_V[:] = u_V * h_V
367   
368    return
369
370def distribute_to_vertices_and_edges_shm(domain):
371    # shm stands for STAGE, HEIGHT, MOMENTUM
372    """Distribution from centroids to vertices specific to the
373    shallow water wave
374    equation.
375
376    It will ensure that h (w-z) is always non-negative even in the
377    presence of steep bed-slopes by taking a weighted average between shallow
378    and deep cases.
379
380    In addition, all conserved quantities get distributed as per either a
381    constant (order==1) or a piecewise linear function (order==2).
382
383    FIXME: more explanation about removal of artificial variability etc
384
385    Precondition:
386      All quantities defined at centroids and bed elevation defined at
387      vertices.
388
389    Postcondition
390      Conserved quantities defined at vertices
391
392    """
393
394    #from config import optimised_gradient_limiter
395
396    #Remove very thin layers of water
397    #protect_against_infinitesimal_and_negative_heights(domain) 
398
399    import sys
400    from Numeric import Float
401    from numpy import array, zeros
402    from config import epsilon, h0
403
404    N = domain.number_of_elements
405
406    #Shortcuts
407    Stage = domain.quantities['stage']
408    Xmom = domain.quantities['xmomentum']
409    Bed = domain.quantities['elevation']
410    Height = domain.quantities['height']
411    Velocity = domain.quantities['velocity']
412
413    #Arrays   
414    w_C   = Stage.centroid_values   
415    uh_C  = Xmom.centroid_values   
416    z_C   = Bed.centroid_values
417    h_C   = Height.centroid_values
418    u_C   = Velocity.centroid_values
419
420   
421    for i in range(N):
422        h_C[i] = w_C[i] - z_C[i]                                               
423        if h_C[i] <= epsilon:
424            uh_C[i] = 0.0
425            u_C[i] = 0.0
426        else:
427            u_C[i] = uh_C[i]/(h_C[i] + h0/h_C[i])
428
429    for name in ['stage', 'height', 'xmomentum']:
430        Q = domain.quantities[name]
431        if domain.order == 1:
432            Q.extrapolate_first_order()
433        elif domain.order == 2:
434            Q.extrapolate_second_order()
435        else:
436            raise 'Unknown order'
437
438    stage_V = domain.quantities['stage'].vertex_values
439    bed_V   = domain.quantities['elevation'].vertex_values
440    h_V     = domain.quantities['height'].vertex_values
441    u_V     = domain.quantities['velocity'].vertex_values               
442    xmom_V  = domain.quantities['xmomentum'].vertex_values
443
444    bed_V[:] = stage_V - h_V
445
446    for i in range(N):                                         
447        if min(h_V[i]) <= 0.0:
448            h_V[i] = array([0.0, 0.0])
449            stage_V[i] = bed_V[i]
450            xmom_V[i] = array([0.0, 0.0])
451   
452    u_V[:]  = xmom_V/(h_V + h0/h_V)   
453   
454    return
455
456
457   
458#
459def protect_against_infinitesimal_and_negative_heights(domain):
460    """Protect against infinitesimal heights and associated high velocities
461    """
462
463    #Shortcuts
464    wc = domain.quantities['stage'].centroid_values
465    zc = domain.quantities['elevation'].centroid_values
466    xmomc = domain.quantities['xmomentum'].centroid_values
467    hc = wc - zc  #Water depths at centroids
468
469    zv = domain.quantities['elevation'].vertex_values
470    wv = domain.quantities['stage'].vertex_values
471    hv = wv-zv
472    xmomv = domain.quantities['xmomentum'].vertex_values
473    #remove the above two lines and corresponding code below
474
475    #Update
476    for k in range(domain.number_of_elements):
477
478        if hc[k] < domain.minimum_allowed_height:
479            #Control stage
480            if hc[k] < domain.epsilon:
481                wc[k] = zc[k] # Contain 'lost mass' error
482                wv[k,0] = zv[k,0]
483                wv[k,1] = zv[k,1]
484               
485            xmomc[k] = 0.0
486           
487   
488
489
490
491#########################
492#Standard forcing terms:
493#
494def gravity(domain):
495    """Apply gravitational pull in the presence of bed slope
496    """
497
498    from util import gradient
499    from Numeric import Float
500    from numpy import zeros, array, sum   
501
502    xmom = domain.quantities['xmomentum'].explicit_update
503    stage = domain.quantities['stage'].explicit_update
504
505    Stage = domain.quantities['stage']
506    Elevation = domain.quantities['elevation']
507    h = Stage.vertex_values - Elevation.vertex_values
508    b = Elevation.vertex_values
509    w = Stage.vertex_values
510
511    x = domain.get_vertex_coordinates()
512    g = domain.g
513
514    for k in range(domain.number_of_elements):
515        avg_h = sum( h[k,:] )/2
516
517        #Compute bed slope
518        x0, x1 = x[k,:]
519        b0, b1 = b[k,:]
520
521        w0, w1 = w[k,:]
522        wx = gradient(x0, x1, w0, w1)
523        bx = gradient(x0, x1, b0, b1)
524       
525        #Update momentum (explicit update is reset to source values)
526        xmom[k] += -g*bx*avg_h
527
528def gravity_F2(domain):
529    """Apply gravitational pull in the presence of bed slope
530    """
531
532    from util import gradient
533    from Numeric import Float
534    from numpy import zeros, array, sum   
535    from parameters import F2#This is an additional friction!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
536
537    xmom = domain.quantities['xmomentum'].explicit_update
538    stage = domain.quantities['stage'].explicit_update
539
540    Stage = domain.quantities['stage']
541    Elevation = domain.quantities['elevation']
542    h = Stage.vertex_values - Elevation.vertex_values
543    b = Elevation.vertex_values
544    w = Stage.vertex_values
545
546    x = domain.get_vertex_coordinates()
547    g = domain.g
548
549    for k in range(domain.number_of_elements):
550        avg_h = sum( h[k,:] )/2
551
552        #Compute bed slope
553        x0, x1 = x[k,:]
554        b0, b1 = b[k,:]
555
556        w0, w1 = w[k,:]
557        wx = gradient(x0, x1, w0, w1)
558        bx = gradient(x0, x1, b0, b1)
559       
560        #Update momentum (explicit update is reset to source values)
561        xmom[k] += -g*bx*avg_h + avg_h*F2#This is an additional friction!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
562
563 
564def manning_friction(domain):
565    """Apply (Manning) friction to water momentum
566    """
567
568    from math import sqrt
569
570    w = domain.quantities['stage'].centroid_values
571    z = domain.quantities['elevation'].centroid_values
572    h = w-z
573
574    uh = domain.quantities['xmomentum'].centroid_values
575    eta = domain.quantities['friction'].centroid_values
576
577    xmom_update = domain.quantities['xmomentum'].semi_implicit_update
578
579    N = domain.number_of_elements
580    eps = domain.minimum_allowed_height
581    g = domain.g
582
583    for k in range(N):
584        if eta[k] >= eps:
585            if h[k] >= eps:
586                S = -g * eta[k]**2 * uh[k]
587                S /= h[k]**(7.0/3)
588
589                #Update momentum
590                xmom_update[k] += S*uh[k]
Note: See TracBrowser for help on using the repository browser.