source: storm_surge/pyvolution/shallow_water.py @ 1

Last change on this file since 1 was 1, checked in by duncan, 20 years ago

initial import

File size: 29.5 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
8FIXME: Write equations here!
9
10
11Conserved quantities are w (water level or stage), uh (x momentum)
12and vh (y momentum).
13
14
15Ole Nielsen, Stephen Roberts, Duncan Gray, Christopher Zoppou
16Geoscience Australia, 2004   
17"""
18
19from domain import *
20Generic_domain = Domain #Rename
21
22class Domain(Generic_domain):
23
24    def __init__(self, coordinates, vertices, boundary = None):
25
26        conserved_quantities = ['level', 'xmomentum', 'ymomentum']
27        other_quantities = ['elevation', 'friction'] 
28       
29        Generic_domain.__init__(self, coordinates, vertices, boundary,
30                                conserved_quantities, other_quantities)
31
32        from config import minimum_allowed_height, g
33        self.minimum_allowed_height = minimum_allowed_height
34        self.g = g
35
36        self.forcing_terms.append(gravity)
37        self.forcing_terms.append(manning_friction)
38
39        #Realtime visualisation
40        self.visualise = False
41
42        #Stored output
43        self.store=False
44        self.format = 'sww'   
45        self.smooth = True
46
47        #Reduction operation for get_vertex_values             
48        #from pytools.stats import mean
49        #self.reduction = mean
50        self.reduction = min  #Looks better near steep slopes
51       
52
53        #Establish shortcuts to relevant quantities (for efficiency)
54        self.w = self.quantities['level']
55        self.uh = self.quantities['xmomentum']                 
56        self.vh = self.quantities['ymomentum']                         
57        self.z = self.quantities['elevation']           
58        self.eta = self.quantities['friction']                 
59
60    def check_integrity(self):
61        Generic_domain.check_integrity(self)
62
63        #Check that we are solving the shallow water wave equation
64
65        msg = 'First conserved quantity must be "level"'
66        assert self.conserved_quantities[0] == 'level', msg
67        msg = 'Second conserved quantity must be "xmomentum"'
68        assert self.conserved_quantities[1] == 'xmomentum', msg
69        msg = 'Third conserved quantity must be "ymomentum"'
70        assert self.conserved_quantities[2] == 'ymomentum', msg
71       
72
73        #Check that levels are >= bed elevation
74        from Numeric import alltrue, greater_equal
75       
76        level = self.quantities['level']
77        bed = self.quantities['elevation']       
78
79        msg = 'All water levels must be greater than the bed elevation'
80        assert alltrue( greater_equal(
81            level.vertex_values, bed.vertex_values )), msg
82
83        assert alltrue( greater_equal(
84            level.edge_values, bed.edge_values )), msg
85
86        assert alltrue( greater_equal(
87            level.centroid_values, bed.centroid_values )), msg       
88
89
90    def compute_fluxes(self):
91        #Call correct module function
92        #(either from this module or C-extension)
93        compute_fluxes(self)
94
95    def distribute_to_vertices_and_edges(self):
96        #Call correct module function
97        #(either from this module or C-extension)       
98        distribute_to_vertices_and_edges(self)
99
100
101    def evolve(self, yieldstep = None, finaltime = None):
102        """Specialisation of basic evolve method from parent class
103        """
104       
105        #Initialise real time viz if requested
106        if self.visualise is True and self.time == 0.0:
107            import realtime_visualisation as visualise
108            visualise.create_surface(self)
109
110        #Store model data, e.g. for visualisation   
111        if self.store is True and self.time == 0.0:
112            self.initialise_storage()   
113           
114
115        #Call basic machinery from parent class
116        for t in Generic_domain.evolve(self, yieldstep, finaltime):
117            #Real time viz
118            if self.visualise is True:
119                visualise.update(self)
120
121            #Store model data, e.g. for subsequent visualisation   
122            if self.store is True:
123                self.store_timestep('level')
124                #FIXME: Could maybe be taken from specified list
125                #of 'store every step' quantities       
126
127
128            #Pass control on to outer loop for more specific actions   
129            yield(t)
130       
131
132    def initialise_storage(self):       
133        """Create and initialise self.writer object for storing data.
134        Also, save x,y and bed elevation
135        """
136
137        import data_manager
138       
139        #Initialise writer
140        self.writer = data_manager.get_dataobject(self, mode = 'w')
141
142        #Store vertices and connectivity
143        self.writer.store_connectivity()                   
144
145    def store_timestep(self, name):
146        """Store named quantity and time.
147
148        Precondition:
149           self.write has been initialised
150        """
151        self.writer.store_timestep(name)
152       
153#Rotation of momentum vector
154def rotate(q, normal, direction = 1):
155    """Rotate the momentum component q (q[1], q[2])
156    from x,y coordinates to coordinates based on normal vector.
157
158    If direction is negative the rotation is inverted.
159   
160    Input vector is preserved
161
162    This function is specific to the shallow water wave equation
163    """
164
165    #FIXME: Needs to be tested
166
167    from Numeric import zeros, Float
168   
169    assert len(q) == 3,\
170           'Vector of conserved quantities must have length 3'\
171           'for 2D shallow water equation'
172
173    try:
174        l = len(normal)
175    except:
176        raise 'Normal vector must be an Numeric array'
177
178    #FIXME: Put this test into C-extension as well
179    assert l == 2, 'Normal vector must have 2 components'
180
181 
182    n1 = normal[0]
183    n2 = normal[1]   
184   
185    r = zeros(len(q), Float) #Rotated quantities
186    r[0] = q[0]              #First quantity, height, is not rotated
187
188    if direction == -1:
189        n2 = -n2
190
191
192    r[1] =  n1*q[1] + n2*q[2]
193    r[2] = -n2*q[1] + n1*q[2]
194   
195    return r
196
197
198
199####################################
200# Flux computation       
201def flux_function(normal, ql, qr, zl, zr): 
202    """Compute fluxes between volumes for the shallow water wave equation
203    cast in terms of w = h+z using the 'central scheme' as described in
204
205    Kurganov, Noelle, Petrova. 'Semidiscrete Central-Upwind Schemes For
206    Hyperbolic Conservation Laws and Hamilton-Jacobi Equations'.
207    Siam J. Sci. Comput. Vol. 23, No. 3, pp. 707-740.
208
209    The implemented formula is given in equation (3.15) on page 714
210
211    Conserved quantities w, uh, vh are stored as elements 0, 1 and 2
212    in the numerical vectors ql an qr.
213
214    Bed elevations zl and zr.
215    """
216
217    from config import g, epsilon
218    from math import sqrt
219    from Numeric import array
220
221    #Align momentums with x-axis
222    q_left  = rotate(ql, normal, direction = 1)
223    q_right = rotate(qr, normal, direction = 1)
224
225    z = (zl+zr)/2 #Take average of field values
226
227    w_left  = q_left[0]   #w=h+z
228    h_left  = w_left-z
229    uh_left = q_left[1]
230
231    if h_left < epsilon:
232        u_left = 0.0  #Could have been negative
233        h_left = 0.0
234    else:   
235        u_left  = uh_left/h_left
236
237
238    w_right  = q_right[0]  #w=h+z
239    h_right  = w_right-z
240    uh_right = q_right[1]
241
242
243    if h_right < epsilon:
244        u_right = 0.0  #Could have been negative
245        h_right = 0.0
246    else:   
247        u_right  = uh_right/h_right
248
249    vh_left  = q_left[2]
250    vh_right = q_right[2]       
251
252    soundspeed_left  = sqrt(g*h_left) 
253    soundspeed_right = sqrt(g*h_right)
254
255    #Maximal wave speed
256    s_max = max(u_left+soundspeed_left, u_right+soundspeed_right, 0)
257   
258    #Minimal wave speed       
259    s_min = min(u_left-soundspeed_left, u_right-soundspeed_right, 0)
260
261    #Flux computation
262    flux_left  = array([u_left*h_left,
263                        u_left*uh_left + 0.5*g*h_left**2,
264                        u_left*vh_left])
265    flux_right = array([u_right*h_right,
266                        u_right*uh_right + 0.5*g*h_right**2,
267                        u_right*vh_right])   
268
269    denom = s_max-s_min
270    if denom == 0.0:
271        edgeflux = array([0.0, 0.0, 0.0])
272        max_speed = 0.0
273    else:   
274        edgeflux = (s_max*flux_left - s_min*flux_right)/denom
275        edgeflux += s_max*s_min*(q_right-q_left)/denom
276
277        edgeflux = rotate(edgeflux, normal, direction=-1)
278        max_speed = max(abs(s_max), abs(s_min))
279
280    return edgeflux, max_speed       
281
282
283def compute_fluxes(domain):
284    """Compute all fluxes and the timestep suitable for all volumes
285    in domain.
286
287    Compute total flux for each conserved quantity using "flux_function"
288
289    Fluxes across each edge are scaled by edgelengths and summed up
290    Resulting flux is then scaled by area and stored in
291    explicit_update for each of the three conserved quantities
292    level, xmomentum and ymomentum   
293
294    The maximal allowable speed computed by the flux_function for each volume
295    is converted to a timestep that must not be exceeded. The minimum of
296    those is computed as the next overall timestep.
297
298    Post conditions:
299      domain.explicit_update is reset to computed flux values
300      domain.timestep is set to the largest step satisfying all volumes.
301    """
302
303    import sys
304    from Numeric import zeros, Float
305
306    N = domain.number_of_elements
307   
308    #Shortcuts
309    Level = domain.quantities['level']
310    Xmom = domain.quantities['xmomentum']
311    Ymom = domain.quantities['ymomentum']
312    Bed = domain.quantities['elevation']       
313
314    #Arrays
315    level = Level.edge_values
316    xmom =  Xmom.edge_values
317    ymom =  Ymom.edge_values
318    bed =   Bed.edge_values   
319
320    level_bdry = Level.boundary_values
321    xmom_bdry =  Xmom.boundary_values
322    ymom_bdry =  Ymom.boundary_values
323   
324    flux = zeros(3, Float) #Work array for summing up fluxes
325
326    #Loop
327    timestep = float(sys.maxint)   
328    for k in range(N):
329
330        flux[:] = 0.  #Reset work array
331        for i in range(3):
332            #Quantities inside volume facing neighbour i
333            ql = [level[k, i], xmom[k, i], ymom[k, i]]
334            zl = bed[k, i]
335
336            #Quantities at neighbour on nearest face
337            n = domain.neighbours[k,i] 
338            if n < 0:
339                m = -n-1 #Convert negative flag to index
340                qr = [level_bdry[m], xmom_bdry[m], ymom_bdry[m]]
341                zr = zl #Extend bed elevation to boundary
342            else:   
343                m = domain.neighbour_edges[k,i]
344                qr = [level[n, m], xmom[n, m], ymom[n, m]]
345                zr = bed[n, m]               
346
347               
348            #Outward pointing normal vector   
349            normal = domain.normals[k, 2*i:2*i+2]
350
351            #Flux computation using provided function
352            edgeflux, max_speed = flux_function(normal, ql, qr, zl, zr)
353            flux -= edgeflux * domain.edgelengths[k,i]
354
355            #Update optimal_timestep
356            try:
357                timestep = min(timestep, domain.radii[k]/max_speed)
358            except ZeroDivisionError:
359                pass
360
361        #Normalise by area and store for when all conserved
362        #quantities get updated
363        flux /= domain.areas[k]
364        Level.explicit_update[k] = flux[0]
365        Xmom.explicit_update[k] = flux[1]
366        Ymom.explicit_update[k] = flux[2]
367
368    domain.timestep = timestep   
369
370
371def compute_fluxes_c(domain):
372    """Wrapper calling C version of compute fluxes
373    """
374
375    import sys
376    from Numeric import zeros, Float
377
378    N = domain.number_of_elements
379   
380    #Shortcuts
381    Level = domain.quantities['level']
382    Xmom = domain.quantities['xmomentum']
383    Ymom = domain.quantities['ymomentum']
384    Bed = domain.quantities['elevation']       
385
386    timestep = float(sys.maxint)   
387    from shallow_water_ext import compute_fluxes
388    domain.timestep = compute_fluxes(timestep, domain.epsilon, domain.g,
389                                     domain.neighbours,
390                                     domain.neighbour_edges,
391                                     domain.normals,
392                                     domain.edgelengths,                       
393                                     domain.radii,
394                                     domain.areas,
395                                     Level.edge_values, 
396                                     Xmom.edge_values, 
397                                     Ymom.edge_values, 
398                                     Bed.edge_values,   
399                                     Level.boundary_values,
400                                     Xmom.boundary_values,
401                                     Ymom.boundary_values,
402                                     Level.explicit_update,
403                                     Xmom.explicit_update,
404                                     Ymom.explicit_update)
405                                     
406
407####################################
408# Module functions for gradient limiting (distribute_to_vertices_and_edges)
409
410def distribute_to_vertices_and_edges(domain):
411    """Distribution from centroids to vertices specific to the
412    shallow water wave
413    equation.
414
415    It will ensure that h (w-z) is always non-negative even in the
416    presence of steep bed-slopes by taking a weighted average between shallow
417    and deep cases.
418   
419    In addition, all conserved quantities get distributed as per either a
420    constant (order==1) or a piecewise linear function (order==2).
421   
422    FIXME: more explanation about removal of artificial variability etc
423
424    Precondition:
425      All quantities defined at centroids and bed elevation defined at
426      vertices.
427
428    Postcondition
429      Conserved quantities defined at vertices
430   
431    """
432
433    #Remove very thin layers of water
434    protect_against_infintesimal_and_negative_heights(domain)
435
436    #Extrapolate all conserved quantities
437    for name in domain.conserved_quantities:
438        Q = domain.quantities[name]
439        if domain.order == 1:
440            Q.extrapolate_first_order()
441        elif domain.order == 2:
442            Q.extrapolate_second_order()
443            Q.limit()                               
444        else:
445            raise 'Unknown order'
446
447    #Take bed slevation into account when heights are small   
448    balance_deep_and_shallow(domain)
449   
450    #Compute edge values by interpolation
451    for name in domain.conserved_quantities:
452        Q = domain.quantities[name]
453        Q.interpolate_from_vertices_to_edges()
454                   
455   
456
457
458def protect_against_infintesimal_and_negative_heights(domain):
459    """Protect against infinitesimal heights and associated high velocities
460    """
461   
462    #Shortcuts
463    wc = domain.quantities['level'].centroid_values
464    zc = domain.quantities['elevation'].centroid_values
465    xmomc = domain.quantities['xmomentum'].centroid_values
466    ymomc = domain.quantities['ymomentum'].centroid_values           
467    hc = wc - zc  #Water depths at centroids   
468
469    #Update
470    for k in range(domain.number_of_elements):
471        if hc[k] < domain.minimum_allowed_height:       
472            if hc[k] < 0.0:
473                #Control level and height
474                wc[k] = zc[k]
475
476            #Control momentum
477            xmomc[k] = ymomc[k] = 0.0
478           
479
480
481def protect_against_infintesimal_and_negative_heights_c(domain):
482    """Protect against infinitesimal heights and associated high velocities
483    """
484   
485    #Shortcuts
486    wc = domain.quantities['level'].centroid_values
487    zc = domain.quantities['elevation'].centroid_values
488    xmomc = domain.quantities['xmomentum'].centroid_values
489    ymomc = domain.quantities['ymomentum'].centroid_values           
490
491    from shallow_water_ext import protect
492
493    protect(domain.minimum_allowed_height, wc, zc, xmomc, ymomc)
494   
495
496def balance_deep_and_shallow(domain):
497    """Compute linear combination between stage as computed by
498    gradient-limiters and stage computed as constant height above bed.
499    The former takes precedence when heights are large compared to the
500    bed slope while the latter takes precedence when heights are
501    relatively small.  Anything in between is computed as a balanced
502    linear combination in order to avoid numerical disturbances which
503    would otherwise appear as a result of hard switching between
504    modes.
505    """
506   
507    #Shortcuts
508    wc = domain.quantities['level'].centroid_values
509    zc = domain.quantities['elevation'].centroid_values
510    hc = wc - zc
511   
512    wv = domain.quantities['level'].vertex_values
513    zv = domain.quantities['elevation'].vertex_values
514    hv = wv-zv
515
516
517    #Computed linear combination between constant levels and and
518    #levels parallel to the bed elevation.     
519    for k in range(domain.number_of_elements):
520        #Compute maximal variation in bed elevation
521        #  This quantitiy is
522        #    dz = max_i abs(z_i - z_c)
523        #  and it is independent of dimension
524        #  In the 1d case zc = (z0+z1)/2
525        #  In the 2d case zc = (z0+z1+z2)/3
526
527        dz = max(abs(zv[k,0]-zc[k]),
528                 abs(zv[k,1]-zc[k]),
529                 abs(zv[k,2]-zc[k]))
530
531       
532        hmin = min( hv[k,:] )
533
534       
535        #Create alpha in [0,1], where alpha==0 means using shallow
536        #first order scheme and alpha==1 means using the stage w as
537        #computed by the gradient limiter (1st or 2nd order)
538        #
539        #If hmin > dz/2 then alpha = 1 and the bed will have no effect
540        #If hmin < 0 then alpha = 0 reverting to constant height above bed.
541     
542        if dz > 0.0:
543            alpha = max( min( 2*hmin/dz, 1.0), 0.0 )
544        else:
545            #Flat bed
546            alpha = 1.0
547
548
549        #Weighted balance between stage parallel to bed elevation
550        #(wvi = zvi + hc) and stage as computed by 1st or 2nd
551        #order gradient limiter
552        #(wvi = zvi + hvi) where i=0,1,2 denotes the vertex ids
553        #
554        #It follows that the updated wvi is
555        #  wvi := (1-alpha)*(zvi+hc) + alpha*(zvi+hvi) =
556        #  zvi + hc + alpha*(hvi - hc)
557        #
558        #Note that hvi = zc+hc-zvi in the first order case (constant).
559
560        if alpha < 1:         
561            for i in range(3):
562                wv[k,i] = zv[k,i] + hc[k] + alpha*(hv[k,i]-hc[k])
563           
564
565            #Momentums at centroids
566            xmomc = domain.quantities['xmomentum'].centroid_values
567            ymomc = domain.quantities['ymomentum'].centroid_values       
568
569            #Momentums at vertices
570            xmomv = domain.quantities['xmomentum'].vertex_values
571            ymomv = domain.quantities['ymomentum'].vertex_values         
572
573            # Update momentum as a linear combination of
574            # xmomc and ymomc (shallow) and momentum
575            # from extrapolator xmomv and ymomv (deep).
576            xmomv[k,:] = (1-alpha)*xmomc[k] + alpha*xmomv[k,:]           
577            ymomv[k,:] = (1-alpha)*ymomc[k] + alpha*ymomv[k,:]
578                                   
579           
580
581def balance_deep_and_shallow_c(domain):
582    """Wrapper for C implementation
583    """
584   
585    #Shortcuts
586    wc = domain.quantities['level'].centroid_values
587    zc = domain.quantities['elevation'].centroid_values
588    hc = wc - zc
589   
590    wv = domain.quantities['level'].vertex_values
591    zv = domain.quantities['elevation'].vertex_values
592    hv = wv-zv
593
594    #Momentums at centroids
595    xmomc = domain.quantities['xmomentum'].centroid_values
596    ymomc = domain.quantities['ymomentum'].centroid_values       
597
598    #Momentums at vertices
599    xmomv = domain.quantities['xmomentum'].vertex_values
600    ymomv = domain.quantities['ymomentum'].vertex_values         
601
602   
603
604    from shallow_water_ext import balance_deep_and_shallow
605    balance_deep_and_shallow(wc, zc, hc, wv, zv, hv,
606                             xmomc, ymomc, xmomv, ymomv)
607
608   
609
610
611###############################################
612#Boundary - specific to the shallow water wave equation
613class Reflective_boundary(Boundary):
614    """Reflective boundary returns same conserved quantities as
615    those present in its neighbour volume but reflected.
616
617    This class is specific to the shallow water equation as it
618    works with the momentum quantities assumed to be the second
619    and third conserved quantities.
620    """
621   
622    def __init__(self, domain = None):       
623        Boundary.__init__(self)
624
625        if domain is None:
626            msg = 'Domain must be specified for reflective boundary'
627            raise msg
628       
629        #Handy shorthands
630        self.level = domain.quantities['level'].edge_values
631        self.xmom = domain.quantities['xmomentum'].edge_values
632        self.ymom = domain.quantities['ymomentum'].edge_values         
633        self.normals = domain.normals
634       
635        from Numeric import zeros, Float       
636        self.conserved_quantities = zeros(3, Float) 
637       
638    def __repr__(self):
639        return 'Reflective_boundary'
640
641
642    def evaluate(self, vol_id, edge_id):
643        """Reflective boundaries reverses the outward momentum
644        of the volume they serve.
645        """
646
647        q = self.conserved_quantities
648        q[0] = self.level[vol_id, edge_id]
649        q[1] = self.xmom[vol_id, edge_id] 
650        q[2] = self.ymom[vol_id, edge_id] 
651       
652        normal = self.normals[vol_id, 2*edge_id:2*edge_id+2]     
653
654       
655        r = rotate(q, normal, direction = 1)
656        r[1] = -r[1]
657        q = rotate(r, normal, direction = -1)
658
659        return q
660
661
662#########################
663#Standard forcing terms:
664#
665def gravity(domain):
666    """Implement forcing function for bed slope working with
667    consecutive data structures of class Volume
668    """
669
670    from util import gradient
671    from Numeric import zeros, Float, array, sum
672
673    xmom = domain.quantities['xmomentum'].explicit_update
674    ymom = domain.quantities['ymomentum'].explicit_update
675
676    Level = domain.quantities['level']
677    Elevation = domain.quantities['elevation']   
678    h = Level.edge_values - Elevation.edge_values
679    v = Elevation.vertex_values
680   
681    x = domain.get_vertex_coordinates()
682    g = domain.g
683   
684    for k in range(domain.number_of_elements):   
685        avg_h = sum( h[k,:] )/3
686       
687        #Compute bed slope
688        x0, y0, x1, y1, x2, y2 = x[k,:]   
689        z0, z1, z2 = v[k,:]
690
691        zx, zy = gradient(x0, y0, x1, y1, x2, y2, z0, z1, z2)
692
693        #Update momentum
694        xmom[k] += -g*zx*avg_h
695        ymom[k] += -g*zy*avg_h       
696
697
698def gravity_c(domain):
699    """Wrapper calling C version
700    """   
701
702    xmom = domain.quantities['xmomentum'].explicit_update
703    ymom = domain.quantities['ymomentum'].explicit_update
704
705    Level = domain.quantities['level']   
706    Elevation = domain.quantities['elevation']   
707    h = Level.edge_values - Elevation.edge_values
708    v = Elevation.vertex_values
709   
710    x = domain.get_vertex_coordinates()
711    g = domain.g
712   
713   
714    from shallow_water_ext import gravity
715    gravity(g, h, v, x, xmom, ymom)
716   
717       
718def manning_friction(domain):
719    """Apply (Manning) friction to water momentum
720    """
721
722    from math import sqrt
723
724    w = domain.quantities['level'].centroid_values
725    uh = domain.quantities['xmomentum'].centroid_values
726    vh = domain.quantities['ymomentum'].centroid_values
727    eta = domain.quantities['friction'].centroid_values   
728   
729    xmom_update = domain.quantities['xmomentum'].semi_implicit_update
730    ymom_update = domain.quantities['ymomentum'].semi_implicit_update   
731   
732    N = domain.number_of_elements
733    eps = domain.minimum_allowed_height
734    g = domain.g
735   
736    for k in range(N):
737        if w[k] >= eps:
738            S = -g * eta[k]**2 * sqrt((uh[k]**2 + vh[k]**2))
739            S /= w[k]**(7.0/3)
740
741            #Update momentum
742            xmom_update[k] += S*uh[k]
743            ymom_update[k] += S*vh[k]
744           
745       
746def manning_friction_c(domain):
747    """Wrapper for c version
748    """
749
750
751    xmom = domain.quantities['xmomentum']
752    ymom = domain.quantities['ymomentum']   
753       
754    w = domain.quantities['level'].centroid_values
755    uh = xmom.centroid_values
756    vh = ymom.centroid_values
757    eta = domain.quantities['friction'].centroid_values   
758   
759    xmom_update = xmom.semi_implicit_update
760    ymom_update = ymom.semi_implicit_update   
761   
762    N = domain.number_of_elements
763    eps = domain.minimum_allowed_height
764    g = domain.g
765
766    from shallow_water_ext import manning_friction
767    manning_friction(g, eps, w, uh, vh, eta, xmom_update, ymom_update)
768       
769
770###########################
771###########################
772#Geometries
773
774
775#FIXME: Rethink this way of creating values.
776
777
778class Weir:
779    """Set a bathymetry for weir with a hole and a downstream gutter
780    x,y are assumed to be in the unit square
781    """
782   
783    def __init__(self, stage):
784        self.inflow_stage = stage
785
786    def __call__(self, x, y):   
787        from Numeric import zeros, Float
788        from math import sqrt
789   
790        N = len(x)
791        assert N == len(y)   
792
793        z = zeros(N, Float)
794        for i in range(N):
795            z[i] = -x[i]/2  #General slope
796
797            #Flattish bit to the left
798            if x[i] < 0.3:
799                z[i] = -x[i]/10 
800           
801            #Weir
802            if x[i] >= 0.3 and x[i] < 0.4:
803                z[i] = -x[i]+0.9
804
805            #Dip
806            x0 = 0.6
807            #depth = -1.3
808            depth = -1.0
809            #plateaux = -0.9
810            plateaux = -0.6           
811            if y[i] < 0.7:
812                if x[i] > x0 and x[i] < 0.9:
813                    z[i] = depth
814
815                #RHS plateaux
816                if x[i] >= 0.9:
817                    z[i] = plateaux
818                   
819                   
820            elif y[i] >= 0.7 and y[i] < 1.5:
821                #Restrict and deepen
822                if x[i] >= x0 and x[i] < 0.8:
823                    z[i] = depth-(y[i]/3-0.3)
824                    #z[i] = depth-y[i]/5
825                    #z[i] = depth
826                elif x[i] >= 0.8:   
827                    #RHS plateaux
828                    z[i] = plateaux
829                   
830            elif y[i] >= 1.5:
831                if x[i] >= x0 and x[i] < 0.8 + (y[i]-1.5)/1.2:
832                    #Widen up and stay at constant depth
833                    z[i] = depth-1.5/5
834                elif x[i] >= 0.8 + (y[i]-1.5)/1.2:                       
835                    #RHS plateaux
836                    z[i] = plateaux                                   
837                   
838
839            #Hole in weir (slightly higher than inflow condition)
840            if x[i] >= 0.3 and x[i] < 0.4 and y[i] > 0.2 and y[i] < 0.4:
841                z[i] = -x[i]+self.inflow_stage + 0.02
842
843            #Channel behind weir
844            x0 = 0.5
845            if x[i] >= 0.4 and x[i] < x0 and y[i] > 0.2 and y[i] < 0.4:
846                z[i] = -x[i]+self.inflow_stage + 0.02
847               
848            if x[i] >= x0 and x[i] < 0.6 and y[i] > 0.2 and y[i] < 0.4:
849                #Flatten it out towards the end
850                z[i] = -x0+self.inflow_stage + 0.02 + (x0-x[i])/5
851
852            #Hole to the east
853            x0 = 1.1; y0 = 0.35
854            #if x[i] < -0.2 and y < 0.5:
855            if sqrt((2*(x[i]-x0))**2 + (2*(y[i]-y0))**2) < 0.2:
856                z[i] = sqrt(((x[i]-x0))**2 + ((y[i]-y0))**2)-1.0
857
858            #Tiny channel draining hole
859            if x[i] >= 1.14 and x[i] < 1.2 and y[i] >= 0.4 and y[i] < 0.6:
860                z[i] = -0.9 #North south
861               
862            if x[i] >= 0.9 and x[i] < 1.18 and y[i] >= 0.58 and y[i] < 0.65:
863                z[i] = -1.0 + (x[i]-0.9)/3 #East west
864           
865           
866
867            #Stuff not in use
868           
869            #Upward slope at inlet to the north west
870            #if x[i] < 0.0: # and y[i] > 0.5:
871            #    #z[i] = -y[i]+0.5  #-x[i]/2
872            #    z[i] = x[i]/4 - y[i]**2 + 0.5
873
874            #Hole to the west
875            #x0 = -0.4; y0 = 0.35 # center
876            #if sqrt((2*(x[i]-x0))**2 + (2*(y[i]-y0))**2) < 0.2:
877            #    z[i] = sqrt(((x[i]-x0))**2 + ((y[i]-y0))**2)-0.2
878
879
880
881
882
883        return z/2
884
885class Weir_simple:
886    """Set a bathymetry for weir with a hole and a downstream gutter
887    x,y are assumed to be in the unit square
888    """
889   
890    def __init__(self, stage):
891        self.inflow_stage = stage
892
893    def __call__(self, x, y):   
894        from Numeric import zeros, Float
895   
896        N = len(x)
897        assert N == len(y)   
898
899        z = zeros(N, Float)
900        for i in range(N):
901            z[i] = -x[i]  #General slope
902
903            #Flat bit to the left
904            if x[i] < 0.3:
905                z[i] = -x[i]/10  #General slope
906           
907            #Weir
908            if x[i] > 0.3 and x[i] < 0.4:
909                z[i] = -x[i]+0.9
910
911            #Dip
912            if x[i] > 0.6 and x[i] < 0.9:
913                z[i] = -x[i]-0.5  #-y[i]/5
914
915            #Hole in weir (slightly higher than inflow condition)
916            if x[i] > 0.3 and x[i] < 0.4 and y[i] > 0.2 and y[i] < 0.4:
917                z[i] = -x[i]+self.inflow_stage + 0.05       
918           
919
920        return z/2
921       
922
923           
924class Constant_height:
925    """Set an initial condition with constant water height, e.g
926    stage s = z+h
927    """
928    def __init__(self, W, h):
929        self.W = W
930        self.h = h
931
932    def __call__(self, x, y):
933        if self.W is None:
934            from Numeric import ones, Float
935            return self.h*ones(len(x), Float)
936        else:
937            return self.W(x,y) + self.h
938
939
940##############################################
941#Initialise module
942
943
944import compile
945if compile.can_use_C_extension('shallow_water_ext.c'):
946    #Replace python version with c implementations
947   
948    from shallow_water_ext import rotate
949    compute_fluxes = compute_fluxes_c
950    gravity = gravity_c
951    manning_friction = manning_friction_c
952    balance_deep_and_shallow = balance_deep_and_shallow_c
953    protect_against_infintesimal_and_negative_heights = protect_against_infintesimal_and_negative_heights_c
954   
955    #distribute_to_vertices_and_edges = distribute_to_vertices_and_edges_c
956    #update_conserved_quantities = update_conserved_quantities_c   
957
958
959#Optimisation with psyco
960from config import use_psyco
961if use_psyco:
962    try:
963        import psyco
964    except:
965        msg = 'WARNING: psyco (speedup) could not import'+\
966              ', you may want to consider installing it'
967        print msg
968    else:
969        psyco.bind(Domain.distribute_to_vertices_and_edges)
970        psyco.bind(Domain.compute_fluxes)
971       
972if __name__ == "__main__":
973    pass
Note: See TracBrowser for help on using the repository browser.