source: inundation/pyvolution/shallow_water.py @ 1797

Last change on this file since 1797 was 1704, checked in by ole, 20 years ago

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