source: inundation/pyvolution/shallow_water.py @ 2648

Last change on this file since 2648 was 2648, checked in by steve, 18 years ago

Looking at island.py example. I changed the order of the implicit and explicit update (in quantity_ext.c) which I think improved the situation a little. But this has left a few fails in the test_all suite which we need to fix up. Mainly small differences in results.
  • Property svn:keywords set to LastChangedDate LastChangedRevision LastChangedBy HeadURL Id
File size: 56.4 KB
Line 
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
9U_t + E_x + G_y = S
10
11where
12
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, ...)
18
19and _t, _x, _y denote the derivative with respect to t, x and y respectively.
20
21The quantities are
22
23symbol    variable name    explanation
24x         x                horizontal distance from origin [m]
25y         y                vertical distance from origin [m]
26z         elevation        elevation of bed on which flow is modelled [m]
27h         height           water height above z [m]
28w         stage            absolute water level, w = z+h [m]
29u                          speed in the x direction [m/s]
30v                          speed in the y direction [m/s]
31uh        xmomentum        momentum in the x direction [m^2/s]
32vh        ymomentum        momentum in the y direction [m^2/s]
33
34eta                        mannings friction coefficient [to appear]
35nu                         wind stress coefficient [to appear]
36
37The conserved quantities are w, uh, vh
38
39
40For details see e.g.
41Christopher Zoppou and Stephen Roberts,
42Catastrophic Collapse of Water Supply Reservoirs in Urban Areas,
43Journal of Hydraulic Engineering, vol. 127, No. 7 July 1999
44
45
46Ole Nielsen, Stephen Roberts, Duncan Gray, Christopher Zoppou
47Geoscience Australia, 2004
48"""
49
50#Subversion keywords:
51#
52#$LastChangedDate: 2006-03-31 07:21:26 +0000 (Fri, 31 Mar 2006) $
53#$LastChangedRevision: 2648 $
54#$LastChangedBy: steve $
55
56
57from domain import Domain as Generic_Domain
58from generic_boundary_conditions import Boundary
59from generic_boundary_conditions import File_boundary
60from generic_boundary_conditions import Dirichlet_boundary
61from generic_boundary_conditions import Time_boundary
62from generic_boundary_conditions import Transmissive_boundary
63
64from utilities.numerical_tools import gradient, mean
65
66
67
68#Shallow water domain
69class Domain(Generic_Domain):
70
71    def __init__(self, coordinates=None, vertices=None, boundary = None,
72                 tagged_elements = None, geo_reference = None,
73                 use_inscribed_circle=False,
74                 mesh_filename=None):
75
76        conserved_quantities = ['stage', 'xmomentum', 'ymomentum']
77        other_quantities = ['elevation', 'friction']
78        Generic_Domain.__init__(self, coordinates, vertices, boundary,
79                                conserved_quantities, other_quantities,
80                                tagged_elements, geo_reference,
81                                use_inscribed_circle, mesh_filename)
82
83        from config import minimum_allowed_height, maximum_allowed_speed, g
84        self.minimum_allowed_height = minimum_allowed_height
85        self.maximum_allowed_speed = maximum_allowed_speed
86        self.g = g
87
88
89        self.forcing_terms.append(manning_friction)
90        self.forcing_terms.append(gravity)
91
92        #Realtime visualisation
93        self.visualiser = None
94        self.visualise  = False
95        self.visualise_color_stage = False
96        self.visualise_stage_range = 1.0
97        self.visualise_timer = True
98        self.visualise_range_z = None
99
100        #Stored output
101        self.store = True
102        self.format = 'sww'
103        self.smooth = True
104
105        #Reduction operation for get_vertex_values
106        self.reduction = mean
107        #self.reduction = min  #Looks better near steep slopes
108
109        self.quantities_to_be_stored = ['stage','xmomentum','ymomentum']
110
111
112    def set_quantities_to_be_stored(self, q):
113        """Specify which quantities will be stored in the sww file.
114
115        q must be either:
116          - the name of a quantity
117          - a list of quantity names
118          - None
119
120        In the two first cases, the named quantities will be stored at each yieldstep
121        (This is in addition to the quantities elevation and friction)
122        If q is None, storage will be switched off altogether.
123        """
124
125
126        if q is None:
127            self.quantities_to_be_stored = []
128            self.store = False
129            return
130
131        if isinstance(q, basestring):
132            q = [q] # Turn argument into a list
133
134        #Check correcness
135        for quantity_name in q:
136            msg = 'Quantity %s is not a valid conserved quantity' %quantity_name
137            assert quantity_name in self.conserved_quantities, msg
138
139        self.quantities_to_be_stored = q
140
141
142    def initialise_visualiser(self,scale_z=1.0,rect=None):
143        #Realtime visualisation
144        if self.visualiser is None:
145            from realtime_visualisation_new import Visualiser
146            self.visualiser = Visualiser(self,scale_z,rect)
147            self.visualiser.setup['elevation']=True
148            self.visualiser.updating['stage']=True
149        self.visualise = True
150        if self.visualise_color_stage == True:
151            self.visualiser.coloring['stage'] = True
152            self.visualiser.qcolor['stage'] = (0.0, 0.0, 0.8)
153        print 'initialise visualiser'
154        print self.visualiser.setup
155        print self.visualiser.updating
156
157    def check_integrity(self):
158        Generic_Domain.check_integrity(self)
159
160        #Check that we are solving the shallow water wave equation
161
162        msg = 'First conserved quantity must be "stage"'
163        assert self.conserved_quantities[0] == 'stage', msg
164        msg = 'Second conserved quantity must be "xmomentum"'
165        assert self.conserved_quantities[1] == 'xmomentum', msg
166        msg = 'Third conserved quantity must be "ymomentum"'
167        assert self.conserved_quantities[2] == 'ymomentum', msg
168
169    def extrapolate_second_order_sw(self):
170        #Call correct module function
171        #(either from this module or C-extension)
172        extrapolate_second_order_sw(self)
173
174    def compute_fluxes(self):
175        #Call correct module function
176        #(either from this module or C-extension)
177        compute_fluxes(self)
178
179    def distribute_to_vertices_and_edges(self):
180        #Call correct module function
181        #(either from this module or C-extension)
182        distribute_to_vertices_and_edges(self)
183
184
185    #FIXME: Under construction
186#     def set_defaults(self):
187#         """Set default values for uninitialised quantities.
188#         This is specific to the shallow water wave equation
189#         Defaults for 'elevation', 'friction', 'xmomentum' and 'ymomentum'
190#         are 0.0. Default for 'stage' is whatever the value of 'elevation'.
191#         """
192
193#         for name in self.other_quantities + self.conserved_quantities:
194#             print name
195#             print self.quantities.keys()
196#             if not self.quantities.has_key(name):
197#                 if name == 'stage':
198
199#                     if self.quantities.has_key('elevation'):
200#                         z = self.quantities['elevation'].vertex_values
201#                         self.set_quantity(name, z)
202#                     else:
203#                         self.set_quantity(name, 0.0)
204#                 else:
205#                     self.set_quantity(name, 0.0)
206
207
208
209#         #Lift negative heights up
210#         #z = self.quantities['elevation'].vertex_values
211#         #w = self.quantities['stage'].vertex_values
212
213#         #h = w-z
214
215#         #for k in range(h.shape[0]):
216#         #    for i in range(3):
217#         #        if h[k, i] < 0.0:
218#         #            w[k, i] = z[k, i]
219
220
221#         #self.quantities['stage'].interpolate()
222
223
224    def evolve(self,
225               yieldstep = None,
226               finaltime = None,
227               duration = None,
228               skip_initial_step = False):
229        """Specialisation of basic evolve method from parent class
230        """
231
232        #Call check integrity here rather than from user scripts
233        #self.check_integrity()
234
235        msg = 'Parameter beta_h must be in the interval [0, 1['
236        assert 0 <= self.beta_h < 1.0, msg
237        msg = 'Parameter beta_w must be in the interval [0, 1['
238        assert 0 <= self.beta_w < 1.0, msg
239
240
241        #Initial update of vertex and edge values before any storage
242        #and or visualisation
243        self.distribute_to_vertices_and_edges()
244
245        #Initialise real time viz if requested
246        if self.visualise is True and self.time == 0.0:
247            if self.visualiser is None:
248                self.initialise_visualiser()
249
250            self.visualiser.update_timer()
251            self.visualiser.setup_all()
252
253        #Store model data, e.g. for visualisation
254        if self.store is True and self.time == 0.0:
255            self.initialise_storage()
256            #print 'Storing results in ' + self.writer.filename
257        else:
258            pass
259            #print 'Results will not be stored.'
260            #print 'To store results set domain.store = True'
261            #FIXME: Diagnostic output should be controlled by
262            # a 'verbose' flag living in domain (or in a parent class)
263
264        #Call basic machinery from parent class
265        for t in Generic_Domain.evolve(self,
266                                       yieldstep=yieldstep,
267                                       finaltime=finaltime,
268                                       duration=duration,
269                                       skip_initial_step=skip_initial_step):
270            #Real time viz
271            if self.visualise is True:
272                self.visualiser.update_all()
273                self.visualiser.update_timer()
274
275
276            #Store model data, e.g. for subsequent visualisation
277            if self.store is True:
278                self.store_timestep(self.quantities_to_be_stored)
279
280            #FIXME: Could maybe be taken from specified list
281            #of 'store every step' quantities
282
283            #Pass control on to outer loop for more specific actions
284            yield(t)
285
286    def initialise_storage(self):
287        """Create and initialise self.writer object for storing data.
288        Also, save x,y and bed elevation
289        """
290
291        import data_manager
292
293        #Initialise writer
294        self.writer = data_manager.get_dataobject(self, mode = 'w')
295
296        #Store vertices and connectivity
297        self.writer.store_connectivity()
298
299
300    def store_timestep(self, name):
301        """Store named quantity and time.
302
303        Precondition:
304           self.write has been initialised
305        """
306        self.writer.store_timestep(name)
307
308
309#=============== End of Shallow Water Domain ===============================
310
311
312
313#Rotation of momentum vector
314def rotate(q, normal, direction = 1):
315    """Rotate the momentum component q (q[1], q[2])
316    from x,y coordinates to coordinates based on normal vector.
317
318    If direction is negative the rotation is inverted.
319
320    Input vector is preserved
321
322    This function is specific to the shallow water wave equation
323    """
324
325    from Numeric import zeros, Float
326
327    assert len(q) == 3,\
328           'Vector of conserved quantities must have length 3'\
329           'for 2D shallow water equation'
330
331    try:
332        l = len(normal)
333    except:
334        raise 'Normal vector must be an Numeric array'
335
336    assert l == 2, 'Normal vector must have 2 components'
337
338
339    n1 = normal[0]
340    n2 = normal[1]
341
342    r = zeros(len(q), Float) #Rotated quantities
343    r[0] = q[0]              #First quantity, height, is not rotated
344
345    if direction == -1:
346        n2 = -n2
347
348
349    r[1] =  n1*q[1] + n2*q[2]
350    r[2] = -n2*q[1] + n1*q[2]
351
352    return r
353
354
355
356####################################
357# Flux computation
358def flux_function(normal, ql, qr, zl, zr):
359    """Compute fluxes between volumes for the shallow water wave equation
360    cast in terms of w = h+z using the 'central scheme' as described in
361
362    Kurganov, Noelle, Petrova. 'Semidiscrete Central-Upwind Schemes For
363    Hyperbolic Conservation Laws and Hamilton-Jacobi Equations'.
364    Siam J. Sci. Comput. Vol. 23, No. 3, pp. 707-740.
365
366    The implemented formula is given in equation (3.15) on page 714
367
368    Conserved quantities w, uh, vh are stored as elements 0, 1 and 2
369    in the numerical vectors ql and qr.
370
371    Bed elevations zl and zr.
372    """
373
374    from config import g, epsilon
375    from math import sqrt
376    from Numeric import array
377
378    #Align momentums with x-axis
379    q_left  = rotate(ql, normal, direction = 1)
380    q_right = rotate(qr, normal, direction = 1)
381
382    z = (zl+zr)/2 #Take average of field values
383
384    w_left  = q_left[0]   #w=h+z
385    h_left  = w_left-z
386    uh_left = q_left[1]
387
388    if h_left < epsilon:
389        u_left = 0.0  #Could have been negative
390        h_left = 0.0
391    else:
392        u_left  = uh_left/h_left
393
394
395    w_right  = q_right[0]  #w=h+z
396    h_right  = w_right-z
397    uh_right = q_right[1]
398
399
400    if h_right < epsilon:
401        u_right = 0.0  #Could have been negative
402        h_right = 0.0
403    else:
404        u_right  = uh_right/h_right
405
406    vh_left  = q_left[2]
407    vh_right = q_right[2]
408
409    soundspeed_left  = sqrt(g*h_left)
410    soundspeed_right = sqrt(g*h_right)
411
412    #Maximal wave speed
413    s_max = max(u_left+soundspeed_left, u_right+soundspeed_right, 0)
414
415    #Minimal wave speed
416    s_min = min(u_left-soundspeed_left, u_right-soundspeed_right, 0)
417
418    #Flux computation
419
420    #FIXME(Ole): Why is it again that we don't
421    #use uh_left and uh_right directly in the first entries?
422    flux_left  = array([u_left*h_left,
423                        u_left*uh_left + 0.5*g*h_left**2,
424                        u_left*vh_left])
425    flux_right = array([u_right*h_right,
426                        u_right*uh_right + 0.5*g*h_right**2,
427                        u_right*vh_right])
428
429    denom = s_max-s_min
430    if denom == 0.0:
431        edgeflux = array([0.0, 0.0, 0.0])
432        max_speed = 0.0
433    else:
434        edgeflux = (s_max*flux_left - s_min*flux_right)/denom
435        edgeflux += s_max*s_min*(q_right-q_left)/denom
436
437        edgeflux = rotate(edgeflux, normal, direction=-1)
438        max_speed = max(abs(s_max), abs(s_min))
439
440    return edgeflux, max_speed
441
442
443def compute_fluxes(domain):
444    """Compute all fluxes and the timestep suitable for all volumes
445    in domain.
446
447    Compute total flux for each conserved quantity using "flux_function"
448
449    Fluxes across each edge are scaled by edgelengths and summed up
450    Resulting flux is then scaled by area and stored in
451    explicit_update for each of the three conserved quantities
452    stage, xmomentum and ymomentum
453
454    The maximal allowable speed computed by the flux_function for each volume
455    is converted to a timestep that must not be exceeded. The minimum of
456    those is computed as the next overall timestep.
457
458    Post conditions:
459      domain.explicit_update is reset to computed flux values
460      domain.timestep is set to the largest step satisfying all volumes.
461    """
462
463    import sys
464    from Numeric import zeros, Float
465
466    N = domain.number_of_elements
467
468    #Shortcuts
469    Stage = domain.quantities['stage']
470    Xmom = domain.quantities['xmomentum']
471    Ymom = domain.quantities['ymomentum']
472    Bed = domain.quantities['elevation']
473
474    #Arrays
475    stage = Stage.edge_values
476    xmom =  Xmom.edge_values
477    ymom =  Ymom.edge_values
478    bed =   Bed.edge_values
479
480    stage_bdry = Stage.boundary_values
481    xmom_bdry =  Xmom.boundary_values
482    ymom_bdry =  Ymom.boundary_values
483
484    flux = zeros(3, Float) #Work array for summing up fluxes
485
486
487    #Loop
488    timestep = float(sys.maxint)
489    for k in range(N):
490
491        flux[:] = 0.  #Reset work array
492        for i in range(3):
493            #Quantities inside volume facing neighbour i
494            ql = [stage[k, i], xmom[k, i], ymom[k, i]]
495            zl = bed[k, i]
496
497            #Quantities at neighbour on nearest face
498            n = domain.neighbours[k,i]
499            if n < 0:
500                m = -n-1 #Convert negative flag to index
501                qr = [stage_bdry[m], xmom_bdry[m], ymom_bdry[m]]
502                zr = zl #Extend bed elevation to boundary
503            else:
504                m = domain.neighbour_edges[k,i]
505                qr = [stage[n, m], xmom[n, m], ymom[n, m]]
506                zr = bed[n, m]
507
508
509            #Outward pointing normal vector
510            normal = domain.normals[k, 2*i:2*i+2]
511
512            #Flux computation using provided function
513            edgeflux, max_speed = flux_function(normal, ql, qr, zl, zr)
514            flux -= edgeflux * domain.edgelengths[k,i]
515
516            #Update optimal_timestep
517            try:
518                timestep = min(timestep, 0.5*domain.radii[k]/max_speed)
519            except ZeroDivisionError:
520                pass
521
522        #Normalise by area and store for when all conserved
523        #quantities get updated
524        flux /= domain.areas[k]
525        Stage.explicit_update[k] = flux[0]
526        Xmom.explicit_update[k] = flux[1]
527        Ymom.explicit_update[k] = flux[2]
528
529
530    domain.timestep = timestep
531
532#MH090605 The following method belongs to the shallow_water domain class
533#see comments in the corresponding method in shallow_water_ext.c
534def extrapolate_second_order_sw_c(domain):
535    """Wrapper calling C version of extrapolate_second_order_sw
536    """
537    import sys
538    from Numeric import zeros, Float
539
540    N = domain.number_of_elements
541
542    #Shortcuts
543    Stage = domain.quantities['stage']
544    Xmom = domain.quantities['xmomentum']
545    Ymom = domain.quantities['ymomentum']
546    from shallow_water_ext import extrapolate_second_order_sw
547    extrapolate_second_order_sw(domain,domain.surrogate_neighbours,
548                                domain.number_of_boundaries,
549                                domain.centroid_coordinates,
550                                Stage.centroid_values,
551                                Xmom.centroid_values,
552                                Ymom.centroid_values,
553                                domain.vertex_coordinates,
554                                Stage.vertex_values,
555                                Xmom.vertex_values,
556                                Ymom.vertex_values)
557
558def compute_fluxes_c(domain):
559    """Wrapper calling C version of compute fluxes
560    """
561
562    import sys
563    from Numeric import zeros, Float
564
565    N = domain.number_of_elements
566
567    #Shortcuts
568    Stage = domain.quantities['stage']
569    Xmom = domain.quantities['xmomentum']
570    Ymom = domain.quantities['ymomentum']
571    Bed = domain.quantities['elevation']
572
573    timestep = float(sys.maxint)
574    from shallow_water_ext import compute_fluxes
575    domain.timestep = compute_fluxes(timestep, domain.epsilon, domain.g,
576                                     domain.neighbours,
577                                     domain.neighbour_edges,
578                                     domain.normals,
579                                     domain.edgelengths,
580                                     domain.radii,
581                                     domain.areas,
582                                     Stage.edge_values,
583                                     Xmom.edge_values,
584                                     Ymom.edge_values,
585                                     Bed.edge_values,
586                                     Stage.boundary_values,
587                                     Xmom.boundary_values,
588                                     Ymom.boundary_values,
589                                     Stage.explicit_update,
590                                     Xmom.explicit_update,
591                                     Ymom.explicit_update,
592                                     domain.already_computed_flux)
593
594
595####################################
596# Module functions for gradient limiting (distribute_to_vertices_and_edges)
597
598def distribute_to_vertices_and_edges(domain):
599    """Distribution from centroids to vertices specific to the
600    shallow water wave
601    equation.
602
603    It will ensure that h (w-z) is always non-negative even in the
604    presence of steep bed-slopes by taking a weighted average between shallow
605    and deep cases.
606
607    In addition, all conserved quantities get distributed as per either a
608    constant (order==1) or a piecewise linear function (order==2).
609
610    FIXME: more explanation about removal of artificial variability etc
611
612    Precondition:
613      All quantities defined at centroids and bed elevation defined at
614      vertices.
615
616    Postcondition
617      Conserved quantities defined at vertices
618
619    """
620
621    from config import optimised_gradient_limiter
622
623    #Remove very thin layers of water
624    protect_against_infinitesimal_and_negative_heights(domain)
625
626    #Extrapolate all conserved quantities
627    if optimised_gradient_limiter:
628        #MH090605 if second order,
629        #perform the extrapolation and limiting on
630        #all of the conserved quantitie
631
632        if (domain.order == 1):
633            for name in domain.conserved_quantities:
634                Q = domain.quantities[name]
635                Q.extrapolate_first_order()
636        elif domain.order == 2:
637            domain.extrapolate_second_order_sw()
638        else:
639            raise 'Unknown order'
640    else:
641        #old code:
642        for name in domain.conserved_quantities:
643            Q = domain.quantities[name]
644            if domain.order == 1:
645                Q.extrapolate_first_order()
646            elif domain.order == 2:
647                Q.extrapolate_second_order()
648                Q.limit()
649            else:
650                raise 'Unknown order'
651
652
653    #Take bed elevation into account when water heights are small
654    balance_deep_and_shallow(domain)
655
656    #Compute edge values by interpolation
657    for name in domain.conserved_quantities:
658        Q = domain.quantities[name]
659        Q.interpolate_from_vertices_to_edges()
660
661
662def protect_against_infinitesimal_and_negative_heights(domain):
663    """Protect against infinitesimal heights and associated high velocities
664    """
665
666    #Shortcuts
667    wc = domain.quantities['stage'].centroid_values
668    zc = domain.quantities['elevation'].centroid_values
669    xmomc = domain.quantities['xmomentum'].centroid_values
670    ymomc = domain.quantities['ymomentum'].centroid_values
671    hc = wc - zc  #Water depths at centroids
672
673    #Update
674    #FIXME: Modify accroditg to c-version - or discard altogether.
675    for k in range(domain.number_of_elements):
676
677        if hc[k] < domain.minimum_allowed_height:
678            #Control stage
679            if hc[k] < domain.epsilon:
680                wc[k] = zc[k] # Contain 'lost mass' error
681
682            #Control momentum
683            xmomc[k] = ymomc[k] = 0.0
684
685
686def protect_against_infinitesimal_and_negative_heights_c(domain):
687    """Protect against infinitesimal heights and associated high velocities
688    """
689
690    #Shortcuts
691    wc = domain.quantities['stage'].centroid_values
692    zc = domain.quantities['elevation'].centroid_values
693    xmomc = domain.quantities['xmomentum'].centroid_values
694    ymomc = domain.quantities['ymomentum'].centroid_values
695
696    from shallow_water_ext import protect
697
698    protect(domain.minimum_allowed_height, domain.maximum_allowed_speed,
699            domain.epsilon, wc, zc, xmomc, ymomc)
700
701
702
703def h_limiter(domain):
704    """Limit slopes for each volume to eliminate artificial variance
705    introduced by e.g. second order extrapolator
706
707    limit on h = w-z
708
709    This limiter depends on two quantities (w,z) so it resides within
710    this module rather than within quantity.py
711    """
712
713    from Numeric import zeros, Float
714
715    N = domain.number_of_elements
716    beta_h = domain.beta_h
717
718    #Shortcuts
719    wc = domain.quantities['stage'].centroid_values
720    zc = domain.quantities['elevation'].centroid_values
721    hc = wc - zc
722
723    wv = domain.quantities['stage'].vertex_values
724    zv = domain.quantities['elevation'].vertex_values
725    hv = wv-zv
726
727    hvbar = zeros(hv.shape, Float) #h-limited values
728
729    #Find min and max of this and neighbour's centroid values
730    hmax = zeros(hc.shape, Float)
731    hmin = zeros(hc.shape, Float)
732
733    for k in range(N):
734        hmax[k] = hmin[k] = hc[k]
735        for i in range(3):
736            n = domain.neighbours[k,i]
737            if n >= 0:
738                hn = hc[n] #Neighbour's centroid value
739
740                hmin[k] = min(hmin[k], hn)
741                hmax[k] = max(hmax[k], hn)
742
743
744    #Diffences between centroids and maxima/minima
745    dhmax = hmax - hc
746    dhmin = hmin - hc
747
748    #Deltas between vertex and centroid values
749    dh = zeros(hv.shape, Float)
750    for i in range(3):
751        dh[:,i] = hv[:,i] - hc
752
753    #Phi limiter
754    for k in range(N):
755
756        #Find the gradient limiter (phi) across vertices
757        phi = 1.0
758        for i in range(3):
759            r = 1.0
760            if (dh[k,i] > 0): r = dhmax[k]/dh[k,i]
761            if (dh[k,i] < 0): r = dhmin[k]/dh[k,i]
762
763            phi = min( min(r*beta_h, 1), phi )
764
765        #Then update using phi limiter
766        for i in range(3):
767            hvbar[k,i] = hc[k] + phi*dh[k,i]
768
769    return hvbar
770
771
772
773def h_limiter_c(domain):
774    """Limit slopes for each volume to eliminate artificial variance
775    introduced by e.g. second order extrapolator
776
777    limit on h = w-z
778
779    This limiter depends on two quantities (w,z) so it resides within
780    this module rather than within quantity.py
781
782    Wrapper for c-extension
783    """
784
785    from Numeric import zeros, Float
786
787    N = domain.number_of_elements
788    beta_h = domain.beta_h
789
790    #Shortcuts
791    wc = domain.quantities['stage'].centroid_values
792    zc = domain.quantities['elevation'].centroid_values
793    hc = wc - zc
794
795    wv = domain.quantities['stage'].vertex_values
796    zv = domain.quantities['elevation'].vertex_values
797    hv = wv - zv
798
799    #Call C-extension
800    from shallow_water_ext import h_limiter_sw as h_limiter
801    hvbar = h_limiter(domain, hc, hv)
802
803    return hvbar
804
805
806def balance_deep_and_shallow(domain):
807    """Compute linear combination between stage as computed by
808    gradient-limiters limiting using w, and stage computed by
809    gradient-limiters limiting using h (h-limiter).
810    The former takes precedence when heights are large compared to the
811    bed slope while the latter takes precedence when heights are
812    relatively small.  Anything in between is computed as a balanced
813    linear combination in order to avoid numerical disturbances which
814    would otherwise appear as a result of hard switching between
815    modes.
816
817    The h-limiter is always applied irrespective of the order.
818    """
819
820    #Shortcuts
821    wc = domain.quantities['stage'].centroid_values
822    zc = domain.quantities['elevation'].centroid_values
823    hc = wc - zc
824
825    wv = domain.quantities['stage'].vertex_values
826    zv = domain.quantities['elevation'].vertex_values
827    hv = wv-zv
828
829    #Limit h
830    hvbar = h_limiter(domain)
831
832    for k in range(domain.number_of_elements):
833        #Compute maximal variation in bed elevation
834        #  This quantitiy is
835        #    dz = max_i abs(z_i - z_c)
836        #  and it is independent of dimension
837        #  In the 1d case zc = (z0+z1)/2
838        #  In the 2d case zc = (z0+z1+z2)/3
839
840        dz = max(abs(zv[k,0]-zc[k]),
841                 abs(zv[k,1]-zc[k]),
842                 abs(zv[k,2]-zc[k]))
843
844
845        hmin = min( hv[k,:] )
846
847        #Create alpha in [0,1], where alpha==0 means using the h-limited
848        #stage and alpha==1 means using the w-limited stage as
849        #computed by the gradient limiter (both 1st or 2nd order)
850
851        #If hmin > dz/2 then alpha = 1 and the bed will have no effect
852        #If hmin < 0 then alpha = 0 reverting to constant height above bed.
853
854        if dz > 0.0:
855            alpha = max( min( 2*hmin/dz, 1.0), 0.0 )
856        else:
857            #Flat bed
858            alpha = 1.0
859
860        #Let
861        #
862        #  wvi be the w-limited stage (wvi = zvi + hvi)
863        #  wvi- be the h-limited state (wvi- = zvi + hvi-)
864        #
865        #
866        #where i=0,1,2 denotes the vertex ids
867        #
868        #Weighted balance between w-limited and h-limited stage is
869        #
870        #  wvi := (1-alpha)*(zvi+hvi-) + alpha*(zvi+hvi)
871        #
872        #It follows that the updated wvi is
873        #  wvi := zvi + (1-alpha)*hvi- + alpha*hvi
874        #
875        # Momentum is balanced between constant and limited
876
877
878        #for i in range(3):
879        #    wv[k,i] = zv[k,i] + hvbar[k,i]
880
881        #return
882
883        if alpha < 1:
884
885            for i in range(3):
886                wv[k,i] = zv[k,i] + (1-alpha)*hvbar[k,i] + alpha*hv[k,i]
887
888            #Momentums at centroids
889            xmomc = domain.quantities['xmomentum'].centroid_values
890            ymomc = domain.quantities['ymomentum'].centroid_values
891
892            #Momentums at vertices
893            xmomv = domain.quantities['xmomentum'].vertex_values
894            ymomv = domain.quantities['ymomentum'].vertex_values
895
896            # Update momentum as a linear combination of
897            # xmomc and ymomc (shallow) and momentum
898            # from extrapolator xmomv and ymomv (deep).
899            xmomv[k,:] = (1-alpha)*xmomc[k] + alpha*xmomv[k,:]
900            ymomv[k,:] = (1-alpha)*ymomc[k] + alpha*ymomv[k,:]
901
902
903def balance_deep_and_shallow_c(domain):
904    """Wrapper for C implementation
905    """
906
907    #Shortcuts
908    wc = domain.quantities['stage'].centroid_values
909    zc = domain.quantities['elevation'].centroid_values
910    hc = wc - zc
911
912    wv = domain.quantities['stage'].vertex_values
913    zv = domain.quantities['elevation'].vertex_values
914    hv = wv - zv
915
916    #Momentums at centroids
917    xmomc = domain.quantities['xmomentum'].centroid_values
918    ymomc = domain.quantities['ymomentum'].centroid_values
919
920    #Momentums at vertices
921    xmomv = domain.quantities['xmomentum'].vertex_values
922    ymomv = domain.quantities['ymomentum'].vertex_values
923
924    #Limit h
925    hvbar = h_limiter(domain)
926
927    #This is how one would make a first order h_limited value
928    #as in the old balancer (pre 17 Feb 2005):
929    #from Numeric import zeros, Float
930    #hvbar = zeros( (len(hc), 3), Float)
931    #for i in range(3):
932    #    hvbar[:,i] = hc[:]
933
934    from shallow_water_ext import balance_deep_and_shallow
935    balance_deep_and_shallow(wc, zc, hc, wv, zv, hv, hvbar,
936                             xmomc, ymomc, xmomv, ymomv)
937
938
939
940
941###############################################
942#Boundaries - specific to the shallow water wave equation
943class Reflective_boundary(Boundary):
944    """Reflective boundary returns same conserved quantities as
945    those present in its neighbour volume but reflected.
946
947    This class is specific to the shallow water equation as it
948    works with the momentum quantities assumed to be the second
949    and third conserved quantities.
950    """
951
952    def __init__(self, domain = None):
953        Boundary.__init__(self)
954
955        if domain is None:
956            msg = 'Domain must be specified for reflective boundary'
957            raise msg
958
959        #Handy shorthands
960        self.stage   = domain.quantities['stage'].edge_values
961        self.xmom    = domain.quantities['xmomentum'].edge_values
962        self.ymom    = domain.quantities['ymomentum'].edge_values
963        self.normals = domain.normals
964
965        from Numeric import zeros, Float
966        self.conserved_quantities = zeros(3, Float)
967
968    def __repr__(self):
969        return 'Reflective_boundary'
970
971
972    def evaluate(self, vol_id, edge_id):
973        """Reflective boundaries reverses the outward momentum
974        of the volume they serve.
975        """
976
977        q = self.conserved_quantities
978        q[0] = self.stage[vol_id, edge_id]
979        q[1] = self.xmom[vol_id, edge_id]
980        q[2] = self.ymom[vol_id, edge_id]
981
982        normal = self.normals[vol_id, 2*edge_id:2*edge_id+2]
983
984
985        r = rotate(q, normal, direction = 1)
986        r[1] = -r[1]
987        q = rotate(r, normal, direction = -1)
988
989        return q
990
991
992
993class Transmissive_Momentum_Set_Stage_boundary(Boundary):
994    """Returns same momentum conserved quantities as
995    those present in its neighbour volume. Sets stage
996
997    Underlying domain must be specified when boundary is instantiated
998    """
999
1000    def __init__(self, domain = None, function=None):
1001        Boundary.__init__(self)
1002
1003        if domain is None:
1004            msg = 'Domain must be specified for this type boundary'
1005            raise msg
1006
1007        if function is None:
1008            msg = 'Function must be specified for this type boundary'
1009            raise msg
1010
1011        self.domain   = domain
1012        self.function = function
1013
1014    def __repr__(self):
1015        return 'Transmissive_Momentum_Set_Stage_boundary(%s)' %self.domain
1016
1017    def evaluate(self, vol_id, edge_id):
1018        """Transmissive Momentum Set Stage boundaries return the edge momentum
1019        values of the volume they serve.
1020        """
1021
1022        q = self.domain.get_conserved_quantities(vol_id, edge = edge_id)
1023        value = self.function(self.domain.time)
1024        q[0] = value[0]
1025        return q
1026
1027
1028        #FIXME: Consider this (taken from File_boundary) to allow
1029        #spatial variation
1030        #if vol_id is not None and edge_id is not None:
1031        #    i = self.boundary_indices[ vol_id, edge_id ]
1032        #    return self.F(t, point_id = i)
1033        #else:
1034        #    return self.F(t)
1035
1036
1037
1038class Dirichlet_Discharge_boundary(Boundary):
1039    """
1040    Sets stage (stage0)
1041    Sets momentum (wh0) in the inward normal direction.
1042
1043    Underlying domain must be specified when boundary is instantiated
1044    """
1045
1046    def __init__(self, domain = None, stage0=None, wh0=None):
1047        Boundary.__init__(self)
1048
1049        if domain is None:
1050            msg = 'Domain must be specified for this type boundary'
1051            raise msg
1052
1053        if stage0 is None:
1054            raise 'set stage'
1055
1056        if wh0 is None:
1057            wh0 = 0.0
1058
1059        self.domain   = domain
1060        self.stage0  = stage0
1061        self.wh0 = wh0
1062
1063    def __repr__(self):
1064        return 'Dirichlet_Discharge_boundary(%s)' %self.domain
1065
1066    def evaluate(self, vol_id, edge_id):
1067        """Set discharge in the (inward) normal direction
1068        """
1069
1070        normal = self.domain.get_normal(vol_id,edge_id)
1071        q = [self.stage0, -self.wh0*normal[0], -self.wh0*normal[1]]
1072        return q
1073
1074
1075        #FIXME: Consider this (taken from File_boundary) to allow
1076        #spatial variation
1077        #if vol_id is not None and edge_id is not None:
1078        #    i = self.boundary_indices[ vol_id, edge_id ]
1079        #    return self.F(t, point_id = i)
1080        #else:
1081        #    return self.F(t)
1082
1083
1084
1085#class Spatio_temporal_boundary(Boundary):
1086#    """The spatio-temporal boundary, reads values for the conserved
1087#    quantities from an sww NetCDF file, and returns interpolated values
1088#    at the midpoints of each associated boundaty segment.
1089#    Time dependency is interpolated linearly as in util.File_function.#
1090#
1091#    Example:
1092#    Bf = Spatio_temporal_boundary('source_file.sww', domain)
1093#
1094#    """
1095Spatio_temporal_boundary = File_boundary
1096
1097
1098
1099
1100#########################
1101#Standard forcing terms:
1102#
1103def gravity(domain):
1104    """Apply gravitational pull in the presence of bed slope
1105    """
1106
1107    from Numeric import zeros, Float, array, sum
1108
1109    xmom = domain.quantities['xmomentum'].explicit_update
1110    ymom = domain.quantities['ymomentum'].explicit_update
1111
1112    Stage = domain.quantities['stage']
1113    Elevation = domain.quantities['elevation']
1114    h = Stage.edge_values - Elevation.edge_values
1115    v = Elevation.vertex_values
1116
1117    x = domain.get_vertex_coordinates()
1118    g = domain.g
1119
1120    for k in range(domain.number_of_elements):
1121        avg_h = sum( h[k,:] )/3
1122
1123        #Compute bed slope
1124        x0, y0, x1, y1, x2, y2 = x[k,:]
1125        z0, z1, z2 = v[k,:]
1126
1127        zx, zy = gradient(x0, y0, x1, y1, x2, y2, z0, z1, z2)
1128
1129        #Update momentum
1130        xmom[k] += -g*zx*avg_h
1131        ymom[k] += -g*zy*avg_h
1132
1133
1134def gravity_c(domain):
1135    """Wrapper calling C version
1136    """
1137
1138    xmom = domain.quantities['xmomentum'].explicit_update
1139    ymom = domain.quantities['ymomentum'].explicit_update
1140
1141    Stage = domain.quantities['stage']
1142    Elevation = domain.quantities['elevation']
1143    h = Stage.edge_values - Elevation.edge_values
1144    v = Elevation.vertex_values
1145
1146    x = domain.get_vertex_coordinates()
1147    g = domain.g
1148
1149
1150    from shallow_water_ext import gravity
1151    gravity(g, h, v, x, xmom, ymom)
1152
1153
1154def manning_friction(domain):
1155    """Apply (Manning) friction to water momentum
1156    """
1157
1158    from math import sqrt
1159
1160    w = domain.quantities['stage'].centroid_values
1161    z = domain.quantities['elevation'].centroid_values
1162    h = w-z
1163
1164    uh = domain.quantities['xmomentum'].centroid_values
1165    vh = domain.quantities['ymomentum'].centroid_values
1166    eta = domain.quantities['friction'].centroid_values
1167
1168    xmom_update = domain.quantities['xmomentum'].semi_implicit_update
1169    ymom_update = domain.quantities['ymomentum'].semi_implicit_update
1170
1171    N = domain.number_of_elements
1172    eps = domain.minimum_allowed_height
1173    g = domain.g
1174
1175    for k in range(N):
1176        if eta[k] >= eps:
1177            if h[k] >= eps:
1178                S = -g * eta[k]**2 * sqrt((uh[k]**2 + vh[k]**2))
1179                S /= h[k]**(7.0/3)
1180
1181                #Update momentum
1182                xmom_update[k] += S*uh[k]
1183                ymom_update[k] += S*vh[k]
1184
1185
1186def manning_friction_implicit_c(domain):
1187    """Wrapper for c version
1188    """
1189
1190
1191    xmom = domain.quantities['xmomentum']
1192    ymom = domain.quantities['ymomentum']
1193
1194    w = domain.quantities['stage'].centroid_values
1195    z = domain.quantities['elevation'].centroid_values
1196
1197    uh = xmom.centroid_values
1198    vh = ymom.centroid_values
1199    eta = domain.quantities['friction'].centroid_values
1200
1201    xmom_update = xmom.semi_implicit_update
1202    ymom_update = ymom.semi_implicit_update
1203
1204    N = domain.number_of_elements
1205    eps = domain.minimum_allowed_height
1206    g = domain.g
1207
1208    from shallow_water_ext import manning_friction
1209    manning_friction(g, eps, w, z, uh, vh, eta, xmom_update, ymom_update)
1210
1211
1212def manning_friction_explicit_c(domain):
1213    """Wrapper for c version
1214    """
1215
1216    xmom = domain.quantities['xmomentum']
1217    ymom = domain.quantities['ymomentum']
1218
1219    w = domain.quantities['stage'].centroid_values
1220    z = domain.quantities['elevation'].centroid_values
1221
1222    uh = xmom.centroid_values
1223    vh = ymom.centroid_values
1224    eta = domain.quantities['friction'].centroid_values
1225
1226    xmom_update = xmom.explicit_update
1227    ymom_update = ymom.explicit_update
1228
1229    N = domain.number_of_elements
1230    eps = domain.minimum_allowed_height
1231    g = domain.g
1232
1233    from shallow_water_ext import manning_friction
1234    manning_friction(g, eps, w, z, uh, vh, eta, xmom_update, ymom_update)
1235
1236
1237def linear_friction(domain):
1238    """Apply linear friction to water momentum
1239
1240    Assumes quantity: 'linear_friction' to be present
1241    """
1242
1243    from math import sqrt
1244
1245    w = domain.quantities['stage'].centroid_values
1246    z = domain.quantities['elevation'].centroid_values
1247    h = w-z
1248
1249    uh = domain.quantities['xmomentum'].centroid_values
1250    vh = domain.quantities['ymomentum'].centroid_values
1251    tau = domain.quantities['linear_friction'].centroid_values
1252
1253    xmom_update = domain.quantities['xmomentum'].semi_implicit_update
1254    ymom_update = domain.quantities['ymomentum'].semi_implicit_update
1255
1256    N = domain.number_of_elements
1257    eps = domain.minimum_allowed_height
1258    g = domain.g #Not necessary? Why was this added?
1259
1260    for k in range(N):
1261        if tau[k] >= eps:
1262            if h[k] >= eps:
1263                S = -tau[k]/h[k]
1264
1265                #Update momentum
1266                xmom_update[k] += S*uh[k]
1267                ymom_update[k] += S*vh[k]
1268
1269
1270
1271def check_forcefield(f):
1272    """Check that f is either
1273    1: a callable object f(t,x,y), where x and y are vectors
1274       and that it returns an array or a list of same length
1275       as x and y
1276    2: a scalar
1277    """
1278
1279    from Numeric import ones, Float, array
1280
1281
1282    if callable(f):
1283        N = 3
1284        x = ones(3, Float)
1285        y = ones(3, Float)
1286        try:
1287            q = f(1.0, x=x, y=y)
1288        except Exception, e:
1289            msg = 'Function %s could not be executed:\n%s' %(f, e)
1290            #FIXME: Reconsider this semantics
1291            raise msg
1292
1293        try:
1294            q = array(q).astype(Float)
1295        except:
1296            msg = 'Return value from vector function %s could ' %f
1297            msg += 'not be converted into a Numeric array of floats.\n'
1298            msg += 'Specified function should return either list or array.'
1299            raise msg
1300
1301        #Is this really what we want?
1302        msg = 'Return vector from function %s ' %f
1303        msg += 'must have same lenght as input vectors'
1304        assert len(q) == N, msg
1305
1306    else:
1307        try:
1308            f = float(f)
1309        except:
1310            msg = 'Force field %s must be either a scalar' %f
1311            msg += ' or a vector function'
1312            raise msg
1313    return f
1314
1315
1316class Wind_stress:
1317    """Apply wind stress to water momentum in terms of
1318    wind speed [m/s] and wind direction [degrees]
1319    """
1320
1321    def __init__(self, *args, **kwargs):
1322        """Initialise windfield from wind speed s [m/s]
1323        and wind direction phi [degrees]
1324
1325        Inputs v and phi can be either scalars or Python functions, e.g.
1326
1327        W = Wind_stress(10, 178)
1328
1329        #FIXME - 'normal' degrees are assumed for now, i.e. the
1330        vector (1,0) has zero degrees.
1331        We may need to convert from 'compass' degrees later on and also
1332        map from True north to grid north.
1333
1334        Arguments can also be Python functions of t,x,y as in
1335
1336        def speed(t,x,y):
1337            ...
1338            return s
1339
1340        def angle(t,x,y):
1341            ...
1342            return phi
1343
1344        where x and y are vectors.
1345
1346        and then pass the functions in
1347
1348        W = Wind_stress(speed, angle)
1349
1350        The instantiated object W can be appended to the list of
1351        forcing_terms as in
1352
1353        Alternatively, one vector valued function for (speed, angle)
1354        can be applied, providing both quantities simultaneously.
1355        As in
1356        W = Wind_stress(F), where returns (speed, angle) for each t.
1357
1358        domain.forcing_terms.append(W)
1359        """
1360
1361        from config import rho_a, rho_w, eta_w
1362        from Numeric import array, Float
1363
1364        if len(args) == 2:
1365            s = args[0]
1366            phi = args[1]
1367        elif len(args) == 1:
1368            #Assume vector function returning (s, phi)(t,x,y)
1369            vector_function = args[0]
1370            s = lambda t,x,y: vector_function(t,x=x,y=y)[0]
1371            phi = lambda t,x,y: vector_function(t,x=x,y=y)[1]
1372        else:
1373           #Assume info is in 2 keyword arguments
1374
1375           if len(kwargs) == 2:
1376               s = kwargs['s']
1377               phi = kwargs['phi']
1378           else:
1379               raise 'Assumes two keyword arguments: s=..., phi=....'
1380
1381        self.speed = check_forcefield(s)
1382        self.phi = check_forcefield(phi)
1383
1384        self.const = eta_w*rho_a/rho_w
1385
1386
1387    def __call__(self, domain):
1388        """Evaluate windfield based on values found in domain
1389        """
1390
1391        from math import pi, cos, sin, sqrt
1392        from Numeric import Float, ones, ArrayType
1393
1394        xmom_update = domain.quantities['xmomentum'].explicit_update
1395        ymom_update = domain.quantities['ymomentum'].explicit_update
1396
1397        N = domain.number_of_elements
1398        t = domain.time
1399
1400        if callable(self.speed):
1401            xc = domain.get_centroid_coordinates()
1402            s_vec = self.speed(t, xc[:,0], xc[:,1])
1403        else:
1404            #Assume s is a scalar
1405
1406            try:
1407                s_vec = self.speed * ones(N, Float)
1408            except:
1409                msg = 'Speed must be either callable or a scalar: %s' %self.s
1410                raise msg
1411
1412
1413        if callable(self.phi):
1414            xc = domain.get_centroid_coordinates()
1415            phi_vec = self.phi(t, xc[:,0], xc[:,1])
1416        else:
1417            #Assume phi is a scalar
1418
1419            try:
1420                phi_vec = self.phi * ones(N, Float)
1421            except:
1422                msg = 'Angle must be either callable or a scalar: %s' %self.phi
1423                raise msg
1424
1425        assign_windfield_values(xmom_update, ymom_update,
1426                                s_vec, phi_vec, self.const)
1427
1428
1429def assign_windfield_values(xmom_update, ymom_update,
1430                            s_vec, phi_vec, const):
1431    """Python version of assigning wind field to update vectors.
1432    A c version also exists (for speed)
1433    """
1434    from math import pi, cos, sin, sqrt
1435
1436    N = len(s_vec)
1437    for k in range(N):
1438        s = s_vec[k]
1439        phi = phi_vec[k]
1440
1441        #Convert to radians
1442        phi = phi*pi/180
1443
1444        #Compute velocity vector (u, v)
1445        u = s*cos(phi)
1446        v = s*sin(phi)
1447
1448        #Compute wind stress
1449        S = const * sqrt(u**2 + v**2)
1450        xmom_update[k] += S*u
1451        ymom_update[k] += S*v
1452
1453
1454
1455##############################
1456#OBSOLETE STUFF
1457
1458def balance_deep_and_shallow_old(domain):
1459    """Compute linear combination between stage as computed by
1460    gradient-limiters and stage computed as constant height above bed.
1461    The former takes precedence when heights are large compared to the
1462    bed slope while the latter takes precedence when heights are
1463    relatively small.  Anything in between is computed as a balanced
1464    linear combination in order to avoid numerical disturbances which
1465    would otherwise appear as a result of hard switching between
1466    modes.
1467    """
1468
1469    #OBSOLETE
1470
1471    #Shortcuts
1472    wc = domain.quantities['stage'].centroid_values
1473    zc = domain.quantities['elevation'].centroid_values
1474    hc = wc - zc
1475
1476    wv = domain.quantities['stage'].vertex_values
1477    zv = domain.quantities['elevation'].vertex_values
1478    hv = wv-zv
1479
1480
1481    #Computed linear combination between constant stages and and
1482    #stages parallel to the bed elevation.
1483    for k in range(domain.number_of_elements):
1484        #Compute maximal variation in bed elevation
1485        #  This quantitiy is
1486        #    dz = max_i abs(z_i - z_c)
1487        #  and it is independent of dimension
1488        #  In the 1d case zc = (z0+z1)/2
1489        #  In the 2d case zc = (z0+z1+z2)/3
1490
1491        dz = max(abs(zv[k,0]-zc[k]),
1492                 abs(zv[k,1]-zc[k]),
1493                 abs(zv[k,2]-zc[k]))
1494
1495
1496        hmin = min( hv[k,:] )
1497
1498        #Create alpha in [0,1], where alpha==0 means using shallow
1499        #first order scheme and alpha==1 means using the stage w as
1500        #computed by the gradient limiter (1st or 2nd order)
1501        #
1502        #If hmin > dz/2 then alpha = 1 and the bed will have no effect
1503        #If hmin < 0 then alpha = 0 reverting to constant height above bed.
1504
1505        if dz > 0.0:
1506            alpha = max( min( 2*hmin/dz, 1.0), 0.0 )
1507        else:
1508            #Flat bed
1509            alpha = 1.0
1510
1511        #Weighted balance between stage parallel to bed elevation
1512        #(wvi = zvi + hc) and stage as computed by 1st or 2nd
1513        #order gradient limiter
1514        #(wvi = zvi + hvi) where i=0,1,2 denotes the vertex ids
1515        #
1516        #It follows that the updated wvi is
1517        #  wvi := (1-alpha)*(zvi+hc) + alpha*(zvi+hvi) =
1518        #  zvi + hc + alpha*(hvi - hc)
1519        #
1520        #Note that hvi = zc+hc-zvi in the first order case (constant).
1521
1522        if alpha < 1:
1523            for i in range(3):
1524                wv[k,i] = zv[k,i] + hc[k] + alpha*(hv[k,i]-hc[k])
1525
1526
1527            #Momentums at centroids
1528            xmomc = domain.quantities['xmomentum'].centroid_values
1529            ymomc = domain.quantities['ymomentum'].centroid_values
1530
1531            #Momentums at vertices
1532            xmomv = domain.quantities['xmomentum'].vertex_values
1533            ymomv = domain.quantities['ymomentum'].vertex_values
1534
1535            # Update momentum as a linear combination of
1536            # xmomc and ymomc (shallow) and momentum
1537            # from extrapolator xmomv and ymomv (deep).
1538            xmomv[k,:] = (1-alpha)*xmomc[k] + alpha*xmomv[k,:]
1539            ymomv[k,:] = (1-alpha)*ymomc[k] + alpha*ymomv[k,:]
1540
1541
1542
1543
1544
1545###########################
1546###########################
1547#Geometries
1548
1549
1550#FIXME: Rethink this way of creating values.
1551
1552
1553class Weir:
1554    """Set a bathymetry for weir with a hole and a downstream gutter
1555    x,y are assumed to be in the unit square
1556    """
1557
1558    def __init__(self, stage):
1559        self.inflow_stage = stage
1560
1561    def __call__(self, x, y):
1562        from Numeric import zeros, Float
1563        from math import sqrt
1564
1565        N = len(x)
1566        assert N == len(y)
1567
1568        z = zeros(N, Float)
1569        for i in range(N):
1570            z[i] = -x[i]/2  #General slope
1571
1572            #Flattish bit to the left
1573            if x[i] < 0.3:
1574                z[i] = -x[i]/10
1575
1576            #Weir
1577            if x[i] >= 0.3 and x[i] < 0.4:
1578                z[i] = -x[i]+0.9
1579
1580            #Dip
1581            x0 = 0.6
1582            #depth = -1.3
1583            depth = -1.0
1584            #plateaux = -0.9
1585            plateaux = -0.6
1586            if y[i] < 0.7:
1587                if x[i] > x0 and x[i] < 0.9:
1588                    z[i] = depth
1589
1590                #RHS plateaux
1591                if x[i] >= 0.9:
1592                    z[i] = plateaux
1593
1594
1595            elif y[i] >= 0.7 and y[i] < 1.5:
1596                #Restrict and deepen
1597                if x[i] >= x0 and x[i] < 0.8:
1598                    z[i] = depth-(y[i]/3-0.3)
1599                    #z[i] = depth-y[i]/5
1600                    #z[i] = depth
1601                elif x[i] >= 0.8:
1602                    #RHS plateaux
1603                    z[i] = plateaux
1604
1605            elif y[i] >= 1.5:
1606                if x[i] >= x0 and x[i] < 0.8 + (y[i]-1.5)/1.2:
1607                    #Widen up and stay at constant depth
1608                    z[i] = depth-1.5/5
1609                elif x[i] >= 0.8 + (y[i]-1.5)/1.2:
1610                    #RHS plateaux
1611                    z[i] = plateaux
1612
1613
1614            #Hole in weir (slightly higher than inflow condition)
1615            if x[i] >= 0.3 and x[i] < 0.4 and y[i] > 0.2 and y[i] < 0.4:
1616                z[i] = -x[i]+self.inflow_stage + 0.02
1617
1618            #Channel behind weir
1619            x0 = 0.5
1620            if x[i] >= 0.4 and x[i] < x0 and y[i] > 0.2 and y[i] < 0.4:
1621                z[i] = -x[i]+self.inflow_stage + 0.02
1622
1623            if x[i] >= x0 and x[i] < 0.6 and y[i] > 0.2 and y[i] < 0.4:
1624                #Flatten it out towards the end
1625                z[i] = -x0+self.inflow_stage + 0.02 + (x0-x[i])/5
1626
1627            #Hole to the east
1628            x0 = 1.1; y0 = 0.35
1629            #if x[i] < -0.2 and y < 0.5:
1630            if sqrt((2*(x[i]-x0))**2 + (2*(y[i]-y0))**2) < 0.2:
1631                z[i] = sqrt(((x[i]-x0))**2 + ((y[i]-y0))**2)-1.0
1632
1633            #Tiny channel draining hole
1634            if x[i] >= 1.14 and x[i] < 1.2 and y[i] >= 0.4 and y[i] < 0.6:
1635                z[i] = -0.9 #North south
1636
1637            if x[i] >= 0.9 and x[i] < 1.18 and y[i] >= 0.58 and y[i] < 0.65:
1638                z[i] = -1.0 + (x[i]-0.9)/3 #East west
1639
1640
1641
1642            #Stuff not in use
1643
1644            #Upward slope at inlet to the north west
1645            #if x[i] < 0.0: # and y[i] > 0.5:
1646            #    #z[i] = -y[i]+0.5  #-x[i]/2
1647            #    z[i] = x[i]/4 - y[i]**2 + 0.5
1648
1649            #Hole to the west
1650            #x0 = -0.4; y0 = 0.35 # center
1651            #if sqrt((2*(x[i]-x0))**2 + (2*(y[i]-y0))**2) < 0.2:
1652            #    z[i] = sqrt(((x[i]-x0))**2 + ((y[i]-y0))**2)-0.2
1653
1654
1655
1656
1657
1658        return z/2
1659
1660class Weir_simple:
1661    """Set a bathymetry for weir with a hole and a downstream gutter
1662    x,y are assumed to be in the unit square
1663    """
1664
1665    def __init__(self, stage):
1666        self.inflow_stage = stage
1667
1668    def __call__(self, x, y):
1669        from Numeric import zeros, Float
1670
1671        N = len(x)
1672        assert N == len(y)
1673
1674        z = zeros(N, Float)
1675        for i in range(N):
1676            z[i] = -x[i]  #General slope
1677
1678            #Flat bit to the left
1679            if x[i] < 0.3:
1680                z[i] = -x[i]/10  #General slope
1681
1682            #Weir
1683            if x[i] > 0.3 and x[i] < 0.4:
1684                z[i] = -x[i]+0.9
1685
1686            #Dip
1687            if x[i] > 0.6 and x[i] < 0.9:
1688                z[i] = -x[i]-0.5  #-y[i]/5
1689
1690            #Hole in weir (slightly higher than inflow condition)
1691            if x[i] > 0.3 and x[i] < 0.4 and y[i] > 0.2 and y[i] < 0.4:
1692                z[i] = -x[i]+self.inflow_stage + 0.05
1693
1694
1695        return z/2
1696
1697
1698
1699class Constant_stage:
1700    """Set an initial condition with constant stage
1701    """
1702    def __init__(self, s):
1703        self.s = s
1704
1705    def __call__(self, x, y):
1706        return self.s
1707
1708class Constant_height:
1709    """Set an initial condition with constant water height, e.g
1710    stage s = z+h
1711    """
1712
1713    def __init__(self, W, h):
1714        self.W = W
1715        self.h = h
1716
1717    def __call__(self, x, y):
1718        if self.W is None:
1719            from Numeric import ones, Float
1720            return self.h*ones(len(x), Float)
1721        else:
1722            return self.W(x,y) + self.h
1723
1724
1725
1726
1727def compute_fluxes_python(domain):
1728    """Compute all fluxes and the timestep suitable for all volumes
1729    in domain.
1730
1731    Compute total flux for each conserved quantity using "flux_function"
1732
1733    Fluxes across each edge are scaled by edgelengths and summed up
1734    Resulting flux is then scaled by area and stored in
1735    explicit_update for each of the three conserved quantities
1736    stage, xmomentum and ymomentum
1737
1738    The maximal allowable speed computed by the flux_function for each volume
1739    is converted to a timestep that must not be exceeded. The minimum of
1740    those is computed as the next overall timestep.
1741
1742    Post conditions:
1743      domain.explicit_update is reset to computed flux values
1744      domain.timestep is set to the largest step satisfying all volumes.
1745    """
1746
1747    import sys
1748    from Numeric import zeros, Float
1749
1750    N = domain.number_of_elements
1751
1752    #Shortcuts
1753    Stage = domain.quantities['stage']
1754    Xmom = domain.quantities['xmomentum']
1755    Ymom = domain.quantities['ymomentum']
1756    Bed = domain.quantities['elevation']
1757
1758    #Arrays
1759    stage = Stage.edge_values
1760    xmom =  Xmom.edge_values
1761    ymom =  Ymom.edge_values
1762    bed =   Bed.edge_values
1763
1764    stage_bdry = Stage.boundary_values
1765    xmom_bdry =  Xmom.boundary_values
1766    ymom_bdry =  Ymom.boundary_values
1767
1768    flux = zeros((N,3), Float) #Work array for summing up fluxes
1769
1770    #Loop
1771    timestep = float(sys.maxint)
1772    for k in range(N):
1773
1774        for i in range(3):
1775            #Quantities inside volume facing neighbour i
1776            ql = [stage[k, i], xmom[k, i], ymom[k, i]]
1777            zl = bed[k, i]
1778
1779            #Quantities at neighbour on nearest face
1780            n = domain.neighbours[k,i]
1781            if n < 0:
1782                m = -n-1 #Convert negative flag to index
1783                qr = [stage_bdry[m], xmom_bdry[m], ymom_bdry[m]]
1784                zr = zl #Extend bed elevation to boundary
1785            else:
1786                m = domain.neighbour_edges[k,i]
1787                qr = [stage[n, m], xmom[n, m], ymom[n, m]]
1788                zr = bed[n, m]
1789
1790
1791            #Outward pointing normal vector
1792            normal = domain.normals[k, 2*i:2*i+2]
1793
1794            #Flux computation using provided function
1795            edgeflux, max_speed = flux_function(normal, ql, qr, zl, zr)
1796
1797            flux[k,:] = edgeflux
1798
1799    return flux
1800
1801
1802
1803
1804
1805
1806
1807##############################################
1808#Initialise module
1809
1810
1811from utilities import compile
1812if compile.can_use_C_extension('shallow_water_ext.c'):
1813    #Replace python version with c implementations
1814
1815    from shallow_water_ext import rotate, assign_windfield_values
1816    compute_fluxes = compute_fluxes_c
1817    extrapolate_second_order_sw=extrapolate_second_order_sw_c
1818    gravity = gravity_c
1819    manning_friction = manning_friction_implicit_c
1820    h_limiter = h_limiter_c
1821    balance_deep_and_shallow = balance_deep_and_shallow_c
1822    protect_against_infinitesimal_and_negative_heights = protect_against_infinitesimal_and_negative_heights_c
1823
1824
1825    #distribute_to_vertices_and_edges = distribute_to_vertices_and_edges_c #(like MH's)
1826
1827
1828
1829#Optimisation with psyco
1830from config import use_psyco
1831if use_psyco:
1832    try:
1833        import psyco
1834    except:
1835        import os
1836        if os.name == 'posix' and os.uname()[4] == 'x86_64':
1837            pass
1838            #Psyco isn't supported on 64 bit systems, but it doesn't matter
1839        else:
1840            msg = 'WARNING: psyco (speedup) could not import'+\
1841                  ', you may want to consider installing it'
1842            print msg
1843    else:
1844        psyco.bind(Domain.distribute_to_vertices_and_edges)
1845        psyco.bind(Domain.compute_fluxes)
1846
1847if __name__ == "__main__":
1848    pass
1849
1850# Profiling stuff
1851import profile
1852profiler = profile.Profile()
Note: See TracBrowser for help on using the repository browser.