source: anuga_work/development/shallow_water_1d/shallow_water_1d.py @ 4960

Last change on this file since 4960 was 4960, checked in by steve, 16 years ago
File size: 39.7 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
45#from domain import *
46#from domain_order2 import *
47from domain_t2 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
55        conserved_quantities = ['stage', 'xmomentum']
56        other_quantities = ['elevation', 'friction']
57        Generic_Domain.__init__(self, coordinates, boundary,
58                                conserved_quantities, other_quantities,
59                                tagged_elements)
60       
61        from config import minimum_allowed_height, g
62        self.minimum_allowed_height = minimum_allowed_height
63        self.g = g
64
65        #forcing terms not included in 1d domain ?WHy?
66        self.forcing_terms.append(gravity)
67        #self.forcing_terms.append(manning_friction)
68        #print "\nI have Removed forcing terms line 64 1dsw"
69
70        #Realtime visualisation
71        self.visualiser = None
72        self.visualise  = False
73        self.visualise_color_stage = False
74        self.visualise_stage_range = 1.0
75        self.visualise_timer = True
76        self.visualise_range_z = None
77       
78        #Stored output
79        self.store = True
80        self.format = 'sww'
81        self.smooth = True
82
83        #Evolve parametrs
84        self.cfl = 1.0
85       
86        #Reduction operation for get_vertex_values
87        from util import mean
88        self.reduction = mean
89        #self.reduction = min  #Looks better near steep slopes
90
91        self.quantities_to_be_stored = ['stage','xmomentum']
92
93
94    def set_quantities_to_be_stored(self, q):
95        """Specify which quantities will be stored in the sww file.
96
97        q must be either:
98          - the name of a quantity
99          - a list of quantity names
100          - None
101
102        In the two first cases, the named quantities will be stored at each
103        yieldstep
104        (This is in addition to the quantities elevation and friction) 
105        If q is None, storage will be switched off altogether.
106        """
107
108
109        if q is None:
110            self.quantities_to_be_stored = []           
111            self.store = False
112            return
113
114        if isinstance(q, basestring):
115            q = [q] # Turn argument into a list
116
117        #Check correcness   
118        for quantity_name in q:
119            msg = 'Quantity %s is not a valid conserved quantity' %quantity_name
120            assert quantity_name in self.conserved_quantities, msg
121       
122        self.quantities_to_be_stored = q
123       
124
125    def initialise_visualiser(self,scale_z=1.0,rect=None):
126        #Realtime visualisation
127        if self.visualiser is None:
128            from realtime_visualisation_new import Visualiser
129            self.visualiser = Visualiser(self,scale_z,rect)
130            self.visualiser.setup['elevation']=True
131            self.visualiser.updating['stage']=True
132        self.visualise = True
133        if self.visualise_color_stage == True:
134            self.visualiser.coloring['stage'] = True
135            self.visualiser.qcolor['stage'] = (0.0, 0.0, 0.8)
136        print 'initialise visualiser'
137        print self.visualiser.setup
138        print self.visualiser.updating
139
140    def check_integrity(self):
141        Generic_Domain.check_integrity(self)
142        #Check that we are solving the shallow water wave equation
143
144        msg = 'First conserved quantity must be "stage"'
145        assert self.conserved_quantities[0] == 'stage', msg
146        msg = 'Second conserved quantity must be "xmomentum"'
147        assert self.conserved_quantities[1] == 'xmomentum', msg
148
149    def extrapolate_second_order_sw(self):
150        #Call correct module function
151        #(either from this module or C-extension)
152        extrapolate_second_order_sw(self)
153
154    def compute_fluxes(self):
155        #Call correct module function
156        #(either from this module or C-extension)
157        compute_fluxes(self)
158
159    def compute_timestep(self):
160        #Call correct module function
161        compute_timestep(self)
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,
169               skip_initial_step = False):
170        """Specialisation of basic evolve method from parent class
171        """
172
173        #Call check integrity here rather than from user scripts
174        #self.check_integrity()
175
176        #msg = 'Parameter beta_h must be in the interval [0, 1)'
177        #assert 0 <= self.beta_h < 1.0, msg
178        #msg = 'Parameter beta_w must be in the interval [0, 1)'
179        #assert 0 <= self.beta_w < 1.0, msg
180
181
182        #Initial update of vertex and edge values before any storage
183        #and or visualisation
184       
185        self.distribute_to_vertices_and_edges()
186       
187        #Initialise real time viz if requested
188        #if self.visualise is True and self.time == 0.0:
189        #    if self.visualiser is None:
190        #        self.initialise_visualiser()
191        #
192        #    self.visualiser.update_timer()
193        #    self.visualiser.setup_all()
194
195        #Store model data, e.g. for visualisation
196        #if self.store is True and self.time == 0.0:
197        #    self.initialise_storage()
198        #    #print 'Storing results in ' + self.writer.filename
199        #else:
200        #    pass
201        #    #print 'Results will not be stored.'
202        #    #print 'To store results set domain.store = True'
203        #    #FIXME: Diagnostic output should be controlled by
204        #    # a 'verbose' flag living in domain (or in a parent class)
205
206        #Call basic machinery from parent class
207        for t in Generic_Domain.evolve(self, yieldstep, finaltime,
208                                       skip_initial_step):
209            #Real time viz
210        #    if self.visualise is True:
211        #        self.visualiser.update_all()
212        #        self.visualiser.update_timer()
213
214
215            #Store model data, e.g. for subsequent visualisation
216        #    if self.store is True:
217        #        self.store_timestep(self.quantities_to_be_stored)
218
219            #FIXME: Could maybe be taken from specified list
220            #of 'store every step' quantities
221
222            #Pass control on to outer loop for more specific actions
223            yield(t)
224
225    def initialise_storage(self):
226        """Create and initialise self.writer object for storing data.
227        Also, save x and bed elevation
228        """
229
230        import data_manager
231
232        #Initialise writer
233        self.writer = data_manager.get_dataobject(self, mode = 'w')
234
235        #Store vertices and connectivity
236        self.writer.store_connectivity()
237
238
239    def store_timestep(self, name):
240        """Store named quantity and time.
241
242        Precondition:
243           self.write has been initialised
244        """
245        self.writer.store_timestep(name)
246
247
248#=============== End of Shallow Water Domain ===============================
249
250#Rotation of momentum vector
251def rotate(q, normal, direction = 1):
252    """Rotate the momentum component q (q[1], q[2])
253    from x,y coordinates to coordinates based on normal vector.
254
255    If direction is negative the rotation is inverted.
256
257    Input vector is preserved
258
259    This function is specific to the shallow water wave equation
260    """
261
262    from Numeric import zeros, Float
263
264    assert len(q) == 3,\
265           'Vector of conserved quantities must have length 3'\
266           'for 2D shallow water equation'
267
268    try:
269        l = len(normal)
270    except:
271        raise 'Normal vector must be an Numeric array'
272
273    assert l == 2, 'Normal vector must have 2 components'
274
275
276    n1 = normal[0]
277    n2 = normal[1]
278
279    r = zeros(len(q), Float) #Rotated quantities
280    r[0] = q[0]              #First quantity, height, is not rotated
281
282    if direction == -1:
283        n2 = -n2
284
285
286    r[1] =  n1*q[1] + n2*q[2]
287    r[2] = -n2*q[1] + n1*q[2]
288
289    return r
290
291
292def flux_function(normal, ql, qr, zl, zr):
293    """Compute fluxes between volumes for the shallow water wave equation
294    cast in terms of w = h+z using the 'central scheme' as described in
295
296    Kurganov, Noelle, Petrova. 'Semidiscrete Central-Upwind Schemes For
297    Hyperbolic Conservation Laws and Hamilton-Jacobi Equations'.
298    Siam J. Sci. Comput. Vol. 23, No. 3, pp. 707-740.
299
300    The implemented formula is given in equation (3.15) on page 714
301
302    Conserved quantities w, uh, are stored as elements 0 and 1
303    in the numerical vectors ql an qr.
304
305    Bed elevations zl and zr.
306    """
307
308    from config import g, epsilon
309    from math import sqrt
310    from Numeric import array
311
312    #print 'ql',ql
313
314    #Align momentums with x-axis
315    #q_left  = rotate(ql, normal, direction = 1)
316    #q_right = rotate(qr, normal, direction = 1)
317    q_left = ql
318    q_left[1] = q_left[1]*normal
319    q_right = qr
320    q_right[1] = q_right[1]*normal
321
322    #z = (zl+zr)/2 #Take average of field values
323    z = 0.5*(zl+zr) #Take average of field values
324
325    w_left  = q_left[0]   #w=h+z
326    h_left  = w_left-z
327    uh_left = q_left[1]
328
329    if h_left < epsilon:
330        u_left = 0.0  #Could have been negative
331        h_left = 0.0
332    else:
333        u_left  = uh_left/h_left
334
335
336    w_right  = q_right[0]  #w=h+z
337    h_right  = w_right-z
338    uh_right = q_right[1]
339
340
341    if h_right < epsilon:
342        u_right = 0.0  #Could have been negative
343        h_right = 0.0
344    else:
345        u_right  = uh_right/h_right
346
347    #vh_left  = q_left[2]
348    #vh_right = q_right[2]
349
350    #print h_right
351    #print u_right
352    #print h_left
353    #print u_right
354
355    soundspeed_left  = sqrt(g*h_left)
356    soundspeed_right = sqrt(g*h_right)
357
358    #Maximal wave speed
359    s_max = max(u_left+soundspeed_left, u_right+soundspeed_right, 0)
360   
361    #Minimal wave speed
362    s_min = min(u_left-soundspeed_left, u_right-soundspeed_right, 0)
363   
364    #Flux computation
365
366    #flux_left  = array([u_left*h_left,
367    #                    u_left*uh_left + 0.5*g*h_left**2])
368    #flux_right = array([u_right*h_right,
369    #                    u_right*uh_right + 0.5*g*h_right**2])
370    flux_left  = array([u_left*h_left,
371                        u_left*uh_left + 0.5*g*h_left*h_left])
372    flux_right = array([u_right*h_right,
373                        u_right*uh_right + 0.5*g*h_right*h_right])
374
375    denom = s_max-s_min
376    if denom == 0.0:
377        edgeflux = array([0.0, 0.0])
378        max_speed = 0.0
379    else:
380        edgeflux = (s_max*flux_left - s_min*flux_right)/denom
381        edgeflux += s_max*s_min*(q_right-q_left)/denom
382       
383        edgeflux[1] = edgeflux[1]*normal
384
385        max_speed = max(abs(s_max), abs(s_min))
386
387    return edgeflux, max_speed
388
389
390def flux_function_split(normal, ql, qr, zl, zr):
391    from config import g, epsilon
392    from math import sqrt
393    from Numeric import array
394
395    #print 'ql',ql
396
397    #Align momentums with x-axis
398    #q_left  = rotate(ql, normal, direction = 1)
399    #q_right = rotate(qr, normal, direction = 1)
400    q_left = ql
401    q_left[1] = q_left[1]*normal
402    q_right = qr
403    q_right[1] = q_right[1]*normal
404
405    #z = (zl+zr)/2 #Take average of field values
406    z = 0.5*(zl+zr) #Take average of field values
407
408    w_left  = q_left[0]   #w=h+z
409    h_left  = w_left-z
410    uh_left = q_left[1]
411
412    if h_left < epsilon:
413        u_left = 0.0  #Could have been negative
414        h_left = 0.0
415    else:
416        u_left  = uh_left/h_left
417
418
419    w_right  = q_right[0]  #w=h+z
420    h_right  = w_right-z
421    uh_right = q_right[1]
422
423
424    if h_right < epsilon:
425        u_right = 0.0  #Could have been negative
426        h_right = 0.0
427    else:
428        u_right  = uh_right/h_right
429
430    #vh_left  = q_left[2]
431    #vh_right = q_right[2]
432
433    #soundspeed_left  = sqrt(g*h_left)
434    #soundspeed_right = sqrt(g*h_right)
435
436    #Maximal wave speed
437    #s_max = max(u_left+soundspeed_left, u_right+soundspeed_right, 0)
438    s_max = max(u_left, u_right, 0)
439   
440    #Minimal wave speed
441    #s_min = min(u_left-soundspeed_left, u_right-soundspeed_right, 0)
442    s_min = min(u_left, u_right, 0)
443   
444    #Flux computation
445
446    #flux_left  = array([u_left*h_left,
447    #                    u_left*uh_left + 0.5*g*h_left*h_left])
448    #flux_right = array([u_right*h_right,
449    #                    u_right*uh_right + 0.5*g*h_right*h_right])
450    flux_left  = array([u_left*h_left,
451                        u_left*uh_left])# + 0.5*g*h_left*h_left])
452    flux_right = array([u_right*h_right,
453                        u_right*uh_right])# + 0.5*g*h_right*h_right])
454   
455    denom = s_max-s_min
456    if denom == 0.0:
457        edgeflux = array([0.0, 0.0])
458        max_speed = 0.0
459    else:
460        edgeflux = (s_max*flux_left - s_min*flux_right)/denom
461        edgeflux += s_max*s_min*(q_right-q_left)/denom
462       
463        edgeflux[1] = edgeflux[1]*normal
464
465        max_speed = max(abs(s_max), abs(s_min))
466
467    return edgeflux, max_speed   
468
469def compute_timestep(domain):
470    import sys
471    from Numeric import zeros, Float
472
473    N = domain.number_of_elements
474
475    #Shortcuts
476    Stage = domain.quantities['stage']
477    Xmom = domain.quantities['xmomentum']
478    Bed = domain.quantities['elevation']
479   
480    stage = Stage.vertex_values
481    xmom =  Xmom.vertex_values
482    bed =   Bed.vertex_values
483
484    stage_bdry = Stage.boundary_values
485    xmom_bdry =  Xmom.boundary_values
486
487    flux = zeros(2, Float) #Work array for summing up fluxes
488    ql = zeros(2, Float)
489    qr = zeros(2, Float)
490
491    #Loop
492    timestep = float(sys.maxint)
493    enter = True
494    for k in range(N):
495
496        flux[:] = 0.  #Reset work array
497        for i in range(2):
498            #Quantities inside volume facing neighbour i
499            ql = [stage[k, i], xmom[k, i]]
500            zl = bed[k, i]
501
502            #Quantities at neighbour on nearest face
503            n = domain.neighbours[k,i]
504            if n < 0:
505                m = -n-1 #Convert negative flag to index
506                qr[0] = stage_bdry[m]
507                qr[1] = xmom_bdry[m]
508                zr = zl #Extend bed elevation to boundary
509            else:
510                #m = domain.neighbour_edges[k,i]
511                m = domain.neighbour_vertices[k,i]
512                qr[0] = stage[n, m]
513                qr[1] =  xmom[n, m]
514                zr = bed[n, m]
515
516
517            #Outward pointing normal vector
518            normal = domain.normals[k, i]
519
520            if domain.split == False:
521                edgeflux, max_speed = flux_function(normal, ql, qr, zl, zr)
522            elif domain.split == True:
523                edgeflux, max_speed = flux_function_split(normal, ql, qr, zl, zr)
524            #Update optimal_timestep
525            try:
526                timestep = min(timestep, domain.cfl*0.5*domain.areas[k]/max_speed)
527            except ZeroDivisionError:
528                pass
529
530    domain.timestep = timestep
531
532def compute_fluxes(domain):
533    """Compute all fluxes and the timestep suitable for all volumes
534    in domain.
535
536    Compute total flux for each conserved quantity using "flux_function"
537
538    Fluxes across each edge are scaled by edgelengths and summed up
539    Resulting flux is then scaled by area and stored in
540    explicit_update for each of the three conserved quantities
541    stage, xmomentum and ymomentum
542
543    The maximal allowable speed computed by the flux_function for each volume
544    is converted to a timestep that must not be exceeded. The minimum of
545    those is computed as the next overall timestep.
546
547    Post conditions:
548      domain.explicit_update is reset to computed flux values
549      domain.timestep is set to the largest step satisfying all volumes.
550    """
551
552    import sys
553    from Numeric import zeros, Float
554
555    N = domain.number_of_elements
556
557    #Shortcuts
558    Stage = domain.quantities['stage']
559    Xmom = domain.quantities['xmomentum']
560#    Ymom = domain.quantities['ymomentum']
561    Bed = domain.quantities['elevation']
562
563    #Arrays
564    #stage = Stage.edge_values
565    #xmom =  Xmom.edge_values
566 #   ymom =  Ymom.edge_values
567    #bed =   Bed.edge_values
568   
569    stage = Stage.vertex_values
570    xmom =  Xmom.vertex_values
571    bed =   Bed.vertex_values
572   
573    #print 'stage edge values', stage
574    #print 'xmom edge values', xmom
575    #print 'bed values', bed
576
577    stage_bdry = Stage.boundary_values
578    xmom_bdry =  Xmom.boundary_values
579    #print 'stage_bdry',stage_bdry
580    #print 'xmom_bdry', xmom_bdry
581#    ymom_bdry =  Ymom.boundary_values
582
583#    flux = zeros(3, Float) #Work array for summing up fluxes
584    flux = zeros(2, Float) #Work array for summing up fluxes
585    ql = zeros(2, Float)
586    qr = zeros(2, Float)
587
588    #Loop
589    timestep = float(sys.maxint)
590    enter = True
591    for k in range(N):
592
593        flux[:] = 0.  #Reset work array
594        #for i in range(3):
595        for i in range(2):
596            #Quantities inside volume facing neighbour i
597            #ql[0] = stage[k, i]
598            #ql[1] = xmom[k, i]
599            ql = [stage[k, i], xmom[k, i]]
600            zl = bed[k, i]
601
602            #Quantities at neighbour on nearest face
603            n = domain.neighbours[k,i]
604            if n < 0:
605                m = -n-1 #Convert negative flag to index
606                qr[0] = stage_bdry[m]
607                qr[1] = xmom_bdry[m]
608                zr = zl #Extend bed elevation to boundary
609            else:
610                #m = domain.neighbour_edges[k,i]
611                m = domain.neighbour_vertices[k,i]
612                #qr = [stage[n, m], xmom[n, m], ymom[n, m]]
613                qr[0] = stage[n, m]
614                qr[1] =  xmom[n, m]
615                zr = bed[n, m]
616
617
618            #Outward pointing normal vector
619            normal = domain.normals[k, i]
620       
621            #Flux computation using provided function
622            #edgeflux, max_speed = flux_function(normal, ql, qr, zl, zr)
623            #print 'ql',ql
624            #print 'qr',qr
625           
626
627            if domain.split == False:
628                edgeflux, max_speed = flux_function(normal, ql, qr, zl, zr)
629            elif domain.split == True:
630                edgeflux, max_speed = flux_function_split(normal, ql, qr, zl, zr)
631            #print 'edgeflux', edgeflux
632
633            # THIS IS THE LINE TO DEAL WITH LEFT AND RIGHT FLUXES
634            # flux = edgefluxleft - edgefluxright
635            flux -= edgeflux #* domain.edgelengths[k,i]
636            #Update optimal_timestep
637            try:
638                #timestep = min(timestep, 0.5*domain.radii[k]/max_speed)
639                timestep = min(timestep, domain.cfl*0.5*domain.areas[k]/max_speed)
640            except ZeroDivisionError:
641                pass
642
643        #Normalise by area and store for when all conserved
644        #quantities get updated
645        flux /= domain.areas[k]
646
647        Stage.explicit_update[k] = flux[0]
648        Xmom.explicit_update[k] = flux[1]
649        #Ymom.explicit_update[k] = flux[2]
650        #print "flux cell",k,flux[0]
651
652    domain.timestep = timestep
653    #print domain.quantities['stage'].centroid_values
654
655####################################
656
657# Module functions for gradient limiting (distribute_to_vertices_and_edges)
658
659def distribute_to_vertices_and_edges(domain):
660    """Distribution from centroids to vertices specific to the
661    shallow water wave
662    equation.
663
664    It will ensure that h (w-z) is always non-negative even in the
665    presence of steep bed-slopes by taking a weighted average between shallow
666    and deep cases.
667
668    In addition, all conserved quantities get distributed as per either a
669    constant (order==1) or a piecewise linear function (order==2).
670
671    FIXME: more explanation about removal of artificial variability etc
672
673    Precondition:
674      All quantities defined at centroids and bed elevation defined at
675      vertices.
676
677    Postcondition
678      Conserved quantities defined at vertices
679
680    """
681
682    #from config import optimised_gradient_limiter
683
684    #Remove very thin layers of water
685    protect_against_infinitesimal_and_negative_heights(domain)
686   
687
688    #Extrapolate all conserved quantities
689    #if optimised_gradient_limiter:
690    #    #MH090605 if second order,
691    #    #perform the extrapolation and limiting on
692    #    #all of the conserved quantities
693
694    #    if (domain.order == 1):
695    #        for name in domain.conserved_quantities:
696    #            Q = domain.quantities[name]
697    #            Q.extrapolate_first_order()
698    #    elif domain.order == 2:
699    #        domain.extrapolate_second_order_sw()
700    #    else:
701    #        raise 'Unknown order'
702    #else:
703        #old code:
704
705    for name in domain.conserved_quantities:
706        Q = domain.quantities[name]
707        if domain.order == 1:
708            Q.extrapolate_first_order()
709        elif domain.order == 2:
710            #print "add extrapolate second order to shallow water"
711            #if name != 'height':
712            Q.extrapolate_second_order()
713            #Q.limit()
714        else:
715            raise 'Unknown order'
716
717    #Take bed elevation into account when water heights are small
718    #balance_deep_and_shallow(domain)
719    #protect_against_infinitesimal_and_negative_heights(domain)
720
721#Compute edge values by interpolation
722    #for name in domain.conserved_quantities:
723    #    Q = domain.quantities[name]
724    #    Q.interpolate_from_vertices_to_edges()   
725   
726def protect_against_infinitesimal_and_negative_heights(domain):
727    """Protect against infinitesimal heights and associated high velocities
728    """
729
730    #Shortcuts
731    wc = domain.quantities['stage'].centroid_values
732    zc = domain.quantities['elevation'].centroid_values
733    xmomc = domain.quantities['xmomentum'].centroid_values
734#    ymomc = domain.quantities['ymomentum'].centroid_values
735    hc = wc - zc  #Water depths at centroids
736
737    zv = domain.quantities['elevation'].vertex_values
738    wv = domain.quantities['stage'].vertex_values
739    hv = wv-zv
740    xmomv = domain.quantities['xmomentum'].vertex_values
741    #remove the above two lines and corresponding code below
742
743    #Update
744    for k in range(domain.number_of_elements):
745
746        if hc[k] < domain.minimum_allowed_height:
747            #Control stage
748            if hc[k] < domain.epsilon:
749                wc[k] = zc[k] # Contain 'lost mass' error
750                wv[k,0] = zv[k,0]
751                wv[k,1] = zv[k,1]
752               
753            xmomc[k] = 0.0
754           
755        #N = domain.number_of_elements
756        #if (k == 0) | (k==N-1):
757        #    wc[k] = zc[k] # Contain 'lost mass' error
758        #    wv[k,0] = zv[k,0]
759        #    wv[k,1] = zv[k,1]
760   
761def h_limiter(domain):
762    """Limit slopes for each volume to eliminate artificial variance
763    introduced by e.g. second order extrapolator
764
765    limit on h = w-z
766
767    This limiter depends on two quantities (w,z) so it resides within
768    this module rather than within quantity.py
769    """
770
771    from Numeric import zeros, Float
772
773    N = domain.number_of_elements
774    beta_h = domain.beta_h
775
776    #Shortcuts
777    wc = domain.quantities['stage'].centroid_values
778    zc = domain.quantities['elevation'].centroid_values
779    hc = wc - zc
780
781    wv = domain.quantities['stage'].vertex_values
782    zv = domain.quantities['elevation'].vertex_values
783    hv = wv-zv
784
785    hvbar = zeros(hv.shape, Float) #h-limited values
786
787    #Find min and max of this and neighbour's centroid values
788    hmax = zeros(hc.shape, Float)
789    hmin = zeros(hc.shape, Float)
790
791    for k in range(N):
792        hmax[k] = hmin[k] = hc[k]
793        #for i in range(3):
794        for i in range(2):   
795            n = domain.neighbours[k,i]
796            if n >= 0:
797                hn = hc[n] #Neighbour's centroid value
798
799                hmin[k] = min(hmin[k], hn)
800                hmax[k] = max(hmax[k], hn)
801
802
803    #Diffences between centroids and maxima/minima
804    dhmax = hmax - hc
805    dhmin = hmin - hc
806
807    #Deltas between vertex and centroid values
808    dh = zeros(hv.shape, Float)
809    #for i in range(3):
810    for i in range(2):
811        dh[:,i] = hv[:,i] - hc
812
813    #Phi limiter
814    for k in range(N):
815
816        #Find the gradient limiter (phi) across vertices
817        phi = 1.0
818        #for i in range(3):
819        for i in range(2):
820            r = 1.0
821            if (dh[k,i] > 0): r = dhmax[k]/dh[k,i]
822            if (dh[k,i] < 0): r = dhmin[k]/dh[k,i]
823
824            phi = min( min(r*beta_h, 1), phi )
825
826        #Then update using phi limiter
827        #for i in range(3):
828        for i in range(2):
829            hvbar[k,i] = hc[k] + phi*dh[k,i]
830
831    return hvbar
832
833def balance_deep_and_shallow(domain):
834    """Compute linear combination between stage as computed by
835    gradient-limiters limiting using w, and stage computed by
836    gradient-limiters limiting using h (h-limiter).
837    The former takes precedence when heights are large compared to the
838    bed slope while the latter takes precedence when heights are
839    relatively small.  Anything in between is computed as a balanced
840    linear combination in order to avoid numerical disturbances which
841    would otherwise appear as a result of hard switching between
842    modes.
843
844    The h-limiter is always applied irrespective of the order.
845    """
846
847    #Shortcuts
848    wc = domain.quantities['stage'].centroid_values
849    zc = domain.quantities['elevation'].centroid_values
850    hc = wc - zc
851
852    wv = domain.quantities['stage'].vertex_values
853    zv = domain.quantities['elevation'].vertex_values
854    hv = wv-zv
855
856    #Limit h
857    hvbar = h_limiter(domain)
858
859    for k in range(domain.number_of_elements):
860        #Compute maximal variation in bed elevation
861        #  This quantitiy is
862        #    dz = max_i abs(z_i - z_c)
863        #  and it is independent of dimension
864        #  In the 1d case zc = (z0+z1)/2
865        #  In the 2d case zc = (z0+z1+z2)/3
866
867        dz = max(abs(zv[k,0]-zc[k]),
868                 abs(zv[k,1]-zc[k]))#,
869#                 abs(zv[k,2]-zc[k]))
870
871
872        hmin = min( hv[k,:] )
873
874        #Create alpha in [0,1], where alpha==0 means using the h-limited
875        #stage and alpha==1 means using the w-limited stage as
876        #computed by the gradient limiter (both 1st or 2nd order)
877
878        #If hmin > dz/2 then alpha = 1 and the bed will have no effect
879        #If hmin < 0 then alpha = 0 reverting to constant height above bed.
880
881        if dz > 0.0:
882            alpha = max( min( 2*hmin/dz, 1.0), 0.0 )
883        else:
884            #Flat bed
885            alpha = 1.0
886
887        alpha  = 0.0
888        #Let
889        #
890        #  wvi be the w-limited stage (wvi = zvi + hvi)
891        #  wvi- be the h-limited state (wvi- = zvi + hvi-)
892        #
893        #
894        #where i=0,1,2 denotes the vertex ids
895        #
896        #Weighted balance between w-limited and h-limited stage is
897        #
898        #  wvi := (1-alpha)*(zvi+hvi-) + alpha*(zvi+hvi)
899        #
900        #It follows that the updated wvi is
901        #  wvi := zvi + (1-alpha)*hvi- + alpha*hvi
902        #
903        # Momentum is balanced between constant and limited
904
905
906        #for i in range(3):
907        #    wv[k,i] = zv[k,i] + hvbar[k,i]
908
909        #return
910
911        if alpha < 1:
912
913            #for i in range(3):
914            for i in range(2):
915                wv[k,i] = zv[k,i] + (1.0-alpha)*hvbar[k,i] + alpha*hv[k,i]
916
917            #Momentums at centroids
918            xmomc = domain.quantities['xmomentum'].centroid_values
919#            ymomc = domain.quantities['ymomentum'].centroid_values
920
921            #Momentums at vertices
922            xmomv = domain.quantities['xmomentum'].vertex_values
923#           ymomv = domain.quantities['ymomentum'].vertex_values
924
925            # Update momentum as a linear combination of
926            # xmomc and ymomc (shallow) and momentum
927            # from extrapolator xmomv and ymomv (deep).
928            xmomv[k,:] = (1.0-alpha)*xmomc[k] + alpha*xmomv[k,:]
929#            ymomv[k,:] = (1-alpha)*ymomc[k] + alpha*ymomv[k,:]
930
931
932###############################################
933#Boundaries - specific to the shallow water wave equation
934class Reflective_boundary(Boundary):
935    """Reflective boundary returns same conserved quantities as
936    those present in its neighbour volume but reflected.
937
938    This class is specific to the shallow water equation as it
939    works with the momentum quantities assumed to be the second
940    and third conserved quantities.
941    """
942
943    def __init__(self, domain = None):
944        Boundary.__init__(self)
945
946        if domain is None:
947            msg = 'Domain must be specified for reflective boundary'
948            raise msg
949
950        #Handy shorthands
951        #self.stage   = domain.quantities['stage'].edge_values
952        #self.xmom    = domain.quantities['xmomentum'].edge_values
953        #self.ymom    = domain.quantities['ymomentum'].edge_values
954        self.normals = domain.normals
955        self.stage   = domain.quantities['stage'].vertex_values
956        self.xmom    = domain.quantities['xmomentum'].vertex_values
957
958        from Numeric import zeros, Float
959        #self.conserved_quantities = zeros(3, Float)
960        self.conserved_quantities = zeros(2, Float)
961
962    def __repr__(self):
963        return 'Reflective_boundary'
964
965
966    def evaluate(self, vol_id, edge_id):
967        """Reflective boundaries reverses the outward momentum
968        of the volume they serve.
969        """
970
971        q = self.conserved_quantities
972        q[0] = self.stage[vol_id, edge_id]
973        q[1] = self.xmom[vol_id, edge_id]
974        #q[2] = self.ymom[vol_id, edge_id]
975        #normal = self.normals[vol_id, 2*edge_id:2*edge_id+2]
976        #normal = self.normals[vol_id, 2*edge_id:2*edge_id+1]
977        normal = self.normals[vol_id,edge_id]
978
979        #r = rotate(q, normal, direction = 1)
980        #r[1] = -r[1]
981        #q = rotate(r, normal, direction = -1)
982        r = q
983        r[1] = normal*r[1]
984        r[1] = -r[1]
985        r[1] = normal*r[1]
986        q = r
987        #For start interval there is no outward momentum so do not need to
988        #reverse direction in this case
989
990        return q
991
992class Dirichlet_boundary(Boundary):
993    """Dirichlet boundary returns constant values for the
994    conserved quantities
995    """
996
997
998    def __init__(self, conserved_quantities=None):
999        Boundary.__init__(self)
1000
1001        if conserved_quantities is None:
1002            msg = 'Must specify one value for each conserved quantity'
1003            raise msg
1004
1005        from Numeric import array, Float
1006        self.conserved_quantities=array(conserved_quantities).astype(Float)
1007
1008    def __repr__(self):
1009        return 'Dirichlet boundary (%s)' %self.conserved_quantities
1010
1011    def evaluate(self, vol_id=None, edge_id=None):
1012        return self.conserved_quantities
1013
1014
1015#########################
1016#Standard forcing terms:
1017#
1018def gravity(domain):
1019    """Apply gravitational pull in the presence of bed slope
1020    """
1021
1022    from util import gradient
1023    from Numeric import zeros, Float, array, sum
1024
1025    xmom = domain.quantities['xmomentum'].explicit_update
1026    stage = domain.quantities['stage'].explicit_update
1027#    ymom = domain.quantities['ymomentum'].explicit_update
1028
1029    Stage = domain.quantities['stage']
1030    Elevation = domain.quantities['elevation']
1031    #h = Stage.edge_values - Elevation.edge_values
1032    h = Stage.vertex_values - Elevation.vertex_values
1033    b = Elevation.vertex_values
1034    w = Stage.vertex_values
1035
1036    x = domain.get_vertex_coordinates()
1037    g = domain.g
1038
1039    for k in range(domain.number_of_elements):
1040#        avg_h = sum( h[k,:] )/3
1041        avg_h = sum( h[k,:] )/2
1042
1043        #Compute bed slope
1044        #x0, y0, x1, y1, x2, y2 = x[k,:]
1045        x0, x1 = x[k,:]
1046        #z0, z1, z2 = v[k,:]
1047        b0, b1 = b[k,:]
1048
1049        w0, w1 = w[k,:]
1050        wx = gradient(x0, x1, w0, w1)
1051
1052        #zx, zy = gradient(x0, y0, x1, y1, x2, y2, z0, z1, z2)
1053        bx = gradient(x0, x1, b0, b1)
1054       
1055        #Update momentum (explicit update is reset to source values)
1056        if domain.split == False:
1057            xmom[k] += -g*bx*avg_h
1058            #xmom[k] = -g*bx*avg_h
1059            #stage[k] = 0.0
1060        elif domain.split == True:
1061            xmom[k] += -g*wx*avg_h
1062            #xmom[k] = -g*wx*avg_h
1063        #ymom[k] += -g*zy*avg_h
1064 
1065def manning_friction(domain):
1066    """Apply (Manning) friction to water momentum
1067    """
1068
1069    from math import sqrt
1070
1071    w = domain.quantities['stage'].centroid_values
1072    z = domain.quantities['elevation'].centroid_values
1073    h = w-z
1074
1075    uh = domain.quantities['xmomentum'].centroid_values
1076    #vh = domain.quantities['ymomentum'].centroid_values
1077    eta = domain.quantities['friction'].centroid_values
1078
1079    xmom_update = domain.quantities['xmomentum'].semi_implicit_update
1080    #ymom_update = domain.quantities['ymomentum'].semi_implicit_update
1081
1082    N = domain.number_of_elements
1083    eps = domain.minimum_allowed_height
1084    g = domain.g
1085
1086    for k in range(N):
1087        if eta[k] >= eps:
1088            if h[k] >= eps:
1089                #S = -g * eta[k]**2 * sqrt((uh[k]**2 + vh[k]**2))
1090                S = -g * eta[k]**2 * uh[k]
1091                S /= h[k]**(7.0/3)
1092
1093                #Update momentum
1094                xmom_update[k] += S*uh[k]
1095                #ymom_update[k] += S*vh[k]
1096
1097def linear_friction(domain):
1098    """Apply linear friction to water momentum
1099
1100    Assumes quantity: 'linear_friction' to be present
1101    """
1102
1103    from math import sqrt
1104
1105    w = domain.quantities['stage'].centroid_values
1106    z = domain.quantities['elevation'].centroid_values
1107    h = w-z
1108
1109    uh = domain.quantities['xmomentum'].centroid_values
1110#    vh = domain.quantities['ymomentum'].centroid_values
1111    tau = domain.quantities['linear_friction'].centroid_values
1112
1113    xmom_update = domain.quantities['xmomentum'].semi_implicit_update
1114#    ymom_update = domain.quantities['ymomentum'].semi_implicit_update
1115
1116    N = domain.number_of_elements
1117    eps = domain.minimum_allowed_height
1118    g = domain.g #Not necessary? Why was this added?
1119
1120    for k in range(N):
1121        if tau[k] >= eps:
1122            if h[k] >= eps:
1123                S = -tau[k]/h[k]
1124
1125                #Update momentum
1126                xmom_update[k] += S*uh[k]
1127 #              ymom_update[k] += S*vh[k]
1128
1129
1130
1131def check_forcefield(f):
1132    """Check that f is either
1133    1: a callable object f(t,x,y), where x and y are vectors
1134       and that it returns an array or a list of same length
1135       as x and y
1136    2: a scalar
1137    """
1138
1139    from Numeric import ones, Float, array
1140
1141
1142    if callable(f):
1143        #N = 3
1144        N = 2
1145        #x = ones(3, Float)
1146        #y = ones(3, Float)
1147        x = ones(2, Float)
1148        #y = ones(2, Float)
1149       
1150        try:
1151            #q = f(1.0, x=x, y=y)
1152            q = f(1.0, x=x)
1153        except Exception, e:
1154            msg = 'Function %s could not be executed:\n%s' %(f, e)
1155            #FIXME: Reconsider this semantics
1156            raise msg
1157
1158        try:
1159            q = array(q).astype(Float)
1160        except:
1161            msg = 'Return value from vector function %s could ' %f
1162            msg += 'not be converted into a Numeric array of floats.\n'
1163            msg += 'Specified function should return either list or array.'
1164            raise msg
1165
1166        #Is this really what we want?
1167        msg = 'Return vector from function %s ' %f
1168        msg += 'must have same lenght as input vectors'
1169        assert len(q) == N, msg
1170
1171    else:
1172        try:
1173            f = float(f)
1174        except:
1175            msg = 'Force field %s must be either a scalar' %f
1176            msg += ' or a vector function'
1177            raise msg
1178    return f
1179
1180class Wind_stress:
1181    """Apply wind stress to water momentum in terms of
1182    wind speed [m/s] and wind direction [degrees]
1183    """
1184
1185    def __init__(self, *args, **kwargs):
1186        """Initialise windfield from wind speed s [m/s]
1187        and wind direction phi [degrees]
1188
1189        Inputs v and phi can be either scalars or Python functions, e.g.
1190
1191        W = Wind_stress(10, 178)
1192
1193        #FIXME - 'normal' degrees are assumed for now, i.e. the
1194        vector (1,0) has zero degrees.
1195        We may need to convert from 'compass' degrees later on and also
1196        map from True north to grid north.
1197
1198        Arguments can also be Python functions of t,x,y as in
1199
1200        def speed(t,x,y):
1201            ...
1202            return s
1203
1204        def angle(t,x,y):
1205            ...
1206            return phi
1207
1208        where x and y are vectors.
1209
1210        and then pass the functions in
1211
1212        W = Wind_stress(speed, angle)
1213
1214        The instantiated object W can be appended to the list of
1215        forcing_terms as in
1216
1217        Alternatively, one vector valued function for (speed, angle)
1218        can be applied, providing both quantities simultaneously.
1219        As in
1220        W = Wind_stress(F), where returns (speed, angle) for each t.
1221
1222        domain.forcing_terms.append(W)
1223        """
1224
1225        from config import rho_a, rho_w, eta_w
1226        from Numeric import array, Float
1227
1228        if len(args) == 2:
1229            s = args[0]
1230            phi = args[1]
1231        elif len(args) == 1:
1232            #Assume vector function returning (s, phi)(t,x,y)
1233            vector_function = args[0]
1234            #s = lambda t,x,y: vector_function(t,x=x,y=y)[0]
1235            #phi = lambda t,x,y: vector_function(t,x=x,y=y)[1]
1236            s = lambda t,x: vector_function(t,x=x)[0]
1237            phi = lambda t,x: vector_function(t,x=x)[1]
1238        else:
1239           #Assume info is in 2 keyword arguments
1240
1241           if len(kwargs) == 2:
1242               s = kwargs['s']
1243               phi = kwargs['phi']
1244           else:
1245               raise 'Assumes two keyword arguments: s=..., phi=....'
1246
1247        print 'phi', phi
1248        self.speed = check_forcefield(s)
1249        self.phi = check_forcefield(phi)
1250
1251        self.const = eta_w*rho_a/rho_w
1252
1253
1254    def __call__(self, domain):
1255        """Evaluate windfield based on values found in domain
1256        """
1257
1258        from math import pi, cos, sin, sqrt
1259        from Numeric import Float, ones, ArrayType
1260
1261        xmom_update = domain.quantities['xmomentum'].explicit_update
1262        #ymom_update = domain.quantities['ymomentum'].explicit_update
1263
1264        N = domain.number_of_elements
1265        t = domain.time
1266
1267        if callable(self.speed):
1268            xc = domain.get_centroid_coordinates()
1269            #s_vec = self.speed(t, xc[:,0], xc[:,1])
1270            s_vec = self.speed(t, xc)
1271        else:
1272            #Assume s is a scalar
1273
1274            try:
1275                s_vec = self.speed * ones(N, Float)
1276            except:
1277                msg = 'Speed must be either callable or a scalar: %s' %self.s
1278                raise msg
1279
1280
1281        if callable(self.phi):
1282            xc = domain.get_centroid_coordinates()
1283            #phi_vec = self.phi(t, xc[:,0], xc[:,1])
1284            phi_vec = self.phi(t, xc)
1285        else:
1286            #Assume phi is a scalar
1287
1288            try:
1289                phi_vec = self.phi * ones(N, Float)
1290            except:
1291                msg = 'Angle must be either callable or a scalar: %s' %self.phi
1292                raise msg
1293
1294        #assign_windfield_values(xmom_update, ymom_update,
1295        #                        s_vec, phi_vec, self.const)
1296        assign_windfield_values(xmom_update, s_vec, phi_vec, self.const)
1297
1298
1299#def assign_windfield_values(xmom_update, ymom_update,
1300#                            s_vec, phi_vec, const):
1301def assign_windfield_values(xmom_update, s_vec, phi_vec, const):
1302    """Python version of assigning wind field to update vectors.
1303    A c version also exists (for speed)
1304    """
1305    from math import pi, cos, sin, sqrt
1306
1307    N = len(s_vec)
1308    for k in range(N):
1309        s = s_vec[k]
1310        phi = phi_vec[k]
1311
1312        #Convert to radians
1313        phi = phi*pi/180
1314
1315        #Compute velocity vector (u, v)
1316        u = s*cos(phi)
1317        v = s*sin(phi)
1318
1319        #Compute wind stress
1320        #S = const * sqrt(u**2 + v**2)
1321        S = const * u
1322        xmom_update[k] += S*u
1323        #ymom_update[k] += S*v
Note: See TracBrowser for help on using the repository browser.