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

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