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

Last change on this file since 5224 was 5186, checked in by ole, 17 years ago

Changed header according to style (Peter Petkovic)
  • Property svn:keywords set to LastChangedDate LastChangedRevision LastChangedBy HeadURL Id
File size: 55.4 KB
Line 
1"""Finite-volume computations of the shallow water wave equation.
2
3Title: ANGUA shallow_water_domain - 2D triangular domains for finite-volume
4       computations of the shallow water wave equation.
5
6
7Author: Ole Nielsen, Ole.Nielsen@ga.gov.au
8        Stephen Roberts, Stephen.Roberts@anu.edu.au
9        Duncan Gray, Duncan.Gray@ga.gov.au
10
11CreationDate: 2004
12
13Description:
14    This module contains a specialisation of class Domain from
15    module domain.py consisting of methods specific to the
16    Shallow Water Wave Equation
17
18    U_t + E_x + G_y = S
19
20    where
21
22    U = [w, uh, vh]
23    E = [uh, u^2h + gh^2/2, uvh]
24    G = [vh, uvh, v^2h + gh^2/2]
25    S represents source terms forcing the system
26    (e.g. gravity, friction, wind stress, ...)
27
28    and _t, _x, _y denote the derivative with respect to t, x and y
29    respectively.
30
31
32    The quantities are
33
34    symbol    variable name    explanation
35    x         x                horizontal distance from origin [m]
36    y         y                vertical distance from origin [m]
37    z         elevation        elevation of bed on which flow is modelled [m]
38    h         height           water height above z [m]
39    w         stage            absolute water level, w = z+h [m]
40    u                          speed in the x direction [m/s]
41    v                          speed in the y direction [m/s]
42    uh        xmomentum        momentum in the x direction [m^2/s]
43    vh        ymomentum        momentum in the y direction [m^2/s]
44
45    eta                        mannings friction coefficient [to appear]
46    nu                         wind stress coefficient [to appear]
47
48    The conserved quantities are w, uh, vh
49
50Reference:
51    Catastrophic Collapse of Water Supply Reservoirs in Urban Areas,
52    Christopher Zoppou and Stephen Roberts,
53    Journal of Hydraulic Engineering, vol. 127, No. 7 July 1999
54
55    Hydrodynamic modelling of coastal inundation.
56    Nielsen, O., S. Roberts, D. Gray, A. McPherson and A. Hitchman
57    In Zerger, A. and Argent, R.M. (eds) MODSIM 2005 International Congress on
58    Modelling and Simulation. Modelling and Simulation Society of Australia and
59    New Zealand, December 2005, pp. 518-523. ISBN: 0-9758400-2-9.
60    http://www.mssanz.org.au/modsim05/papers/nielsen.pdf
61
62
63SeeAlso:
64    TRAC administration of ANUGA (User Manuals etc) at
65    https://datamining.anu.edu.au/anuga and Subversion repository at
66    $HeadURL: anuga_core/source/anuga/shallow_water/shallow_water_domain.py $
67
68Constraints: See GPL license in the user guide
69Version: 1.0 ($Revision: 5186 $)
70ModifiedBy:
71    $Author: ole $
72    $Date: 2008-04-01 23:36:55 +0000 (Tue, 01 Apr 2008) $
73
74"""
75
76# Subversion keywords:
77#
78# $LastChangedDate: 2008-04-01 23:36:55 +0000 (Tue, 01 Apr 2008) $
79# $LastChangedRevision: 5186 $
80# $LastChangedBy: ole $
81
82from Numeric import zeros, ones, Float, array, sum, size
83from Numeric import compress, arange
84
85
86from anuga.abstract_2d_finite_volumes.domain import Domain as Generic_Domain
87from anuga.abstract_2d_finite_volumes.generic_boundary_conditions\
88     import Boundary
89from anuga.abstract_2d_finite_volumes.generic_boundary_conditions\
90     import File_boundary
91from anuga.abstract_2d_finite_volumes.generic_boundary_conditions\
92     import Dirichlet_boundary
93from anuga.abstract_2d_finite_volumes.generic_boundary_conditions\
94     import Time_boundary
95from anuga.abstract_2d_finite_volumes.generic_boundary_conditions\
96     import Transmissive_boundary
97
98from anuga.utilities.numerical_tools import gradient, mean
99from anuga.config import minimum_storable_height
100from anuga.config import minimum_allowed_height, maximum_allowed_speed
101from anuga.config import g, epsilon, beta_h, beta_w, beta_w_dry,\
102     beta_uh, beta_uh_dry, beta_vh, beta_vh_dry, tight_slope_limiters
103from anuga.config import alpha_balance
104from anuga.config import optimise_dry_cells
105from anuga.config import optimised_gradient_limiter
106from anuga.config import use_edge_limiter
107
108#---------------------
109# Shallow water domain
110#---------------------
111class Domain(Generic_Domain):
112
113    conserved_quantities = ['stage', 'xmomentum', 'ymomentum']
114    other_quantities = ['elevation', 'friction']
115   
116    def __init__(self,
117                 coordinates=None,
118                 vertices=None,
119                 boundary=None,
120                 tagged_elements=None,
121                 geo_reference=None,
122                 use_inscribed_circle=False,
123                 mesh_filename=None,
124                 use_cache=False,
125                 verbose=False,
126                 full_send_dict=None,
127                 ghost_recv_dict=None,
128                 processor=0,
129                 numproc=1,
130                 number_of_full_nodes=None,
131                 number_of_full_triangles=None):
132
133
134        other_quantities = ['elevation', 'friction']
135        Generic_Domain.__init__(self,
136                                coordinates,
137                                vertices,
138                                boundary,
139                                Domain.conserved_quantities,
140                                Domain.other_quantities,
141                                tagged_elements,
142                                geo_reference,
143                                use_inscribed_circle,
144                                mesh_filename,
145                                use_cache,
146                                verbose,
147                                full_send_dict,
148                                ghost_recv_dict,
149                                processor,
150                                numproc,
151                                number_of_full_nodes=number_of_full_nodes,
152                                number_of_full_triangles=number_of_full_triangles) 
153
154
155        self.set_minimum_allowed_height(minimum_allowed_height)
156       
157        self.maximum_allowed_speed = maximum_allowed_speed
158        self.g = g
159        self.beta_w      = beta_w
160        self.beta_w_dry  = beta_w_dry
161        self.beta_uh     = beta_uh
162        self.beta_uh_dry = beta_uh_dry
163        self.beta_vh     = beta_vh
164        self.beta_vh_dry = beta_vh_dry
165        self.beta_h      = beta_h
166        self.alpha_balance = alpha_balance
167
168        self.tight_slope_limiters = tight_slope_limiters
169        self.optimise_dry_cells = optimise_dry_cells
170
171        self.forcing_terms.append(manning_friction_implicit)
172        self.forcing_terms.append(gravity)
173
174        # Stored output
175        self.store = True
176        self.format = 'sww'
177        self.set_store_vertices_uniquely(False)
178        self.minimum_storable_height = minimum_storable_height
179        self.quantities_to_be_stored = ['stage','xmomentum','ymomentum']
180
181        # Limiters
182        self.use_edge_limiter = use_edge_limiter
183
184        self.optimised_gradient_limiter = optimised_gradient_limiter
185               
186
187    def set_all_limiters(self, beta):
188        """Shorthand to assign one constant value [0,1[ to all limiters.
189        0 Corresponds to first order, where as larger values make use of
190        the second order scheme.
191        """
192
193        self.beta_w      = beta
194        self.beta_w_dry  = beta
195        self.quantities['stage'].beta = beta
196       
197        self.beta_uh     = beta
198        self.beta_uh_dry = beta
199        self.quantities['xmomentum'].beta = beta
200       
201        self.beta_vh     = beta
202        self.beta_vh_dry = beta
203        self.quantities['ymomentum'].beta = beta
204       
205        self.beta_h      = beta
206       
207
208    def set_store_vertices_uniquely(self, flag, reduction=None):
209        """Decide whether vertex values should be stored uniquely as
210        computed in the model or whether they should be reduced to one
211        value per vertex using self.reduction.
212        """
213
214        # FIXME (Ole): how about using the word continuous vertex values?
215        self.smooth = not flag
216
217        # Reduction operation for get_vertex_values
218        if reduction is None:
219            self.reduction = mean
220            #self.reduction = min  #Looks better near steep slopes
221
222
223    def set_minimum_storable_height(self, minimum_storable_height):
224        """
225        Set the minimum depth that will be recognised when writing
226        to an sww file. This is useful for removing thin water layers
227        that seems to be caused by friction creep.
228
229        The minimum allowed sww depth is in meters.
230        """
231        self.minimum_storable_height = minimum_storable_height
232
233
234    def set_minimum_allowed_height(self, minimum_allowed_height):
235        """
236        Set the minimum depth that will be recognised in the numerical
237        scheme
238
239        The minimum allowed depth is in meters.
240
241        The parameter H0 (Minimal height for flux computation)
242        is also set by this function
243        """
244
245        #FIXME (Ole): rename H0 to minimum_allowed_height_in_flux_computation
246
247        #FIXME (Ole): Maybe use histogram to identify isolated extreme speeds
248        #and deal with them adaptively similarly to how we used to use 1 order
249        #steps to recover.
250        self.minimum_allowed_height = minimum_allowed_height
251        self.H0 = minimum_allowed_height   
252       
253
254    def set_maximum_allowed_speed(self, maximum_allowed_speed):
255        """
256        Set the maximum particle speed that is allowed in water
257        shallower than minimum_allowed_height. This is useful for
258        controlling speeds in very thin layers of water and at the same time
259        allow some movement avoiding pooling of water.
260
261        """
262        self.maximum_allowed_speed = maximum_allowed_speed
263
264
265    def set_points_file_block_line_size(self,points_file_block_line_size):
266        """
267        Set the minimum depth that will be recognised when writing
268        to an sww file. This is useful for removing thin water layers
269        that seems to be caused by friction creep.
270
271        The minimum allowed sww depth is in meters.
272        """
273        self.points_file_block_line_size = points_file_block_line_size
274       
275       
276    def set_quantities_to_be_stored(self, q):
277        """Specify which quantities will be stored in the sww file.
278
279        q must be either:
280          - the name of a quantity
281          - a list of quantity names
282          - None
283
284        In the two first cases, the named quantities will be stored at
285        each yieldstep (This is in addition to the quantities elevation
286        and friction)
287       
288        If q is None, storage will be switched off altogether.
289        """
290
291
292        if q is None:
293            self.quantities_to_be_stored = []
294            self.store = False
295            return
296
297        if isinstance(q, basestring):
298            q = [q] # Turn argument into a list
299
300        # Check correcness
301        for quantity_name in q:
302            msg = 'Quantity %s is not a valid conserved quantity'\
303                  %quantity_name
304           
305            assert quantity_name in self.conserved_quantities, msg
306
307        self.quantities_to_be_stored = q
308
309
310
311    def get_wet_elements(self, indices=None):
312        """Return indices for elements where h > minimum_allowed_height
313
314        Optional argument:
315            indices is the set of element ids that the operation applies to.
316
317        Usage:
318            indices = get_wet_elements()
319
320        Note, centroid values are used for this operation           
321        """
322
323        # Water depth below which it is considered to be 0 in the model
324        # FIXME (Ole): Allow this to be specified as a keyword argument as well
325        from anuga.config import minimum_allowed_height
326
327       
328        elevation = self.get_quantity('elevation').\
329                    get_values(location='centroids', indices=indices)
330        stage = self.get_quantity('stage').\
331                get_values(location='centroids', indices=indices)       
332        depth = stage - elevation
333
334        # Select indices for which depth > 0
335        wet_indices = compress(depth > minimum_allowed_height,
336                               arange(len(depth)))
337        return wet_indices
338
339
340    def get_maximum_inundation_elevation(self, indices=None):
341        """Return highest elevation where h > 0
342
343        Optional argument:
344            indices is the set of element ids that the operation applies to.
345
346        Usage:
347            q = get_maximum_inundation_elevation()
348
349        Note, centroid values are used for this operation           
350        """
351
352        wet_elements = self.get_wet_elements(indices)
353        return self.get_quantity('elevation').\
354               get_maximum_value(indices=wet_elements)
355
356
357    def get_maximum_inundation_location(self, indices=None):
358        """Return location of highest elevation where h > 0
359
360        Optional argument:
361            indices is the set of element ids that the operation applies to.
362
363        Usage:
364            q = get_maximum_inundation_location()
365
366        Note, centroid values are used for this operation           
367        """
368
369        wet_elements = self.get_wet_elements(indices)
370        return self.get_quantity('elevation').\
371               get_maximum_location(indices=wet_elements)   
372
373    def check_integrity(self):
374        Generic_Domain.check_integrity(self)
375
376        #Check that we are solving the shallow water wave equation
377
378        msg = 'First conserved quantity must be "stage"'
379        assert self.conserved_quantities[0] == 'stage', msg
380        msg = 'Second conserved quantity must be "xmomentum"'
381        assert self.conserved_quantities[1] == 'xmomentum', msg
382        msg = 'Third conserved quantity must be "ymomentum"'
383        assert self.conserved_quantities[2] == 'ymomentum', msg
384
385    def extrapolate_second_order_sw(self):
386        #Call correct module function
387        #(either from this module or C-extension)
388        extrapolate_second_order_sw(self)
389
390    def compute_fluxes(self):
391        #Call correct module function
392        #(either from this module or C-extension)
393        compute_fluxes(self)
394
395    def distribute_to_vertices_and_edges(self):
396        # Call correct module function
397        # (either from this module or C-extension)
398        if self.use_edge_limiter:
399            protect_against_infinitesimal_and_negative_heights(self)
400            for name in self.conserved_quantities:
401                Q = self.quantities[name]
402                if self._order_ == 1:
403                    Q.extrapolate_first_order()
404                elif self._order_ == 2:
405                    Q.extrapolate_second_order_and_limit()
406                else:
407                    raise 'Unknown order'
408            balance_deep_and_shallow(self)
409
410            #Compute edge values by interpolation
411            for name in self.conserved_quantities:
412                Q = self.quantities[name]
413                Q.interpolate_from_vertices_to_edges()
414        else:
415            distribute_to_vertices_and_edges(self)
416
417
418
419
420    def evolve(self,
421               yieldstep = None,
422               finaltime = None,
423               duration = None,
424               skip_initial_step = False):
425        """Specialisation of basic evolve method from parent class
426        """
427
428        # Call check integrity here rather than from user scripts
429        # self.check_integrity()
430
431        msg = 'Parameter beta_h must be in the interval [0, 2['
432        assert 0 <= self.beta_h <= 2.0, msg
433        msg = 'Parameter beta_w must be in the interval [0, 2['
434        assert 0 <= self.beta_w <= 2.0, msg
435
436
437        # Initial update of vertex and edge values before any STORAGE
438        # and or visualisation
439        # This is done again in the initialisation of the Generic_Domain
440        # evolve loop but we do it here to ensure the values are ok for storage
441        self.distribute_to_vertices_and_edges()
442
443        if self.store is True and self.time == 0.0:
444            self.initialise_storage()
445            # print 'Storing results in ' + self.writer.filename
446        else:
447            pass
448            # print 'Results will not be stored.'
449            # print 'To store results set domain.store = True'
450            # FIXME: Diagnostic output should be controlled by
451            # a 'verbose' flag living in domain (or in a parent class)
452
453        # Call basic machinery from parent class
454        for t in Generic_Domain.evolve(self,
455                                       yieldstep=yieldstep,
456                                       finaltime=finaltime,
457                                       duration=duration,
458                                       skip_initial_step=skip_initial_step):
459
460            # Store model data, e.g. for subsequent visualisation
461            if self.store is True:
462                self.store_timestep(self.quantities_to_be_stored)
463
464            # FIXME: Could maybe be taken from specified list
465            # of 'store every step' quantities
466
467            # Pass control on to outer loop for more specific actions
468            yield(t)
469
470
471    def initialise_storage(self):
472        """Create and initialise self.writer object for storing data.
473        Also, save x,y and bed elevation
474        """
475
476        from anuga.shallow_water.data_manager import get_dataobject
477
478        # Initialise writer
479        self.writer = get_dataobject(self, mode = 'w')
480
481        # Store vertices and connectivity
482        self.writer.store_connectivity()
483
484
485    def store_timestep(self, name):
486        """Store named quantity and time.
487
488        Precondition:
489           self.write has been initialised
490        """
491        self.writer.store_timestep(name)
492
493       
494    def timestepping_statistics(self,
495                                track_speeds=False,
496                                triangle_id=None):       
497        """Return string with time stepping statistics for printing or logging
498
499        Optional boolean keyword track_speeds decides whether to report
500        location of smallest timestep as well as a histogram and percentile
501        report.
502        """
503
504        from Numeric import sqrt
505        from anuga.config import epsilon, g               
506
507
508        # Call basic machinery from parent class
509        msg = Generic_Domain.timestepping_statistics(self,
510                                                     track_speeds,
511                                                     triangle_id)
512
513        if track_speeds is True:
514
515            # qwidth determines the text field used for quantities
516            qwidth = self.qwidth
517       
518            # Selected triangle
519            k = self.k
520
521            # Report some derived quantities at vertices, edges and centroid
522            # specific to the shallow water wave equation
523
524            z = self.quantities['elevation']
525            w = self.quantities['stage']           
526
527            Vw = w.get_values(location='vertices', indices=[k])[0]
528            Ew = w.get_values(location='edges', indices=[k])[0]
529            Cw = w.get_values(location='centroids', indices=[k])
530
531            Vz = z.get_values(location='vertices', indices=[k])[0]
532            Ez = z.get_values(location='edges', indices=[k])[0]
533            Cz = z.get_values(location='centroids', indices=[k])                             
534               
535
536            name = 'depth'
537            Vh = Vw-Vz
538            Eh = Ew-Ez
539            Ch = Cw-Cz
540           
541            s  = '    %s: vertex_values =  %.4f,\t %.4f,\t %.4f\n'\
542                 %(name.ljust(qwidth), Vh[0], Vh[1], Vh[2])
543           
544            s += '    %s: edge_values =    %.4f,\t %.4f,\t %.4f\n'\
545                 %(name.ljust(qwidth), Eh[0], Eh[1], Eh[2])
546           
547            s += '    %s: centroid_value = %.4f\n'\
548                 %(name.ljust(qwidth), Ch[0])                               
549           
550            msg += s
551
552            uh = self.quantities['xmomentum']
553            vh = self.quantities['ymomentum']
554
555            Vuh = uh.get_values(location='vertices', indices=[k])[0]
556            Euh = uh.get_values(location='edges', indices=[k])[0]
557            Cuh = uh.get_values(location='centroids', indices=[k])
558           
559            Vvh = vh.get_values(location='vertices', indices=[k])[0]
560            Evh = vh.get_values(location='edges', indices=[k])[0]
561            Cvh = vh.get_values(location='centroids', indices=[k])
562
563            # Speeds in each direction
564            Vu = Vuh/(Vh + epsilon)
565            Eu = Euh/(Eh + epsilon)
566            Cu = Cuh/(Ch + epsilon)           
567            name = 'U'
568            s  = '    %s: vertex_values =  %.4f,\t %.4f,\t %.4f\n'\
569                 %(name.ljust(qwidth), Vu[0], Vu[1], Vu[2])
570           
571            s += '    %s: edge_values =    %.4f,\t %.4f,\t %.4f\n'\
572                 %(name.ljust(qwidth), Eu[0], Eu[1], Eu[2])
573           
574            s += '    %s: centroid_value = %.4f\n'\
575                 %(name.ljust(qwidth), Cu[0])                               
576           
577            msg += s
578
579            Vv = Vvh/(Vh + epsilon)
580            Ev = Evh/(Eh + epsilon)
581            Cv = Cvh/(Ch + epsilon)           
582            name = 'V'
583            s  = '    %s: vertex_values =  %.4f,\t %.4f,\t %.4f\n'\
584                 %(name.ljust(qwidth), Vv[0], Vv[1], Vv[2])
585           
586            s += '    %s: edge_values =    %.4f,\t %.4f,\t %.4f\n'\
587                 %(name.ljust(qwidth), Ev[0], Ev[1], Ev[2])
588           
589            s += '    %s: centroid_value = %.4f\n'\
590                 %(name.ljust(qwidth), Cv[0])                               
591           
592            msg += s
593
594
595            # Froude number in each direction
596            name = 'Froude (x)'
597            Vfx = Vu/(sqrt(g*Vh) + epsilon)
598            Efx = Eu/(sqrt(g*Eh) + epsilon)
599            Cfx = Cu/(sqrt(g*Ch) + epsilon)
600           
601            s  = '    %s: vertex_values =  %.4f,\t %.4f,\t %.4f\n'\
602                 %(name.ljust(qwidth), Vfx[0], Vfx[1], Vfx[2])
603           
604            s += '    %s: edge_values =    %.4f,\t %.4f,\t %.4f\n'\
605                 %(name.ljust(qwidth), Efx[0], Efx[1], Efx[2])
606           
607            s += '    %s: centroid_value = %.4f\n'\
608                 %(name.ljust(qwidth), Cfx[0])                               
609           
610            msg += s
611
612
613            name = 'Froude (y)'
614            Vfy = Vv/(sqrt(g*Vh) + epsilon)
615            Efy = Ev/(sqrt(g*Eh) + epsilon)
616            Cfy = Cv/(sqrt(g*Ch) + epsilon)
617           
618            s  = '    %s: vertex_values =  %.4f,\t %.4f,\t %.4f\n'\
619                 %(name.ljust(qwidth), Vfy[0], Vfy[1], Vfy[2])
620           
621            s += '    %s: edge_values =    %.4f,\t %.4f,\t %.4f\n'\
622                 %(name.ljust(qwidth), Efy[0], Efy[1], Efy[2])
623           
624            s += '    %s: centroid_value = %.4f\n'\
625                 %(name.ljust(qwidth), Cfy[0])                               
626           
627            msg += s           
628
629               
630
631        return msg
632       
633       
634
635#=============== End of class Shallow Water Domain ===============================
636
637
638#-----------------
639# Flux computation
640#-----------------
641
642def compute_fluxes(domain):
643    """Compute all fluxes and the timestep suitable for all volumes
644    in domain.
645
646    Compute total flux for each conserved quantity using "flux_function"
647
648    Fluxes across each edge are scaled by edgelengths and summed up
649    Resulting flux is then scaled by area and stored in
650    explicit_update for each of the three conserved quantities
651    stage, xmomentum and ymomentum
652
653    The maximal allowable speed computed by the flux_function for each volume
654    is converted to a timestep that must not be exceeded. The minimum of
655    those is computed as the next overall timestep.
656
657    Post conditions:
658      domain.explicit_update is reset to computed flux values
659      domain.timestep is set to the largest step satisfying all volumes.
660   
661
662    This wrapper calls the underlying C version of compute fluxes
663    """
664
665    import sys
666
667    N = len(domain) # number_of_triangles
668
669    # Shortcuts
670    Stage = domain.quantities['stage']
671    Xmom = domain.quantities['xmomentum']
672    Ymom = domain.quantities['ymomentum']
673    Bed = domain.quantities['elevation']
674
675    timestep = float(sys.maxint)
676    from shallow_water_ext import\
677         compute_fluxes_ext_central as compute_fluxes_ext
678
679
680    flux_timestep = compute_fluxes_ext(timestep,
681                                       domain.epsilon,
682                                       domain.H0,
683                                       domain.g,
684                                       domain.neighbours,
685                                       domain.neighbour_edges,
686                                       domain.normals,
687                                       domain.edgelengths,
688                                       domain.radii,
689                                       domain.areas,
690                                       domain.tri_full_flag,
691                                       Stage.edge_values,
692                                       Xmom.edge_values,
693                                       Ymom.edge_values,
694                                       Bed.edge_values,
695                                       Stage.boundary_values,
696                                       Xmom.boundary_values,
697                                       Ymom.boundary_values,
698                                       Stage.explicit_update,
699                                       Xmom.explicit_update,
700                                       Ymom.explicit_update,
701                                       domain.already_computed_flux,
702                                       domain.max_speed,
703                                       int(domain.optimise_dry_cells))
704
705    domain.flux_timestep = flux_timestep
706
707
708
709#---------------------------------------
710# Module functions for gradient limiting
711#---------------------------------------
712
713
714# MH090605 The following method belongs to the shallow_water domain class
715# see comments in the corresponding method in shallow_water_ext.c
716def extrapolate_second_order_sw(domain):
717    """Wrapper calling C version of extrapolate_second_order_sw
718    """
719    import sys
720
721    N = len(domain) # number_of_triangles
722
723    # Shortcuts
724    Stage = domain.quantities['stage']
725    Xmom = domain.quantities['xmomentum']
726    Ymom = domain.quantities['ymomentum']
727    Elevation = domain.quantities['elevation']
728
729    from shallow_water_ext import extrapolate_second_order_sw as extrapol2
730    extrapol2(domain,
731              domain.surrogate_neighbours,
732              domain.number_of_boundaries,
733              domain.centroid_coordinates,
734              Stage.centroid_values,
735              Xmom.centroid_values,
736              Ymom.centroid_values,
737              Elevation.centroid_values,
738              domain.vertex_coordinates,
739              Stage.vertex_values,
740              Xmom.vertex_values,
741              Ymom.vertex_values,
742              Elevation.vertex_values,
743              int(domain.optimise_dry_cells))
744
745
746def distribute_to_vertices_and_edges(domain):
747    """Distribution from centroids to vertices specific to the
748    shallow water wave
749    equation.
750
751    It will ensure that h (w-z) is always non-negative even in the
752    presence of steep bed-slopes by taking a weighted average between shallow
753    and deep cases.
754
755    In addition, all conserved quantities get distributed as per either a
756    constant (order==1) or a piecewise linear function (order==2).
757
758    FIXME: more explanation about removal of artificial variability etc
759
760    Precondition:
761      All quantities defined at centroids and bed elevation defined at
762      vertices.
763
764    Postcondition
765      Conserved quantities defined at vertices
766
767    """
768
769   
770
771    # Remove very thin layers of water
772    protect_against_infinitesimal_and_negative_heights(domain)
773
774    # Extrapolate all conserved quantities
775    if domain.optimised_gradient_limiter:
776        # MH090605 if second order,
777        # perform the extrapolation and limiting on
778        # all of the conserved quantities
779
780        if (domain._order_ == 1):
781            for name in domain.conserved_quantities:
782                Q = domain.quantities[name]
783                Q.extrapolate_first_order()
784        elif domain._order_ == 2:
785            domain.extrapolate_second_order_sw()
786        else:
787            raise 'Unknown order'
788    else:
789        # Old code:
790        for name in domain.conserved_quantities:
791            Q = domain.quantities[name]
792
793            if domain._order_ == 1:
794                Q.extrapolate_first_order()
795            elif domain._order_ == 2:
796                Q.extrapolate_second_order()
797                Q.limit()
798            else:
799                raise 'Unknown order'
800
801
802    #Take bed elevation into account when water heights are small
803    balance_deep_and_shallow(domain)
804
805    #Compute edge values by interpolation
806    for name in domain.conserved_quantities:
807        Q = domain.quantities[name]
808        Q.interpolate_from_vertices_to_edges()
809
810
811def protect_against_infinitesimal_and_negative_heights(domain):
812    """Protect against infinitesimal heights and associated high velocities
813    """
814
815    # Shortcuts
816    wc = domain.quantities['stage'].centroid_values
817    zc = domain.quantities['elevation'].centroid_values
818    xmomc = domain.quantities['xmomentum'].centroid_values
819    ymomc = domain.quantities['ymomentum'].centroid_values
820
821    from shallow_water_ext import protect
822
823    protect(domain.minimum_allowed_height, domain.maximum_allowed_speed,
824            domain.epsilon, wc, zc, xmomc, ymomc)
825
826
827def h_limiter(domain):
828    """Limit slopes for each volume to eliminate artificial variance
829    introduced by e.g. second order extrapolator
830
831    limit on h = w-z
832
833    This limiter depends on two quantities (w,z) so it resides within
834    this module rather than within quantity.py
835
836    Wrapper for c-extension
837    """
838
839    N = len(domain) # number_of_triangles
840    beta_h = domain.beta_h
841
842    # Shortcuts
843    wc = domain.quantities['stage'].centroid_values
844    zc = domain.quantities['elevation'].centroid_values
845    hc = wc - zc
846
847    wv = domain.quantities['stage'].vertex_values
848    zv = domain.quantities['elevation'].vertex_values
849    hv = wv - zv
850
851    #Call C-extension
852    from shallow_water_ext import h_limiter_sw
853    hvbar = h_limiter_sw(domain, hc, hv)
854
855    return hvbar
856
857
858def balance_deep_and_shallow(domain):
859    """Compute linear combination between stage as computed by
860    gradient-limiters limiting using w, and stage computed by
861    gradient-limiters limiting using h (h-limiter).
862    The former takes precedence when heights are large compared to the
863    bed slope while the latter takes precedence when heights are
864    relatively small.  Anything in between is computed as a balanced
865    linear combination in order to avoid numerical disturbances which
866    would otherwise appear as a result of hard switching between
867    modes.
868
869    Wrapper for C implementation
870    """
871
872    # FIXME (Ole): I reckon this can be simplified significantly:
873    #
874    # Always use beta_h == 0, and phase it out.
875    # Compute hc and hv in the c-code
876    # Omit updating xmomv
877    #
878    from shallow_water_ext import balance_deep_and_shallow as balance_deep_and_shallow_c
879
880
881    #print 'calling balance depth and shallow'
882    # Shortcuts
883    wc = domain.quantities['stage'].centroid_values
884    zc = domain.quantities['elevation'].centroid_values
885
886    wv = domain.quantities['stage'].vertex_values
887    zv = domain.quantities['elevation'].vertex_values
888
889    # Momentums at centroids
890    xmomc = domain.quantities['xmomentum'].centroid_values
891    ymomc = domain.quantities['ymomentum'].centroid_values
892
893    # Momentums at vertices
894    xmomv = domain.quantities['xmomentum'].vertex_values
895    ymomv = domain.quantities['ymomentum'].vertex_values
896
897    # Limit h
898    if domain.beta_h > 0:
899        hvbar = h_limiter(domain)
900       
901        balance_deep_and_shallow_c(domain, domain.beta_h,
902                                   wc, zc, wv, zv, hvbar,
903                                   xmomc, ymomc, xmomv, ymomv)       
904    else:
905        # print 'Using first order h-limiter'
906        # FIXME: Pass wc in for now - it will be ignored.
907       
908        # This is how one would make a first order h_limited value
909        # as in the old balancer (pre 17 Feb 2005):
910        #  If we wish to hard wire this, one should modify the C-code
911        # from Numeric import zeros, Float
912        # hvbar = zeros( (len(wc), 3), Float)
913        # for i in range(3):
914        #     hvbar[:,i] = wc[:] - zc[:]
915
916        balance_deep_and_shallow_c(domain, domain.beta_h,
917                                   wc, zc, wv, zv, wc, 
918                                   xmomc, ymomc, xmomv, ymomv)
919
920
921
922
923#------------------------------------------------------------------
924# Boundary conditions - specific to the shallow water wave equation
925#------------------------------------------------------------------
926class Reflective_boundary(Boundary):
927    """Reflective boundary returns same conserved quantities as
928    those present in its neighbour volume but reflected.
929
930    This class is specific to the shallow water equation as it
931    works with the momentum quantities assumed to be the second
932    and third conserved quantities.
933    """
934
935    def __init__(self, domain = None):
936        Boundary.__init__(self)
937
938        if domain is None:
939            msg = 'Domain must be specified for reflective boundary'
940            raise msg
941
942        # Handy shorthands
943        self.stage   = domain.quantities['stage'].edge_values
944        self.xmom    = domain.quantities['xmomentum'].edge_values
945        self.ymom    = domain.quantities['ymomentum'].edge_values
946        self.normals = domain.normals
947
948        self.conserved_quantities = zeros(3, Float)
949
950    def __repr__(self):
951        return 'Reflective_boundary'
952
953
954    def evaluate(self, vol_id, edge_id):
955        """Reflective boundaries reverses the outward momentum
956        of the volume they serve.
957        """
958
959        q = self.conserved_quantities
960        q[0] = self.stage[vol_id, edge_id]
961        q[1] = self.xmom[vol_id, edge_id]
962        q[2] = self.ymom[vol_id, edge_id]
963
964        normal = self.normals[vol_id, 2*edge_id:2*edge_id+2]
965
966
967        r = rotate(q, normal, direction = 1)
968        r[1] = -r[1]
969        q = rotate(r, normal, direction = -1)
970
971        return q
972
973
974
975class Transmissive_Momentum_Set_Stage_boundary(Boundary):
976    """Returns same momentum conserved quantities as
977    those present in its neighbour volume.
978    Sets stage by specifying a function f of time which may either be a
979    vector function or a scalar function
980
981    Example:
982
983    def waveform(t):
984        return sea_level + normalized_amplitude/cosh(t-25)**2
985
986    Bts = Transmissive_Momentum_Set_Stage_boundary(domain, waveform)
987   
988
989    Underlying domain must be specified when boundary is instantiated
990    """
991
992    def __init__(self, domain = None, function=None):
993        Boundary.__init__(self)
994
995        if domain is None:
996            msg = 'Domain must be specified for this type boundary'
997            raise msg
998
999        if function is None:
1000            msg = 'Function must be specified for this type boundary'
1001            raise msg
1002
1003        self.domain   = domain
1004        self.function = function
1005
1006    def __repr__(self):
1007        return 'Transmissive_Momentum_Set_Stage_boundary(%s)' %self.domain
1008
1009    def evaluate(self, vol_id, edge_id):
1010        """Transmissive Momentum Set Stage boundaries return the edge momentum
1011        values of the volume they serve.
1012        """
1013
1014        q = self.domain.get_conserved_quantities(vol_id, edge = edge_id)
1015
1016
1017        t = self.domain.time
1018
1019        if hasattr(self.function, 'time'):
1020            # Roll boundary over if time exceeds           
1021            while t > self.function.time[-1]:
1022                msg = 'WARNING: domain time %.2f has exceeded' %t
1023                msg += 'time provided in '
1024                msg += 'transmissive_momentum_set_stage boundary object.\n'
1025                msg += 'I will continue, reusing the object from t==0'
1026                print msg
1027                t -= self.function.time[-1]
1028
1029
1030        value = self.function(t)
1031        try:
1032            x = float(value)
1033        except:
1034            x = float(value[0])
1035           
1036        q[0] = x
1037        return q
1038
1039
1040        # FIXME: Consider this (taken from File_boundary) to allow
1041        # spatial variation
1042        # if vol_id is not None and edge_id is not None:
1043        #     i = self.boundary_indices[ vol_id, edge_id ]
1044        #     return self.F(t, point_id = i)
1045        # else:
1046        #     return self.F(t)
1047
1048
1049
1050class Dirichlet_Discharge_boundary(Boundary):
1051    """
1052    Sets stage (stage0)
1053    Sets momentum (wh0) in the inward normal direction.
1054
1055    Underlying domain must be specified when boundary is instantiated
1056    """
1057
1058    def __init__(self, domain = None, stage0=None, wh0=None):
1059        Boundary.__init__(self)
1060
1061        if domain is None:
1062            msg = 'Domain must be specified for this type boundary'
1063            raise msg
1064
1065        if stage0 is None:
1066            raise 'set stage'
1067
1068        if wh0 is None:
1069            wh0 = 0.0
1070
1071        self.domain   = domain
1072        self.stage0  = stage0
1073        self.wh0 = wh0
1074
1075    def __repr__(self):
1076        return 'Dirichlet_Discharge_boundary(%s)' %self.domain
1077
1078    def evaluate(self, vol_id, edge_id):
1079        """Set discharge in the (inward) normal direction
1080        """
1081
1082        normal = self.domain.get_normal(vol_id,edge_id)
1083        q = [self.stage0, -self.wh0*normal[0], -self.wh0*normal[1]]
1084        return q
1085
1086
1087        # FIXME: Consider this (taken from File_boundary) to allow
1088        # spatial variation
1089        # if vol_id is not None and edge_id is not None:
1090        #     i = self.boundary_indices[ vol_id, edge_id ]
1091        #     return self.F(t, point_id = i)
1092        # else:
1093        #     return self.F(t)
1094
1095
1096class Field_boundary(Boundary):
1097    """Set boundary from given field represented in an sww file containing values
1098    for stage, xmomentum and ymomentum.
1099    Optionally, the user can specify mean_stage to offset the stage provided in the
1100    sww file.
1101
1102    This function is a thin wrapper around the generic File_boundary. The
1103    difference between the file_boundary and field_boundary is only that the
1104    field_boundary will allow you to change the level of the stage height when
1105    you read in the boundary condition. This is very useful when running
1106    different tide heights in the same area as you need only to convert one
1107    boundary condition to a SWW file, ideally for tide height of 0 m
1108    (saving disk space). Then you can use field_boundary to read this SWW file
1109    and change the stage height (tide) on the fly depending on the scenario.
1110   
1111    """
1112
1113
1114    def __init__(self, filename, domain,
1115                 mean_stage=0.0,
1116                 time_thinning=1, 
1117                 use_cache=False,
1118                 verbose=False):
1119        """Constructor
1120
1121        filename: Name of sww file
1122        domain: pointer to shallow water domain for which the boundary applies
1123        mean_stage: The mean water level which will be added to stage derived
1124                    from the sww file
1125        time_thinning: Will set how many time steps from the sww file read in
1126                       will be interpolated to the boundary. For example if
1127                       the sww file has 1 second time steps and is 24 hours
1128                       in length it has 86400 time steps. If you set
1129                       time_thinning to 1 it will read all these steps.
1130                       If you set it to 100 it will read every 100th step eg
1131                       only 864 step. This parameter is very useful to increase
1132                       the speed of a model run that you are setting up
1133                       and testing.
1134        use_cache:
1135        verbose:
1136       
1137        """
1138
1139        # Create generic file_boundary object
1140        self.file_boundary = File_boundary(filename, domain,
1141                                           time_thinning=time_thinning,
1142                                           use_cache=use_cache,
1143                                           verbose=verbose)
1144       
1145        # Record information from File_boundary
1146        self.F = self.file_boundary.F
1147        self.domain = self.file_boundary.domain
1148       
1149        # Record mean stage
1150        self.mean_stage = mean_stage
1151
1152
1153    def __repr__(self):
1154        return 'Field boundary'
1155
1156
1157    def evaluate(self, vol_id=None, edge_id=None):
1158        """Return linearly interpolated values based on domain.time
1159
1160        vol_id and edge_id are ignored
1161        """
1162
1163        # Evaluate file boundary
1164        q = self.file_boundary.evaluate(vol_id, edge_id)
1165
1166        # Adjust stage
1167        for j, name in enumerate(self.domain.conserved_quantities):
1168            if name == 'stage':
1169                q[j] += self.mean_stage
1170        return q
1171
1172   
1173
1174#-----------------------
1175# Standard forcing terms
1176#-----------------------
1177
1178def gravity(domain):
1179    """Apply gravitational pull in the presence of bed slope
1180    Wrapper calls underlying C implementation
1181    """
1182
1183    xmom = domain.quantities['xmomentum'].explicit_update
1184    ymom = domain.quantities['ymomentum'].explicit_update
1185
1186    stage = domain.quantities['stage']
1187    elevation = domain.quantities['elevation']
1188
1189    h = stage.centroid_values - elevation.centroid_values
1190    z = elevation.vertex_values
1191
1192    x = domain.get_vertex_coordinates()
1193    g = domain.g
1194   
1195
1196    from shallow_water_ext import gravity as gravity_c
1197    gravity_c(g, h, z, x, xmom, ymom) #, 1.0e-6)
1198
1199
1200
1201def manning_friction_implicit(domain):
1202    """Apply (Manning) friction to water momentum   
1203    Wrapper for c version
1204    """
1205
1206
1207    #print 'Implicit friction'
1208
1209    xmom = domain.quantities['xmomentum']
1210    ymom = domain.quantities['ymomentum']
1211
1212    w = domain.quantities['stage'].centroid_values
1213    z = domain.quantities['elevation'].centroid_values
1214
1215    uh = xmom.centroid_values
1216    vh = ymom.centroid_values
1217    eta = domain.quantities['friction'].centroid_values
1218
1219    xmom_update = xmom.semi_implicit_update
1220    ymom_update = ymom.semi_implicit_update
1221
1222    N = len(domain)
1223    eps = domain.minimum_allowed_height
1224    g = domain.g
1225
1226    from shallow_water_ext import manning_friction as manning_friction_c
1227    manning_friction_c(g, eps, w, z, uh, vh, eta, xmom_update, ymom_update)
1228
1229
1230def manning_friction_explicit(domain):
1231    """Apply (Manning) friction to water momentum   
1232    Wrapper for c version
1233    """
1234
1235    # print 'Explicit friction'
1236
1237    xmom = domain.quantities['xmomentum']
1238    ymom = domain.quantities['ymomentum']
1239
1240    w = domain.quantities['stage'].centroid_values
1241    z = domain.quantities['elevation'].centroid_values
1242
1243    uh = xmom.centroid_values
1244    vh = ymom.centroid_values
1245    eta = domain.quantities['friction'].centroid_values
1246
1247    xmom_update = xmom.explicit_update
1248    ymom_update = ymom.explicit_update
1249
1250    N = len(domain)
1251    eps = domain.minimum_allowed_height
1252    g = domain.g
1253
1254    from shallow_water_ext import manning_friction as manning_friction_c
1255    manning_friction_c(g, eps, w, z, uh, vh, eta, xmom_update, ymom_update)
1256
1257
1258# FIXME (Ole): This was implemented for use with one of the analytical solutions (Sampson?)
1259# Is it still needed (30 Oct 2007)?
1260def linear_friction(domain):
1261    """Apply linear friction to water momentum
1262
1263    Assumes quantity: 'linear_friction' to be present
1264    """
1265
1266    from math import sqrt
1267
1268    w = domain.quantities['stage'].centroid_values
1269    z = domain.quantities['elevation'].centroid_values
1270    h = w-z
1271
1272    uh = domain.quantities['xmomentum'].centroid_values
1273    vh = domain.quantities['ymomentum'].centroid_values
1274    tau = domain.quantities['linear_friction'].centroid_values
1275
1276    xmom_update = domain.quantities['xmomentum'].semi_implicit_update
1277    ymom_update = domain.quantities['ymomentum'].semi_implicit_update
1278
1279    N = len(domain) # number_of_triangles
1280    eps = domain.minimum_allowed_height
1281    g = domain.g #Not necessary? Why was this added?
1282
1283    for k in range(N):
1284        if tau[k] >= eps:
1285            if h[k] >= eps:
1286                S = -tau[k]/h[k]
1287
1288                #Update momentum
1289                xmom_update[k] += S*uh[k]
1290                ymom_update[k] += S*vh[k]
1291
1292
1293
1294#---------------------------------
1295# Experimental auxiliary functions
1296#---------------------------------
1297def check_forcefield(f):
1298    """Check that f is either
1299    1: a callable object f(t,x,y), where x and y are vectors
1300       and that it returns an array or a list of same length
1301       as x and y
1302    2: a scalar
1303    """
1304
1305    if callable(f):
1306        N = 3
1307        x = ones(3, Float)
1308        y = ones(3, Float)
1309        try:
1310            q = f(1.0, x=x, y=y)
1311        except Exception, e:
1312            msg = 'Function %s could not be executed:\n%s' %(f, e)
1313            # FIXME: Reconsider this semantics
1314            raise msg
1315
1316        try:
1317            q = array(q).astype(Float)
1318        except:
1319            msg = 'Return value from vector function %s could ' %f
1320            msg += 'not be converted into a Numeric array of floats.\n'
1321            msg += 'Specified function should return either list or array.'
1322            raise msg
1323
1324        # Is this really what we want?
1325        msg = 'Return vector from function %s ' %f
1326        msg += 'must have same lenght as input vectors'
1327        assert len(q) == N, msg
1328
1329    else:
1330        try:
1331            f = float(f)
1332        except:
1333            msg = 'Force field %s must be either a scalar' %f
1334            msg += ' or a vector function'
1335            raise Exception(msg)
1336    return f
1337
1338
1339class Wind_stress:
1340    """Apply wind stress to water momentum in terms of
1341    wind speed [m/s] and wind direction [degrees]
1342    """
1343
1344    def __init__(self, *args, **kwargs):
1345        """Initialise windfield from wind speed s [m/s]
1346        and wind direction phi [degrees]
1347
1348        Inputs v and phi can be either scalars or Python functions, e.g.
1349
1350        W = Wind_stress(10, 178)
1351
1352        #FIXME - 'normal' degrees are assumed for now, i.e. the
1353        vector (1,0) has zero degrees.
1354        We may need to convert from 'compass' degrees later on and also
1355        map from True north to grid north.
1356
1357        Arguments can also be Python functions of t,x,y as in
1358
1359        def speed(t,x,y):
1360            ...
1361            return s
1362
1363        def angle(t,x,y):
1364            ...
1365            return phi
1366
1367        where x and y are vectors.
1368
1369        and then pass the functions in
1370
1371        W = Wind_stress(speed, angle)
1372
1373        The instantiated object W can be appended to the list of
1374        forcing_terms as in
1375
1376        Alternatively, one vector valued function for (speed, angle)
1377        can be applied, providing both quantities simultaneously.
1378        As in
1379        W = Wind_stress(F), where returns (speed, angle) for each t.
1380
1381        domain.forcing_terms.append(W)
1382        """
1383
1384        from anuga.config import rho_a, rho_w, eta_w
1385        from Numeric import array, Float
1386
1387        if len(args) == 2:
1388            s = args[0]
1389            phi = args[1]
1390        elif len(args) == 1:
1391            # Assume vector function returning (s, phi)(t,x,y)
1392            vector_function = args[0]
1393            s = lambda t,x,y: vector_function(t,x=x,y=y)[0]
1394            phi = lambda t,x,y: vector_function(t,x=x,y=y)[1]
1395        else:
1396           # Assume info is in 2 keyword arguments
1397
1398           if len(kwargs) == 2:
1399               s = kwargs['s']
1400               phi = kwargs['phi']
1401           else:
1402               raise 'Assumes two keyword arguments: s=..., phi=....'
1403
1404        self.speed = check_forcefield(s)
1405        self.phi = check_forcefield(phi)
1406
1407        self.const = eta_w*rho_a/rho_w
1408
1409
1410    def __call__(self, domain):
1411        """Evaluate windfield based on values found in domain
1412        """
1413
1414        from math import pi, cos, sin, sqrt
1415        from Numeric import Float, ones, ArrayType
1416
1417        xmom_update = domain.quantities['xmomentum'].explicit_update
1418        ymom_update = domain.quantities['ymomentum'].explicit_update
1419
1420        N = len(domain) # number_of_triangles
1421        t = domain.time
1422
1423        if callable(self.speed):
1424            xc = domain.get_centroid_coordinates()
1425            s_vec = self.speed(t, xc[:,0], xc[:,1])
1426        else:
1427            # Assume s is a scalar
1428
1429            try:
1430                s_vec = self.speed * ones(N, Float)
1431            except:
1432                msg = 'Speed must be either callable or a scalar: %s' %self.s
1433                raise msg
1434
1435
1436        if callable(self.phi):
1437            xc = domain.get_centroid_coordinates()
1438            phi_vec = self.phi(t, xc[:,0], xc[:,1])
1439        else:
1440            # Assume phi is a scalar
1441
1442            try:
1443                phi_vec = self.phi * ones(N, Float)
1444            except:
1445                msg = 'Angle must be either callable or a scalar: %s' %self.phi
1446                raise msg
1447
1448        assign_windfield_values(xmom_update, ymom_update,
1449                                s_vec, phi_vec, self.const)
1450
1451
1452def assign_windfield_values(xmom_update, ymom_update,
1453                            s_vec, phi_vec, const):
1454    """Python version of assigning wind field to update vectors.
1455    A c version also exists (for speed)
1456    """
1457    from math import pi, cos, sin, sqrt
1458
1459    N = len(s_vec)
1460    for k in range(N):
1461        s = s_vec[k]
1462        phi = phi_vec[k]
1463
1464        # Convert to radians
1465        phi = phi*pi/180
1466
1467        # Compute velocity vector (u, v)
1468        u = s*cos(phi)
1469        v = s*sin(phi)
1470
1471        # Compute wind stress
1472        S = const * sqrt(u**2 + v**2)
1473        xmom_update[k] += S*u
1474        ymom_update[k] += S*v
1475
1476
1477
1478class Rainfall:
1479    """Class Rainfall - general 'rain over entire domain' forcing term.
1480   
1481    Used for implementing Rainfall over the entire domain.
1482       
1483        Current Limited to only One Gauge..
1484       
1485        Need to add Spatial Varying Capability
1486        (This module came from copying and amending the Inflow Code)
1487   
1488    Rainfall(rain)
1489       
1490    rain [mm/s]:  Total rain rate over the specified domain. 
1491                  NOTE: Raingauge Data needs to reflect the time step.
1492                  IE: if Gauge is mm read at a time step, then the input
1493                  here is as mm/(timeStep) so 10mm in 5minutes becomes
1494                  10/(5x60) = 0.0333mm/s.
1495       
1496       
1497                  This parameter can be either a constant or a
1498                  function of time. Positive values indicate inflow,
1499                  negative values indicate outflow.
1500                  (and be used for Infiltration - Write Seperate Module)
1501                  The specified flow will be divided by the area of
1502                  the inflow region and then applied to update the
1503                  quantity in question.
1504
1505    polygon: Specifies a polygon to restrict the rainfall.
1506   
1507    Examples
1508    How to put them in a run File...
1509       
1510    #--------------------------------------------------------------------------
1511    # Setup specialised forcing terms
1512    #--------------------------------------------------------------------------
1513    # This is the new element implemented by Ole and Rudy to allow direct
1514    # input of Inflow in mm/s
1515
1516    catchmentrainfall = Rainfall(rain=file_function('Q100_2hr_Rain.tms')) 
1517                        # Note need path to File in String.
1518                        # Else assumed in same directory
1519
1520    domain.forcing_terms.append(catchmentrainfall)
1521    """
1522
1523    # FIXME (OLE): Add a polygon as an alternative.
1524    # FIXME (AnyOne) : Add various methods to allow spatial variations
1525    # FIXME (OLE): Generalise to all quantities
1526
1527   
1528    def __init__(self,
1529                 rain=0.0,
1530                 quantity_name='stage',
1531                 polygon=None):
1532
1533        self.rain = rain
1534        self.quantity_name = quantity_name
1535        self.polygon = polygon
1536   
1537    def __call__(self, domain):
1538       
1539        # FIXME(Ole): Move this to top of file eventually
1540        from anuga.utilities.polygon import inside_polygon       
1541
1542        # Update rainfall
1543        if callable(self.rain):
1544            rain = self.rain(domain.get_time())
1545        else:
1546            rain = self.rain
1547
1548        # Now rain is a number
1549        # Converting mm/s to m/s to apply in ANUGA
1550        # 1mm of rainfall is equivalent to 1 litre /m2
1551        # Flow is expressed as m3/s converted to a stage height in (m)
1552       
1553        # Note 1m3 = 1x10^9mm3 (mls)
1554        # or is that m3 to Litres ??? Check this how is it applied !!!
1555        rain_fall = rain/1000
1556
1557        # Find indices for coordinates inside polygon if specified
1558        indices = None
1559        if self.polygon is not None:
1560            points = domain.get_centroid_coordinates(absolute=True)
1561            indices = inside_polygon(points, self.polygon)
1562           
1563
1564        # Assign rainfall value to explicit update vector to be used as
1565        # a forcing term
1566        quantity = domain.quantities[self.quantity_name].explicit_update
1567        if indices is None:
1568            quantity[:] += rain_fall
1569        else:
1570            # Brute force assignment of restricted rainfall
1571            for i in indices:
1572                quantity[i] += rain_fall
1573
1574       
1575
1576
1577
1578
1579class Inflow:
1580    """Class Inflow - general 'rain and drain' forcing term.
1581   
1582    Useful for implementing flows in and out of the domain.
1583   
1584    Inflow(center, radius, flow)
1585   
1586    center [m]: Coordinates at center of flow point
1587    radius [m]: Size of circular area
1588    flow [m^3/s]: Total flow rate over the specified area. 
1589                  This parameter can be either a constant or a
1590                  function of time. Positive values indicate inflow,
1591                  negative values indicate outflow.
1592                  The specified flow will be divided by the area of
1593                  the inflow region and then applied to update the
1594                  quantity in question.
1595   
1596    Examples
1597
1598    # Constant drain at 0.003 m^3/s.
1599    # The outflow area is 0.07**2*pi=0.0154 m^2
1600    # This corresponds to a rate of change of 0.003/0.0154 = 0.2 m/s
1601    #                                     
1602    Inflow((0.7, 0.4), 0.07, -0.003)
1603
1604
1605    # Tap turning up to a maximum inflow of 0.0142 m^3/s.
1606    # The inflow area is 0.03**2*pi = 0.00283 m^2
1607    # This corresponds to a rate of change of 0.0142/0.00283 = 5 m/s     
1608    # over the specified area
1609    Inflow((0.5, 0.5), 0.03, lambda t: min(0.01*t, 0.0142))
1610
1611    #--------------------------------------------------------------------------
1612    # Setup specialised forcing terms
1613    #--------------------------------------------------------------------------
1614    # This is the new element implemented by Ole to allow direct input
1615    # of Inflow in m^3/s
1616
1617    hydrograph = Inflow(center=(320, 300), radius=10,
1618                        flow=file_function('Q/QPMF_Rot_Sub13.tms'))
1619
1620    domain.forcing_terms.append(hydrograph)
1621   
1622    """
1623
1624    # FIXME (OLE): Add a polygon as an alternative.
1625    # FIXME (OLE): Generalise to all quantities
1626
1627    def __init__(self,
1628                 center=None, radius=None,
1629                 flow=0.0,
1630                 quantity_name='stage',
1631                 verbose=False):
1632
1633        from math import pi
1634
1635       
1636                 
1637        if center is not None and radius is not None:
1638            assert len(center) == 2
1639        else:
1640            msg = 'Both center and radius must be specified'
1641            raise Exception, msg
1642   
1643        self.center = center
1644        self.radius = radius
1645        self.area = radius**2*pi
1646        self.flow = flow
1647        self.quantity_name = quantity_name
1648        self.verbose = verbose
1649   
1650    def __call__(self, domain):
1651
1652        # Determine indices in flow area
1653        if not hasattr(self, 'indices'):
1654            center = self.center
1655            radius = self.radius
1656           
1657            N = len(domain)   
1658            self.indices = []
1659            coordinates = domain.get_centroid_coordinates()     
1660            for k in range(N):
1661                x, y = coordinates[k,:] # Centroid
1662                if ((x-center[0])**2+(y-center[1])**2) < radius**2:
1663                    self.indices.append(k)   
1664
1665        # Update inflow
1666        if callable(self.flow):
1667            flow = self.flow(domain.get_time())
1668        else:
1669            flow = self.flow
1670
1671        # Now flow is a number
1672        if self.verbose is True:
1673            print 'Inflow at time = %.2f = %f' %(domain.get_time(),
1674                                                 flow)
1675       
1676        quantity = domain.quantities[self.quantity_name].explicit_update
1677        for k in self.indices:
1678            quantity[k] += flow/self.area                       
1679
1680
1681#------------------
1682# Initialise module
1683#------------------
1684
1685
1686from anuga.utilities import compile
1687if compile.can_use_C_extension('shallow_water_ext.c'):
1688    # Underlying C implementations can be accessed
1689
1690    from shallow_water_ext import rotate, assign_windfield_values
1691else:
1692    msg = 'C implementations could not be accessed by %s.\n ' %__file__
1693    msg += 'Make sure compile_all.py has been run as described in '
1694    msg += 'the ANUGA installation guide.'
1695    raise Exception, msg
1696
1697
1698# Optimisation with psyco
1699from anuga.config import use_psyco
1700if use_psyco:
1701    try:
1702        import psyco
1703    except:
1704        import os
1705        if os.name == 'posix' and os.uname()[4] == 'x86_64':
1706            pass
1707            #Psyco isn't supported on 64 bit systems, but it doesn't matter
1708        else:
1709            msg = 'WARNING: psyco (speedup) could not import'+\
1710                  ', you may want to consider installing it'
1711            print msg
1712    else:
1713        psyco.bind(Domain.distribute_to_vertices_and_edges)
1714        psyco.bind(Domain.compute_fluxes)
1715
1716if __name__ == "__main__":
1717    pass
1718
1719
Note: See TracBrowser for help on using the repository browser.