source: inundation/ga/storm_surge/pyvolution/shallow_water.py @ 771

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