source: anuga_core/source/anuga/abstract_2d_finite_volumes/shallow_water.py @ 3560

Last change on this file since 3560 was 3560, checked in by ole, 18 years ago

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