source: anuga_core/source/anuga/shallow_water/shallow_water_domain.py @ 3678

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

Added more limiting to cells near dry cells, use beta_*_dry

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