source: inundation/pyvolution/shallow_water_kinetic.py @ 2259

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

Better actually add the new files!

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