source: branches/numpy/anuga/shallow_water/shallow_water_domain.py @ 6410

Last change on this file since 6410 was 6410, checked in by rwilson, 16 years ago

numpy changes.

  • Property svn:keywords set to LastChangedDate LastChangedRevision LastChangedBy HeadURL Id
File size: 78.2 KB
RevLine 
[6410]1"""
2Finite-volume computations of the shallow water wave equation.
[4004]3
[4005]4Title: ANGUA shallow_water_domain - 2D triangular domains for finite-volume
5       computations of the shallow water wave equation.
[3804]6
7
[5186]8Author: Ole Nielsen, Ole.Nielsen@ga.gov.au
9        Stephen Roberts, Stephen.Roberts@anu.edu.au
10        Duncan Gray, Duncan.Gray@ga.gov.au
[4004]11
12CreationDate: 2004
13
14Description:
[4005]15    This module contains a specialisation of class Domain from
16    module domain.py consisting of methods specific to the
17    Shallow Water Wave Equation
[4004]18
[4005]19    U_t + E_x + G_y = S
[3804]20
[4005]21    where
[3804]22
[4005]23    U = [w, uh, vh]
24    E = [uh, u^2h + gh^2/2, uvh]
25    G = [vh, uvh, v^2h + gh^2/2]
26    S represents source terms forcing the system
27    (e.g. gravity, friction, wind stress, ...)
[3804]28
[4005]29    and _t, _x, _y denote the derivative with respect to t, x and y
30    respectively.
[3804]31
32
[4005]33    The quantities are
[3804]34
[4005]35    symbol    variable name    explanation
36    x         x                horizontal distance from origin [m]
37    y         y                vertical distance from origin [m]
38    z         elevation        elevation of bed on which flow is modelled [m]
39    h         height           water height above z [m]
40    w         stage            absolute water level, w = z+h [m]
41    u                          speed in the x direction [m/s]
42    v                          speed in the y direction [m/s]
43    uh        xmomentum        momentum in the x direction [m^2/s]
44    vh        ymomentum        momentum in the y direction [m^2/s]
[3804]45
[4005]46    eta                        mannings friction coefficient [to appear]
47    nu                         wind stress coefficient [to appear]
[3804]48
[4005]49    The conserved quantities are w, uh, vh
[3804]50
[4004]51Reference:
[4005]52    Catastrophic Collapse of Water Supply Reservoirs in Urban Areas,
53    Christopher Zoppou and Stephen Roberts,
54    Journal of Hydraulic Engineering, vol. 127, No. 7 July 1999
[3804]55
[6410]56    Hydrodynamic modelling of coastal inundation.
[4005]57    Nielsen, O., S. Roberts, D. Gray, A. McPherson and A. Hitchman
58    In Zerger, A. and Argent, R.M. (eds) MODSIM 2005 International Congress on
59    Modelling and Simulation. Modelling and Simulation Society of Australia and
60    New Zealand, December 2005, pp. 518-523. ISBN: 0-9758400-2-9.
61    http://www.mssanz.org.au/modsim05/papers/nielsen.pdf
[3804]62
63
[4005]64SeeAlso:
65    TRAC administration of ANUGA (User Manuals etc) at
66    https://datamining.anu.edu.au/anuga and Subversion repository at
67    $HeadURL: branches/numpy/anuga/shallow_water/shallow_water_domain.py $
[3804]68
[4004]69Constraints: See GPL license in the user guide
70Version: 1.0 ($Revision: 6410 $)
71ModifiedBy:
[4005]72    $Author: rwilson $
73    $Date: 2009-02-24 22:37:22 +0000 (Tue, 24 Feb 2009) $
[3804]74"""
75
[4769]76# Subversion keywords:
[3804]77#
[4769]78# $LastChangedDate: 2009-02-24 22:37:22 +0000 (Tue, 24 Feb 2009) $
79# $LastChangedRevision: 6410 $
80# $LastChangedBy: rwilson $
[3804]81
[6410]82
[6304]83import numpy as num
[3804]84
[5729]85from anuga.abstract_2d_finite_volumes.neighbour_mesh import segment_midpoints
[3804]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
[6178]97from anuga.pmesh.mesh_interface import create_mesh_from_regions
[5294]98from anuga.utilities.numerical_tools import gradient, mean, ensure_numeric
[5730]99from anuga.geospatial_data.geospatial_data import ensure_geospatial
100
[3804]101from anuga.config import minimum_storable_height
102from anuga.config import minimum_allowed_height, maximum_allowed_speed
[5442]103from anuga.config import g, epsilon, beta_w, beta_w_dry,\
[4631]104     beta_uh, beta_uh_dry, beta_vh, beta_vh_dry, tight_slope_limiters
[3876]105from anuga.config import alpha_balance
[4685]106from anuga.config import optimise_dry_cells
[5162]107from anuga.config import optimised_gradient_limiter
[5175]108from anuga.config import use_edge_limiter
[5290]109from anuga.config import use_centroid_velocities
[6086]110from anuga.config import netcdf_mode_r, netcdf_mode_w, netcdf_mode_a
[3804]111
[6410]112from anuga.fit_interpolate.interpolate import Modeltime_too_late, \
113                                              Modeltime_too_early
[5290]114
[6410]115from anuga.utilities.polygon import inside_polygon, polygon_area, \
116                                    is_inside_polygon
[5294]117
[5873]118from types import IntType, FloatType
[5874]119from warnings import warn
[5294]120
121
[6410]122################################################################################
123# Shallow water domain
124################################################################################
[6178]125
[6410]126##
127# @brief Class for a shallow water domain.
[3804]128class Domain(Generic_Domain):
129
[4454]130    conserved_quantities = ['stage', 'xmomentum', 'ymomentum']
131    other_quantities = ['elevation', 'friction']
[6410]132
133    ##
134    # @brief Instantiate a shallow water domain.
135    # @param coordinates
136    # @param vertices
137    # @param boundary
138    # @param tagged_elements
139    # @param geo_reference
140    # @param use_inscribed_circle
141    # @param mesh_filename
142    # @param use_cache
143    # @param verbose
144    # @param full_send_dict
145    # @param ghost_recv_dict
146    # @param processor
147    # @param numproc
148    # @param number_of_full_nodes
149    # @param number_of_full_triangles
[3804]150    def __init__(self,
151                 coordinates=None,
152                 vertices=None,
153                 boundary=None,
154                 tagged_elements=None,
155                 geo_reference=None,
156                 use_inscribed_circle=False,
157                 mesh_filename=None,
158                 use_cache=False,
159                 verbose=False,
160                 full_send_dict=None,
161                 ghost_recv_dict=None,
162                 processor=0,
[3926]163                 numproc=1,
[3928]164                 number_of_full_nodes=None,
165                 number_of_full_triangles=None):
[3804]166
167        other_quantities = ['elevation', 'friction']
168        Generic_Domain.__init__(self,
169                                coordinates,
170                                vertices,
171                                boundary,
[4454]172                                Domain.conserved_quantities,
173                                Domain.other_quantities,
[3804]174                                tagged_elements,
175                                geo_reference,
176                                use_inscribed_circle,
177                                mesh_filename,
178                                use_cache,
179                                verbose,
180                                full_send_dict,
181                                ghost_recv_dict,
182                                processor,
[3926]183                                numproc,
184                                number_of_full_nodes=number_of_full_nodes,
[6410]185                                number_of_full_triangles=number_of_full_triangles)
[3804]186
[6410]187        self.set_minimum_allowed_height(minimum_allowed_height)
[4769]188
[3804]189        self.maximum_allowed_speed = maximum_allowed_speed
190        self.g = g
[6410]191        self.beta_w = beta_w
192        self.beta_w_dry = beta_w_dry
193        self.beta_uh = beta_uh
[3804]194        self.beta_uh_dry = beta_uh_dry
[6410]195        self.beta_vh = beta_vh
[3804]196        self.beta_vh_dry = beta_vh_dry
[3876]197        self.alpha_balance = alpha_balance
[3804]198
[4631]199        self.tight_slope_limiters = tight_slope_limiters
[4685]200        self.optimise_dry_cells = optimise_dry_cells
[4239]201
[4769]202        self.forcing_terms.append(manning_friction_implicit)
[3804]203        self.forcing_terms.append(gravity)
204
[4701]205        # Stored output
[3804]206        self.store = True
207        self.format = 'sww'
208        self.set_store_vertices_uniquely(False)
209        self.minimum_storable_height = minimum_storable_height
[6410]210        self.quantities_to_be_stored = ['stage', 'xmomentum', 'ymomentum']
[5162]211
212        # Limiters
[5175]213        self.use_edge_limiter = use_edge_limiter
[5162]214        self.optimised_gradient_limiter = optimised_gradient_limiter
[5290]215        self.use_centroid_velocities = use_centroid_velocities
[3804]216
[6410]217    ##
218    # @brief
219    # @param beta
[3848]220    def set_all_limiters(self, beta):
[6410]221        """Shorthand to assign one constant value [0,1] to all limiters.
[3847]222        0 Corresponds to first order, where as larger values make use of
[6410]223        the second order scheme.
[3847]224        """
[3804]225
[6410]226        self.beta_w = beta
227        self.beta_w_dry = beta
[5162]228        self.quantities['stage'].beta = beta
[6410]229
230        self.beta_uh = beta
[3847]231        self.beta_uh_dry = beta
[5162]232        self.quantities['xmomentum'].beta = beta
[6410]233
234        self.beta_vh = beta
[3847]235        self.beta_vh_dry = beta
[5162]236        self.quantities['ymomentum'].beta = beta
[3847]237
[6410]238    ##
239    # @brief
240    # @param flag
241    # @param reduction
[3804]242    def set_store_vertices_uniquely(self, flag, reduction=None):
243        """Decide whether vertex values should be stored uniquely as
244        computed in the model or whether they should be reduced to one
245        value per vertex using self.reduction.
246        """
[3954]247
248        # FIXME (Ole): how about using the word continuous vertex values?
[3804]249        self.smooth = not flag
250
[4733]251        # Reduction operation for get_vertex_values
[3804]252        if reduction is None:
253            self.reduction = mean
254            #self.reduction = min  #Looks better near steep slopes
255
[6410]256    ##
257    # @brief Set the minimum depth that will be written to an SWW file.
258    # @param minimum_storable_height The minimum stored height (in m).
[3804]259    def set_minimum_storable_height(self, minimum_storable_height):
[6410]260        """Set the minimum depth that will be recognised when writing
[3804]261        to an sww file. This is useful for removing thin water layers
262        that seems to be caused by friction creep.
263
264        The minimum allowed sww depth is in meters.
265        """
[6410]266
[3804]267        self.minimum_storable_height = minimum_storable_height
[4258]268
[6410]269    ##
270    # @brief
271    # @param minimum_allowed_height
[4258]272    def set_minimum_allowed_height(self, minimum_allowed_height):
[6410]273        """Set the minimum depth recognised in the numerical scheme.
[4258]274
275        The minimum allowed depth is in meters.
276
277        The parameter H0 (Minimal height for flux computation)
278        is also set by this function
279        """
[4438]280
281        #FIXME (Ole): rename H0 to minimum_allowed_height_in_flux_computation
[4701]282
283        #FIXME (Ole): Maybe use histogram to identify isolated extreme speeds
284        #and deal with them adaptively similarly to how we used to use 1 order
285        #steps to recover.
[6410]286
[4258]287        self.minimum_allowed_height = minimum_allowed_height
[6410]288        self.H0 = minimum_allowed_height
[3804]289
[6410]290    ##
291    # @brief
292    # @param maximum_allowed_speed
[3804]293    def set_maximum_allowed_speed(self, maximum_allowed_speed):
[6410]294        """Set the maximum particle speed that is allowed in water
[3804]295        shallower than minimum_allowed_height. This is useful for
296        controlling speeds in very thin layers of water and at the same time
297        allow some movement avoiding pooling of water.
[6410]298        """
[3804]299
300        self.maximum_allowed_speed = maximum_allowed_speed
301
[6410]302    ##
303    # @brief
304    # @param points_file_block_line_size
305    def set_points_file_block_line_size(self, points_file_block_line_size):
306        """Set the minimum depth that will be recognised when writing
[4254]307        to an sww file. This is useful for removing thin water layers
308        that seems to be caused by friction creep.
309
310        The minimum allowed sww depth is in meters.
311        """
312        self.points_file_block_line_size = points_file_block_line_size
[6410]313
314    ##
315    # @brief Set the quantities that will be written to an SWW file.
316    # @param q The quantities to be written.
317    # @note Param 'q' may be None, single quantity or list of quantity strings.
318    # @note If 'q' is None, no quantities will be stored in the SWW file.
[3804]319    def set_quantities_to_be_stored(self, q):
320        """Specify which quantities will be stored in the sww file.
321
322        q must be either:
323          - the name of a quantity
324          - a list of quantity names
325          - None
326
327        In the two first cases, the named quantities will be stored at
328        each yieldstep (This is in addition to the quantities elevation
329        and friction)
[6410]330
[3804]331        If q is None, storage will be switched off altogether.
332        """
333
334        if q is None:
335            self.quantities_to_be_stored = []
336            self.store = False
337            return
338
339        if isinstance(q, basestring):
340            q = [q] # Turn argument into a list
341
[4701]342        # Check correcness
[3804]343        for quantity_name in q:
[6410]344            msg = ('Quantity %s is not a valid conserved quantity'
345                   % quantity_name)
[3804]346            assert quantity_name in self.conserved_quantities, msg
347
348        self.quantities_to_be_stored = q
349
[6410]350    ##
351    # @brief
352    # @param indices
[3804]353    def get_wet_elements(self, indices=None):
354        """Return indices for elements where h > minimum_allowed_height
355
356        Optional argument:
357            indices is the set of element ids that the operation applies to.
358
359        Usage:
360            indices = get_wet_elements()
361
[6410]362        Note, centroid values are used for this operation
[3804]363        """
364
365        # Water depth below which it is considered to be 0 in the model
366        # FIXME (Ole): Allow this to be specified as a keyword argument as well
367        from anuga.config import minimum_allowed_height
368
369        elevation = self.get_quantity('elevation').\
[6410]370                        get_values(location='centroids', indices=indices)
371        stage = self.get_quantity('stage').\
[3804]372                    get_values(location='centroids', indices=indices)
373        depth = stage - elevation
374
375        # Select indices for which depth > 0
[6157]376        wet_indices = num.compress(depth > minimum_allowed_height,
377                                   num.arange(len(depth)))
[6410]378        return wet_indices
[3804]379
[6410]380    ##
381    # @brief
382    # @param indices
[3804]383    def get_maximum_inundation_elevation(self, indices=None):
384        """Return highest elevation where h > 0
385
386        Optional argument:
387            indices is the set of element ids that the operation applies to.
388
389        Usage:
390            q = get_maximum_inundation_elevation()
391
[6410]392        Note, centroid values are used for this operation
[3804]393        """
394
395        wet_elements = self.get_wet_elements(indices)
396        return self.get_quantity('elevation').\
[6410]397                   get_maximum_value(indices=wet_elements)
[3804]398
[6410]399    ##
400    # @brief
401    # @param indices
[3804]402    def get_maximum_inundation_location(self, indices=None):
[4554]403        """Return location of highest elevation where h > 0
[3804]404
405        Optional argument:
406            indices is the set of element ids that the operation applies to.
407
408        Usage:
[4554]409            q = get_maximum_inundation_location()
[3804]410
[6410]411        Note, centroid values are used for this operation
[3804]412        """
413
414        wet_elements = self.get_wet_elements(indices)
415        return self.get_quantity('elevation').\
[6410]416                   get_maximum_location(indices=wet_elements)
417
418    ##
419    # @brief Get the total flow through an arbitrary poly line.
420    # @param polyline Representation of desired cross section.
421    # @param verbose True if this method is to be verbose.
422    # @note 'polyline' may contain multiple sections allowing complex shapes.
423    # @note Assume absolute UTM coordinates.
424    def get_flow_through_cross_section(self, polyline, verbose=False):
425        """Get the total flow through an arbitrary poly line.
426
427        This is a run-time equivalent of the function with same name
[5736]428        in data_manager.py
[6410]429
[5729]430        Input:
[6410]431            polyline: Representation of desired cross section - it may contain
432                      multiple sections allowing for complex shapes. Assume
[5729]433                      absolute UTM coordinates.
[6410]434                      Format [[x0, y0], [x1, y1], ...]
435
436        Output:
[5729]437            Q: Total flow [m^3/s] across given segments.
[6410]438        """
439
[5729]440        # Find all intersections and associated triangles.
[6410]441        segments = self.get_intersecting_segments(polyline, use_cache=True,
[5862]442                                                  verbose=verbose)
[3804]443
[5729]444        # Get midpoints
[6410]445        midpoints = segment_midpoints(segments)
446
[5730]447        # Make midpoints Geospatial instances
[6410]448        midpoints = ensure_geospatial(midpoints, self.geo_reference)
449
450        # Compute flow
[5729]451        if verbose: print 'Computing flow through specified cross section'
[6410]452
[5729]453        # Get interpolated values
454        xmomentum = self.get_quantity('xmomentum')
[6410]455        ymomentum = self.get_quantity('ymomentum')
456
457        uh = xmomentum.get_values(interpolation_points=midpoints,
458                                  use_cache=True)
459        vh = ymomentum.get_values(interpolation_points=midpoints,
460                                  use_cache=True)
461
[5729]462        # Compute and sum flows across each segment
[6410]463        total_flow = 0
[5729]464        for i in range(len(uh)):
[6410]465            # Inner product of momentum vector with segment normal [m^2/s]
[5729]466            normal = segments[i].normal
[6410]467            normal_momentum = uh[i]*normal[0] + vh[i]*normal[1]
468
[5729]469            # Flow across this segment [m^3/s]
470            segment_flow = normal_momentum*segments[i].length
471
472            # Accumulate
473            total_flow += segment_flow
[6410]474
[5729]475        return total_flow
[5736]476
[6410]477    ##
478    # @brief
479    # @param polyline Representation of desired cross section.
480    # @param kind Select energy type to compute ('specific' or 'total').
481    # @param verbose True if this method is to be verbose.
482    # @note 'polyline' may contain multiple sections allowing complex shapes.
483    # @note Assume absolute UTM coordinates.
[5736]484    def get_energy_through_cross_section(self, polyline,
[5773]485                                         kind='total',
[6410]486                                         verbose=False):
[5736]487        """Obtain average energy head [m] across specified cross section.
[5729]488
[5736]489        Inputs:
[6410]490            polyline: Representation of desired cross section - it may contain
491                      multiple sections allowing for complex shapes. Assume
[5736]492                      absolute UTM coordinates.
493                      Format [[x0, y0], [x1, y1], ...]
[6410]494            kind:     Select which energy to compute.
[5736]495                      Options are 'specific' and 'total' (default)
496
497        Output:
498            E: Average energy [m] across given segments for all stored times.
499
[6410]500        The average velocity is computed for each triangle intersected by
501        the polyline and averaged weighted by segment lengths.
[5736]502
[6410]503        The typical usage of this function would be to get average energy of
504        flow in a channel, and the polyline would then be a cross section
[5736]505        perpendicular to the flow.
506
[6410]507        #FIXME (Ole) - need name for this energy reflecting that its dimension
[5736]508        is [m].
509        """
510
[6410]511        from anuga.config import g, epsilon, velocity_protection as h0
512
[5736]513        # Find all intersections and associated triangles.
[6410]514        segments = self.get_intersecting_segments(polyline, use_cache=True,
[5862]515                                                  verbose=verbose)
[5736]516
517        # Get midpoints
[6410]518        midpoints = segment_midpoints(segments)
519
[5736]520        # Make midpoints Geospatial instances
[6410]521        midpoints = ensure_geospatial(midpoints, self.geo_reference)
522
523        # Compute energy
524        if verbose: print 'Computing %s energy' %kind
525
[5736]526        # Get interpolated values
[6410]527        stage = self.get_quantity('stage')
528        elevation = self.get_quantity('elevation')
[5736]529        xmomentum = self.get_quantity('xmomentum')
[6410]530        ymomentum = self.get_quantity('ymomentum')
[5736]531
[5866]532        w = stage.get_values(interpolation_points=midpoints, use_cache=True)
[6410]533        z = elevation.get_values(interpolation_points=midpoints, use_cache=True)
534        uh = xmomentum.get_values(interpolation_points=midpoints,
535                                  use_cache=True)
536        vh = ymomentum.get_values(interpolation_points=midpoints,
537                                  use_cache=True)
538        h = w-z                # Depth
539
[5736]540        # Compute total length of polyline for use with weighted averages
541        total_line_length = 0.0
542        for segment in segments:
543            total_line_length += segment.length
[6410]544
[5736]545        # Compute and sum flows across each segment
[6410]546        average_energy = 0.0
[5736]547        for i in range(len(w)):
548            # Average velocity across this segment
549            if h[i] > epsilon:
550                # Use protection against degenerate velocities
551                u = uh[i]/(h[i] + h0/h[i])
552                v = vh[i]/(h[i] + h0/h[i])
553            else:
554                u = v = 0.0
[6410]555
556            speed_squared = u*u + v*v
[5736]557            kinetic_energy = 0.5*speed_squared/g
[6410]558
[5736]559            if kind == 'specific':
560                segment_energy = h[i] + kinetic_energy
561            elif kind == 'total':
[6410]562                segment_energy = w[i] + kinetic_energy
[5736]563            else:
564                msg = 'Energy kind must be either "specific" or "total".'
565                msg += ' I got %s' %kind
566
567            # Add to weighted average
568            weigth = segments[i].length/total_line_length
569            average_energy += segment_energy*weigth
[6410]570
[5736]571        return average_energy
572
[6410]573    ##
574    # @brief Run integrity checks on shallow water domain.
[3804]575    def check_integrity(self):
[6410]576        Generic_Domain.check_integrity(self)    # super integrity check
[3804]577
578        #Check that we are solving the shallow water wave equation
579        msg = 'First conserved quantity must be "stage"'
580        assert self.conserved_quantities[0] == 'stage', msg
581        msg = 'Second conserved quantity must be "xmomentum"'
582        assert self.conserved_quantities[1] == 'xmomentum', msg
583        msg = 'Third conserved quantity must be "ymomentum"'
584        assert self.conserved_quantities[2] == 'ymomentum', msg
585
[6410]586    ##
587    # @brief
[3804]588    def extrapolate_second_order_sw(self):
[6410]589        #Call correct module function (either from this module or C-extension)
[3804]590        extrapolate_second_order_sw(self)
591
[6410]592    ##
593    # @brief
[3804]594    def compute_fluxes(self):
[6410]595        #Call correct module function (either from this module or C-extension)
[3804]596        compute_fluxes(self)
597
[6410]598    ##
599    # @brief
[3804]600    def distribute_to_vertices_and_edges(self):
[4769]601        # Call correct module function
[5176]602        if self.use_edge_limiter:
[6410]603            distribute_using_edge_limiter(self)
[5175]604        else:
[5306]605            distribute_using_vertex_limiter(self)
[3804]606
[6410]607    ##
608    # @brief Evolve the model by one step.
609    # @param yieldstep
610    # @param finaltime
611    # @param duration
612    # @param skip_initial_step
[3804]613    def evolve(self,
[6410]614               yieldstep=None,
615               finaltime=None,
616               duration=None,
617               skip_initial_step=False):
618        """Specialisation of basic evolve method from parent class"""
[3804]619
[4769]620        # Call check integrity here rather than from user scripts
621        # self.check_integrity()
[3804]622
[6410]623        msg = 'Attribute self.beta_w must be in the interval [0, 2]'
[5162]624        assert 0 <= self.beta_w <= 2.0, msg
[3804]625
[4769]626        # Initial update of vertex and edge values before any STORAGE
[6410]627        # and or visualisation.
[4769]628        # This is done again in the initialisation of the Generic_Domain
[6410]629        # evolve loop but we do it here to ensure the values are ok for storage.
[3804]630        self.distribute_to_vertices_and_edges()
631
632        if self.store is True and self.time == 0.0:
633            self.initialise_storage()
634        else:
635            pass
[4769]636            # print 'Results will not be stored.'
637            # print 'To store results set domain.store = True'
638            # FIXME: Diagnostic output should be controlled by
[3804]639            # a 'verbose' flag living in domain (or in a parent class)
640
[4769]641        # Call basic machinery from parent class
[6410]642        for t in Generic_Domain.evolve(self, yieldstep=yieldstep,
643                                       finaltime=finaltime, duration=duration,
[3804]644                                       skip_initial_step=skip_initial_step):
[4769]645            # Store model data, e.g. for subsequent visualisation
[3804]646            if self.store is True:
647                self.store_timestep(self.quantities_to_be_stored)
648
[4769]649            # FIXME: Could maybe be taken from specified list
650            # of 'store every step' quantities
[3804]651
[4769]652            # Pass control on to outer loop for more specific actions
[3804]653            yield(t)
654
[6410]655    ##
656    # @brief
[3804]657    def initialise_storage(self):
658        """Create and initialise self.writer object for storing data.
659        Also, save x,y and bed elevation
660        """
661
662        from anuga.shallow_water.data_manager import get_dataobject
663
[4769]664        # Initialise writer
[6086]665        self.writer = get_dataobject(self, mode=netcdf_mode_w)
[3804]666
[4769]667        # Store vertices and connectivity
[3804]668        self.writer.store_connectivity()
669
[6410]670    ##
671    # @brief
672    # @param name
[3804]673    def store_timestep(self, name):
674        """Store named quantity and time.
675
676        Precondition:
677           self.write has been initialised
678        """
[6410]679
[3804]680        self.writer.store_timestep(name)
681
[6410]682    ##
683    # @brief Get time stepping statistics string for printing.
684    # @param track_speeds
685    # @param triangle_id
[4836]686    def timestepping_statistics(self,
687                                track_speeds=False,
[6410]688                                triangle_id=None):
[4827]689        """Return string with time stepping statistics for printing or logging
[3804]690
[4827]691        Optional boolean keyword track_speeds decides whether to report
692        location of smallest timestep as well as a histogram and percentile
693        report.
694        """
695
[6410]696        from anuga.config import epsilon, g
[4827]697
698        # Call basic machinery from parent class
[6410]699        msg = Generic_Domain.timestepping_statistics(self, track_speeds,
[4836]700                                                     triangle_id)
[4827]701
702        if track_speeds is True:
703            # qwidth determines the text field used for quantities
704            qwidth = self.qwidth
[6410]705
[4836]706            # Selected triangle
[4827]707            k = self.k
708
709            # Report some derived quantities at vertices, edges and centroid
710            # specific to the shallow water wave equation
711            z = self.quantities['elevation']
[6410]712            w = self.quantities['stage']
[4827]713
714            Vw = w.get_values(location='vertices', indices=[k])[0]
715            Ew = w.get_values(location='edges', indices=[k])[0]
716            Cw = w.get_values(location='centroids', indices=[k])
717
718            Vz = z.get_values(location='vertices', indices=[k])[0]
719            Ez = z.get_values(location='edges', indices=[k])[0]
[6410]720            Cz = z.get_values(location='centroids', indices=[k])
[4827]721
722            name = 'depth'
723            Vh = Vw-Vz
724            Eh = Ew-Ez
725            Ch = Cw-Cz
[6410]726
[4827]727            s  = '    %s: vertex_values =  %.4f,\t %.4f,\t %.4f\n'\
728                 %(name.ljust(qwidth), Vh[0], Vh[1], Vh[2])
[6410]729
[4827]730            s += '    %s: edge_values =    %.4f,\t %.4f,\t %.4f\n'\
731                 %(name.ljust(qwidth), Eh[0], Eh[1], Eh[2])
[6410]732
[4827]733            s += '    %s: centroid_value = %.4f\n'\
[6410]734                 %(name.ljust(qwidth), Ch[0])
735
[4827]736            msg += s
737
738            uh = self.quantities['xmomentum']
739            vh = self.quantities['ymomentum']
740
741            Vuh = uh.get_values(location='vertices', indices=[k])[0]
742            Euh = uh.get_values(location='edges', indices=[k])[0]
743            Cuh = uh.get_values(location='centroids', indices=[k])
[6410]744
[4827]745            Vvh = vh.get_values(location='vertices', indices=[k])[0]
746            Evh = vh.get_values(location='edges', indices=[k])[0]
747            Cvh = vh.get_values(location='centroids', indices=[k])
748
749            # Speeds in each direction
750            Vu = Vuh/(Vh + epsilon)
751            Eu = Euh/(Eh + epsilon)
[6410]752            Cu = Cuh/(Ch + epsilon)
[4827]753            name = 'U'
754            s  = '    %s: vertex_values =  %.4f,\t %.4f,\t %.4f\n'\
755                 %(name.ljust(qwidth), Vu[0], Vu[1], Vu[2])
[6410]756
[4827]757            s += '    %s: edge_values =    %.4f,\t %.4f,\t %.4f\n'\
758                 %(name.ljust(qwidth), Eu[0], Eu[1], Eu[2])
[6410]759
[4827]760            s += '    %s: centroid_value = %.4f\n'\
[6410]761                 %(name.ljust(qwidth), Cu[0])
762
[4827]763            msg += s
764
765            Vv = Vvh/(Vh + epsilon)
766            Ev = Evh/(Eh + epsilon)
[6410]767            Cv = Cvh/(Ch + epsilon)
[4827]768            name = 'V'
769            s  = '    %s: vertex_values =  %.4f,\t %.4f,\t %.4f\n'\
770                 %(name.ljust(qwidth), Vv[0], Vv[1], Vv[2])
[6410]771
[4827]772            s += '    %s: edge_values =    %.4f,\t %.4f,\t %.4f\n'\
773                 %(name.ljust(qwidth), Ev[0], Ev[1], Ev[2])
[6410]774
[4827]775            s += '    %s: centroid_value = %.4f\n'\
[6410]776                 %(name.ljust(qwidth), Cv[0])
777
[4827]778            msg += s
779
780            # Froude number in each direction
781            name = 'Froude (x)'
[6157]782            Vfx = Vu/(num.sqrt(g*Vh) + epsilon)
783            Efx = Eu/(num.sqrt(g*Eh) + epsilon)
784            Cfx = Cu/(num.sqrt(g*Ch) + epsilon)
[6410]785
[4827]786            s  = '    %s: vertex_values =  %.4f,\t %.4f,\t %.4f\n'\
787                 %(name.ljust(qwidth), Vfx[0], Vfx[1], Vfx[2])
[6410]788
[4827]789            s += '    %s: edge_values =    %.4f,\t %.4f,\t %.4f\n'\
790                 %(name.ljust(qwidth), Efx[0], Efx[1], Efx[2])
[6410]791
[4827]792            s += '    %s: centroid_value = %.4f\n'\
[6410]793                 %(name.ljust(qwidth), Cfx[0])
794
[4827]795            msg += s
796
797            name = 'Froude (y)'
[6157]798            Vfy = Vv/(num.sqrt(g*Vh) + epsilon)
799            Efy = Ev/(num.sqrt(g*Eh) + epsilon)
800            Cfy = Cv/(num.sqrt(g*Ch) + epsilon)
[6410]801
[4827]802            s  = '    %s: vertex_values =  %.4f,\t %.4f,\t %.4f\n'\
803                 %(name.ljust(qwidth), Vfy[0], Vfy[1], Vfy[2])
[6410]804
[4827]805            s += '    %s: edge_values =    %.4f,\t %.4f,\t %.4f\n'\
806                 %(name.ljust(qwidth), Efy[0], Efy[1], Efy[2])
[6410]807
[4827]808            s += '    %s: centroid_value = %.4f\n'\
[6410]809                 %(name.ljust(qwidth), Cfy[0])
[4827]810
[6410]811            msg += s
[4827]812
813        return msg
814
[6410]815################################################################################
816# End of class Shallow Water Domain
817################################################################################
[3804]818
[4769]819#-----------------
[3804]820# Flux computation
[4769]821#-----------------
[3804]822
[6410]823## @brief Compute fluxes and timestep suitable for all volumes in domain.
824# @param domain The domain to calculate fluxes for.
[3804]825def compute_fluxes(domain):
[6410]826    """Compute fluxes and timestep suitable for all volumes in domain.
[3804]827
828    Compute total flux for each conserved quantity using "flux_function"
829
830    Fluxes across each edge are scaled by edgelengths and summed up
831    Resulting flux is then scaled by area and stored in
832    explicit_update for each of the three conserved quantities
833    stage, xmomentum and ymomentum
834
835    The maximal allowable speed computed by the flux_function for each volume
836    is converted to a timestep that must not be exceeded. The minimum of
837    those is computed as the next overall timestep.
838
839    Post conditions:
840      domain.explicit_update is reset to computed flux values
841      domain.timestep is set to the largest step satisfying all volumes.
[4769]842
843    This wrapper calls the underlying C version of compute fluxes
[3804]844    """
845
846    import sys
[6410]847    from shallow_water_ext import compute_fluxes_ext_central \
848                                  as compute_fluxes_ext
[3804]849
[6410]850    N = len(domain)    # number_of_triangles
[3804]851
[4769]852    # Shortcuts
[3804]853    Stage = domain.quantities['stage']
854    Xmom = domain.quantities['xmomentum']
855    Ymom = domain.quantities['ymomentum']
856    Bed = domain.quantities['elevation']
857
858    timestep = float(sys.maxint)
859
[4769]860    flux_timestep = compute_fluxes_ext(timestep,
861                                       domain.epsilon,
862                                       domain.H0,
863                                       domain.g,
864                                       domain.neighbours,
865                                       domain.neighbour_edges,
866                                       domain.normals,
867                                       domain.edgelengths,
868                                       domain.radii,
869                                       domain.areas,
870                                       domain.tri_full_flag,
871                                       Stage.edge_values,
872                                       Xmom.edge_values,
873                                       Ymom.edge_values,
874                                       Bed.edge_values,
875                                       Stage.boundary_values,
876                                       Xmom.boundary_values,
877                                       Ymom.boundary_values,
878                                       Stage.explicit_update,
879                                       Xmom.explicit_update,
880                                       Ymom.explicit_update,
881                                       domain.already_computed_flux,
882                                       domain.max_speed,
883                                       int(domain.optimise_dry_cells))
[3804]884
[4769]885    domain.flux_timestep = flux_timestep
[3804]886
[6410]887################################################################################
[4769]888# Module functions for gradient limiting
[6410]889################################################################################
[3804]890
[6410]891##
892# @brief Wrapper for C version of extrapolate_second_order_sw.
893# @param domain The domain to operate on.
894# @note MH090605 The following method belongs to the shallow_water domain class
895#       see comments in the corresponding method in shallow_water_ext.c
896def extrapolate_second_order_sw(domain):
897    """Wrapper calling C version of extrapolate_second_order_sw"""
[3804]898
899    import sys
[6410]900    from shallow_water_ext import extrapolate_second_order_sw as extrapol2
[3804]901
[3928]902    N = len(domain) # number_of_triangles
[3804]903
[4710]904    # Shortcuts
[3804]905    Stage = domain.quantities['stage']
906    Xmom = domain.quantities['xmomentum']
907    Ymom = domain.quantities['ymomentum']
908    Elevation = domain.quantities['elevation']
[4710]909
[4769]910    extrapol2(domain,
911              domain.surrogate_neighbours,
912              domain.number_of_boundaries,
913              domain.centroid_coordinates,
914              Stage.centroid_values,
915              Xmom.centroid_values,
916              Ymom.centroid_values,
917              Elevation.centroid_values,
918              domain.vertex_coordinates,
919              Stage.vertex_values,
920              Xmom.vertex_values,
921              Ymom.vertex_values,
922              Elevation.vertex_values,
[5315]923              int(domain.optimise_dry_cells))
[3804]924
[6410]925##
926# @brief Distribution from centroids to vertices specific to the SWW eqn.
927# @param domain The domain to operate on.
[5306]928def distribute_using_vertex_limiter(domain):
[6410]929    """Distribution from centroids to vertices specific to the SWW equation.
[3804]930
931    It will ensure that h (w-z) is always non-negative even in the
932    presence of steep bed-slopes by taking a weighted average between shallow
933    and deep cases.
934
935    In addition, all conserved quantities get distributed as per either a
936    constant (order==1) or a piecewise linear function (order==2).
937
938    FIXME: more explanation about removal of artificial variability etc
939
940    Precondition:
941      All quantities defined at centroids and bed elevation defined at
942      vertices.
943
944    Postcondition
945      Conserved quantities defined at vertices
946    """
947
[4769]948    # Remove very thin layers of water
[3804]949    protect_against_infinitesimal_and_negative_heights(domain)
950
[4769]951    # Extrapolate all conserved quantities
[5162]952    if domain.optimised_gradient_limiter:
[4769]953        # MH090605 if second order,
954        # perform the extrapolation and limiting on
955        # all of the conserved quantities
[3804]956
957        if (domain._order_ == 1):
958            for name in domain.conserved_quantities:
959                Q = domain.quantities[name]
960                Q.extrapolate_first_order()
961        elif domain._order_ == 2:
962            domain.extrapolate_second_order_sw()
963        else:
964            raise 'Unknown order'
965    else:
[4769]966        # Old code:
[3804]967        for name in domain.conserved_quantities:
968            Q = domain.quantities[name]
[4701]969
[3804]970            if domain._order_ == 1:
971                Q.extrapolate_first_order()
972            elif domain._order_ == 2:
[5306]973                Q.extrapolate_second_order_and_limit_by_vertex()
[3804]974            else:
975                raise 'Unknown order'
976
[5290]977    # Take bed elevation into account when water heights are small
[3804]978    balance_deep_and_shallow(domain)
979
[5290]980    # Compute edge values by interpolation
[3804]981    for name in domain.conserved_quantities:
982        Q = domain.quantities[name]
983        Q.interpolate_from_vertices_to_edges()
984
[6410]985##
986# @brief Distribution from centroids to edges specific to the SWW eqn.
987# @param domain The domain to operate on.
[5306]988def distribute_using_edge_limiter(domain):
[6410]989    """Distribution from centroids to edges specific to the SWW eqn.
[5306]990
991    It will ensure that h (w-z) is always non-negative even in the
992    presence of steep bed-slopes by taking a weighted average between shallow
993    and deep cases.
994
995    In addition, all conserved quantities get distributed as per either a
996    constant (order==1) or a piecewise linear function (order==2).
997
998
999    Precondition:
1000      All quantities defined at centroids and bed elevation defined at
1001      vertices.
1002
1003    Postcondition
1004      Conserved quantities defined at vertices
1005    """
1006
1007    # Remove very thin layers of water
1008    protect_against_infinitesimal_and_negative_heights(domain)
1009
1010    for name in domain.conserved_quantities:
1011        Q = domain.quantities[name]
1012        if domain._order_ == 1:
1013            Q.extrapolate_first_order()
1014        elif domain._order_ == 2:
1015            Q.extrapolate_second_order_and_limit_by_edge()
1016        else:
1017            raise 'Unknown order'
1018
1019    balance_deep_and_shallow(domain)
1020
1021    # Compute edge values by interpolation
1022    for name in domain.conserved_quantities:
1023        Q = domain.quantities[name]
1024        Q.interpolate_from_vertices_to_edges()
1025
[6410]1026##
1027# @brief  Protect against infinitesimal heights and associated high velocities.
1028# @param domain The domain to operate on.
[3804]1029def protect_against_infinitesimal_and_negative_heights(domain):
[6410]1030    """Protect against infinitesimal heights and associated high velocities"""
[3804]1031
[6410]1032    from shallow_water_ext import protect
1033
[4769]1034    # Shortcuts
[3804]1035    wc = domain.quantities['stage'].centroid_values
1036    zc = domain.quantities['elevation'].centroid_values
1037    xmomc = domain.quantities['xmomentum'].centroid_values
1038    ymomc = domain.quantities['ymomentum'].centroid_values
1039
1040    protect(domain.minimum_allowed_height, domain.maximum_allowed_speed,
1041            domain.epsilon, wc, zc, xmomc, ymomc)
1042
[6410]1043##
1044# @brief Wrapper for C function balance_deep_and_shallow_c().
1045# @param domain The domain to operate on.
[3804]1046def balance_deep_and_shallow(domain):
1047    """Compute linear combination between stage as computed by
1048    gradient-limiters limiting using w, and stage computed by
1049    gradient-limiters limiting using h (h-limiter).
1050    The former takes precedence when heights are large compared to the
1051    bed slope while the latter takes precedence when heights are
1052    relatively small.  Anything in between is computed as a balanced
1053    linear combination in order to avoid numerical disturbances which
1054    would otherwise appear as a result of hard switching between
1055    modes.
1056
[4769]1057    Wrapper for C implementation
[3804]1058    """
1059
[6410]1060    from shallow_water_ext import balance_deep_and_shallow \
1061                                  as balance_deep_and_shallow_c
[5175]1062
[4733]1063    # Shortcuts
[3804]1064    wc = domain.quantities['stage'].centroid_values
1065    zc = domain.quantities['elevation'].centroid_values
1066    wv = domain.quantities['stage'].vertex_values
1067    zv = domain.quantities['elevation'].vertex_values
1068
[4733]1069    # Momentums at centroids
[3804]1070    xmomc = domain.quantities['xmomentum'].centroid_values
1071    ymomc = domain.quantities['ymomentum'].centroid_values
1072
[4733]1073    # Momentums at vertices
[3804]1074    xmomv = domain.quantities['xmomentum'].vertex_values
1075    ymomv = domain.quantities['ymomentum'].vertex_values
1076
[5442]1077    balance_deep_and_shallow_c(domain,
[6410]1078                               wc, zc, wv, zv, wc,
[5442]1079                               xmomc, ymomc, xmomv, ymomv)
[3804]1080
1081
[6410]1082################################################################################
1083# Boundary conditions - specific to the shallow water wave equation
1084################################################################################
[3804]1085
[6410]1086##
1087# @brief Class for a reflective boundary.
1088# @note Inherits from Boundary.
[3804]1089class Reflective_boundary(Boundary):
1090    """Reflective boundary returns same conserved quantities as
1091    those present in its neighbour volume but reflected.
1092
1093    This class is specific to the shallow water equation as it
1094    works with the momentum quantities assumed to be the second
1095    and third conserved quantities.
1096    """
1097
[6410]1098    ##
1099    # @brief Instantiate a Reflective_boundary.
1100    # @param domain
1101    def __init__(self, domain=None):
[3804]1102        Boundary.__init__(self)
1103
1104        if domain is None:
1105            msg = 'Domain must be specified for reflective boundary'
[6410]1106            raise Exception, msg
[3804]1107
[4769]1108        # Handy shorthands
[6410]1109        self.stage = domain.quantities['stage'].edge_values
1110        self.xmom = domain.quantities['xmomentum'].edge_values
1111        self.ymom = domain.quantities['ymomentum'].edge_values
[3804]1112        self.normals = domain.normals
1113
[6304]1114        self.conserved_quantities = num.zeros(3, num.float)
[3804]1115
[6410]1116    ##
1117    # @brief Return a representation of this instance.
[3804]1118    def __repr__(self):
1119        return 'Reflective_boundary'
1120
[6410]1121    ##
1122    # @brief Calculate reflections (reverse outward momentum).
1123    # @param vol_id
1124    # @param edge_id
[3804]1125    def evaluate(self, vol_id, edge_id):
1126        """Reflective boundaries reverses the outward momentum
1127        of the volume they serve.
1128        """
1129
1130        q = self.conserved_quantities
1131        q[0] = self.stage[vol_id, edge_id]
1132        q[1] = self.xmom[vol_id, edge_id]
1133        q[2] = self.ymom[vol_id, edge_id]
1134
1135        normal = self.normals[vol_id, 2*edge_id:2*edge_id+2]
1136
1137        r = rotate(q, normal, direction = 1)
1138        r[1] = -r[1]
1139        q = rotate(r, normal, direction = -1)
1140
1141        return q
1142
1143
[6410]1144##
1145# @brief Class for a transmissive boundary.
1146# @note Inherits from Boundary.
[6045]1147class Transmissive_momentum_set_stage_boundary(Boundary):
[3804]1148    """Returns same momentum conserved quantities as
1149    those present in its neighbour volume.
1150    Sets stage by specifying a function f of time which may either be a
1151    vector function or a scalar function
1152
1153    Example:
1154
[6410]1155    def waveform(t):
[3804]1156        return sea_level + normalized_amplitude/cosh(t-25)**2
1157
[6045]1158    Bts = Transmissive_momentum_set_stage_boundary(domain, waveform)
[3804]1159
1160    Underlying domain must be specified when boundary is instantiated
1161    """
1162
[6410]1163    ##
1164    # @brief Instantiate a Reflective_boundary.
1165    # @param domain
1166    # @param function
1167    def __init__(self, domain=None, function=None):
[3804]1168        Boundary.__init__(self)
1169
1170        if domain is None:
1171            msg = 'Domain must be specified for this type boundary'
[6410]1172            raise Exception, msg
[3804]1173
1174        if function is None:
1175            msg = 'Function must be specified for this type boundary'
[6410]1176            raise Exception, msg
[3804]1177
[6410]1178        self.domain = domain
[3804]1179        self.function = function
1180
[6410]1181    ##
1182    # @brief Return a representation of this instance.
[3804]1183    def __repr__(self):
[6045]1184        return 'Transmissive_momentum_set_stage_boundary(%s)' %self.domain
[3804]1185
[6410]1186    ##
1187    # @brief Calculate transmissive results.
1188    # @param vol_id
1189    # @param edge_id
[3804]1190    def evaluate(self, vol_id, edge_id):
[6045]1191        """Transmissive momentum set stage boundaries return the edge momentum
[3804]1192        values of the volume they serve.
1193        """
1194
1195        q = self.domain.get_conserved_quantities(vol_id, edge = edge_id)
1196
[5991]1197        t = self.domain.get_time()
[5081]1198
[5089]1199        if hasattr(self.function, 'time'):
[6410]1200            # Roll boundary over if time exceeds
[5089]1201            while t > self.function.time[-1]:
[6410]1202                msg = 'WARNING: domain time %.2f has exceeded' % t
[5089]1203                msg += 'time provided in '
1204                msg += 'transmissive_momentum_set_stage boundary object.\n'
1205                msg += 'I will continue, reusing the object from t==0'
1206                print msg
1207                t -= self.function.time[-1]
[5081]1208
1209        value = self.function(t)
[3804]1210        try:
1211            x = float(value)
[5081]1212        except:
[3804]1213            x = float(value[0])
[6410]1214
[3804]1215        q[0] = x
1216        return q
1217
[4769]1218        # FIXME: Consider this (taken from File_boundary) to allow
1219        # spatial variation
1220        # if vol_id is not None and edge_id is not None:
1221        #     i = self.boundary_indices[ vol_id, edge_id ]
1222        #     return self.F(t, point_id = i)
1223        # else:
1224        #     return self.F(t)
[3804]1225
1226
[6410]1227##
1228# @brief Deprecated boundary class.
[6045]1229class Transmissive_Momentum_Set_Stage_boundary(Transmissive_momentum_set_stage_boundary):
1230    pass
[3804]1231
[6045]1232
[6410]1233##
1234# @brief A transmissive boundary, momentum set to zero.
1235# @note Inherits from Bouondary.
[6045]1236class Transmissive_stage_zero_momentum_boundary(Boundary):
[6410]1237    """Return same stage as those present in its neighbour volume.
1238    Set momentum to zero.
[6045]1239
1240    Underlying domain must be specified when boundary is instantiated
[3804]1241    """
[6045]1242
[6410]1243    ##
1244    # @brief Instantiate a Transmissive (zero momentum) boundary.
1245    # @param domain
[6045]1246    def __init__(self, domain=None):
1247        Boundary.__init__(self)
1248
1249        if domain is None:
[6410]1250            msg = ('Domain must be specified for '
1251                   'Transmissive_stage_zero_momentum boundary')
[6045]1252            raise Exception, msg
1253
1254        self.domain = domain
1255
[6410]1256    ##
1257    # @brief Return a representation of this instance.
[6045]1258    def __repr__(self):
[6410]1259        return 'Transmissive_stage_zero_momentum_boundary(%s)' % self.domain
[6045]1260
[6410]1261    ##
1262    # @brief Calculate transmissive (zero momentum) results.
1263    # @param vol_id
1264    # @param edge_id
[6045]1265    def evaluate(self, vol_id, edge_id):
1266        """Transmissive boundaries return the edge values
1267        of the volume they serve.
1268        """
1269
1270        q = self.domain.get_conserved_quantities(vol_id, edge=edge_id)
[6410]1271
[6045]1272        q[1] = q[2] = 0.0
1273        return q
1274
1275
[6410]1276##
1277# @brief Class for a Dirichlet discharge boundary.
1278# @note Inherits from Boundary.
[6045]1279class Dirichlet_discharge_boundary(Boundary):
1280    """
[3804]1281    Sets stage (stage0)
1282    Sets momentum (wh0) in the inward normal direction.
1283
1284    Underlying domain must be specified when boundary is instantiated
1285    """
1286
[6410]1287    ##
1288    # @brief Instantiate a Dirichlet discharge boundary.
1289    # @param domain
1290    # @param stage0
1291    # @param wh0
[6045]1292    def __init__(self, domain=None, stage0=None, wh0=None):
[3804]1293        Boundary.__init__(self)
1294
1295        if domain is None:
[6410]1296            msg = 'Domain must be specified for this type of boundary'
1297            raise Exception, msg
[3804]1298
1299        if stage0 is None:
[6410]1300            raise Exception, 'Stage must be specified for this type of boundary'
[3804]1301
1302        if wh0 is None:
1303            wh0 = 0.0
1304
[6410]1305        self.domain = domain
1306        self.stage0 = stage0
[3804]1307        self.wh0 = wh0
1308
[6410]1309    ##
1310    # @brief Return a representation of this instance.
[3804]1311    def __repr__(self):
[6410]1312        return 'Dirichlet_Discharge_boundary(%s)' % self.domain
[3804]1313
[6410]1314    ##
1315    # @brief Calculate Dirichlet discharge boundary results.
1316    # @param vol_id
1317    # @param edge_id
[3804]1318    def evaluate(self, vol_id, edge_id):
[6410]1319        """Set discharge in the (inward) normal direction"""
[3804]1320
1321        normal = self.domain.get_normal(vol_id,edge_id)
1322        q = [self.stage0, -self.wh0*normal[0], -self.wh0*normal[1]]
1323        return q
1324
[4769]1325        # FIXME: Consider this (taken from File_boundary) to allow
1326        # spatial variation
1327        # if vol_id is not None and edge_id is not None:
1328        #     i = self.boundary_indices[ vol_id, edge_id ]
1329        #     return self.F(t, point_id = i)
1330        # else:
1331        #     return self.F(t)
[3804]1332
1333
[6410]1334# Backward compatibility
[6045]1335# FIXME(Ole): Deprecate
[6410]1336##
1337# @brief Deprecated
[6045]1338class Dirichlet_Discharge_boundary(Dirichlet_discharge_boundary):
1339    pass
1340
[6410]1341
1342##
1343# @brief A class for a 'field' boundary.
1344# @note Inherits from Boundary.
[4312]1345class Field_boundary(Boundary):
[6410]1346    """Set boundary from given field represented in an sww file containing
1347    values for stage, xmomentum and ymomentum.
[3804]1348
[6410]1349    Optionally, the user can specify mean_stage to offset the stage provided
1350    in the sww file.
1351
1352    This function is a thin wrapper around the generic File_boundary. The
[4754]1353    difference between the file_boundary and field_boundary is only that the
1354    field_boundary will allow you to change the level of the stage height when
[6410]1355    you read in the boundary condition. This is very useful when running
1356    different tide heights in the same area as you need only to convert one
1357    boundary condition to a SWW file, ideally for tide height of 0 m
[4754]1358    (saving disk space). Then you can use field_boundary to read this SWW file
1359    and change the stage height (tide) on the fly depending on the scenario.
[4312]1360    """
[3804]1361
[6410]1362    ##
1363    # @brief Construct an instance of a 'field' boundary.
1364    # @param filename Name of SWW file containing stage, x and ymomentum.
1365    # @param domain Shallow water domain for which the boundary applies.
1366    # @param mean_stage Mean water level added to stage derived from SWW file.
1367    # @param time_thinning Time step thinning factor.
1368    # @param time_limit
1369    # @param boundary_polygon
1370    # @param default_boundary None or an instance of Boundary.
1371    # @param use_cache True if caching is to be used.
1372    # @param verbose True if this method is to be verbose.
1373    def __init__(self,
1374                 filename,
1375                 domain,
[4312]1376                 mean_stage=0.0,
[5657]1377                 time_thinning=1,
[6173]1378                 time_limit=None,
[6410]1379                 boundary_polygon=None,
1380                 default_boundary=None,
[4312]1381                 use_cache=False,
[5657]1382                 verbose=False):
[4312]1383        """Constructor
1384
1385        filename: Name of sww file
1386        domain: pointer to shallow water domain for which the boundary applies
[4754]1387        mean_stage: The mean water level which will be added to stage derived
1388                    from the sww file
1389        time_thinning: Will set how many time steps from the sww file read in
[6410]1390                       will be interpolated to the boundary. For example if
[4754]1391                       the sww file has 1 second time steps and is 24 hours
[6410]1392                       in length it has 86400 time steps. If you set
1393                       time_thinning to 1 it will read all these steps.
[4754]1394                       If you set it to 100 it will read every 100th step eg
1395                       only 864 step. This parameter is very useful to increase
[6410]1396                       the speed of a model run that you are setting up
[4754]1397                       and testing.
[6410]1398
1399        default_boundary: Must be either None or an instance of a
[5657]1400                          class descending from class Boundary.
[6410]1401                          This will be used in case model time exceeds
[5657]1402                          that available in the underlying data.
[6410]1403
[4312]1404        use_cache:
1405        verbose:
1406        """
1407
1408        # Create generic file_boundary object
[6410]1409        self.file_boundary = File_boundary(filename,
1410                                           domain,
[4312]1411                                           time_thinning=time_thinning,
[6173]1412                                           time_limit=time_limit,
[5657]1413                                           boundary_polygon=boundary_polygon,
1414                                           default_boundary=default_boundary,
[4312]1415                                           use_cache=use_cache,
[5657]1416                                           verbose=verbose)
1417
[4312]1418        # Record information from File_boundary
1419        self.F = self.file_boundary.F
1420        self.domain = self.file_boundary.domain
[6410]1421
[4312]1422        # Record mean stage
1423        self.mean_stage = mean_stage
1424
[6410]1425    ##
1426    # @note Generate a string representation of this instance.
[4312]1427    def __repr__(self):
1428        return 'Field boundary'
1429
[6410]1430    ##
1431    # @brief Calculate 'field' boundary results.
1432    # @param vol_id
1433    # @param edge_id
[4312]1434    def evaluate(self, vol_id=None, edge_id=None):
1435        """Return linearly interpolated values based on domain.time
1436
1437        vol_id and edge_id are ignored
1438        """
[6410]1439
[4312]1440        # Evaluate file boundary
1441        q = self.file_boundary.evaluate(vol_id, edge_id)
1442
1443        # Adjust stage
1444        for j, name in enumerate(self.domain.conserved_quantities):
1445            if name == 'stage':
1446                q[j] += self.mean_stage
1447        return q
1448
1449
[6410]1450################################################################################
[4769]1451# Standard forcing terms
[6410]1452################################################################################
[4769]1453
[6410]1454##
1455# @brief Apply gravitational pull in the presence of bed slope.
1456# @param domain The domain to apply gravity to.
1457# @note Wrapper for C function gravity_c().
[3804]1458def gravity(domain):
1459    """Apply gravitational pull in the presence of bed slope
[4769]1460    Wrapper calls underlying C implementation
[3804]1461    """
1462
[6410]1463    from shallow_water_ext import gravity as gravity_c
1464
[3804]1465    xmom = domain.quantities['xmomentum'].explicit_update
1466    ymom = domain.quantities['ymomentum'].explicit_update
1467
[4687]1468    stage = domain.quantities['stage']
1469    elevation = domain.quantities['elevation']
[3804]1470
[4687]1471    h = stage.centroid_values - elevation.centroid_values
1472    z = elevation.vertex_values
1473
[3804]1474    x = domain.get_vertex_coordinates()
1475    g = domain.g
1476
[6410]1477    gravity_c(g, h, z, x, xmom, ymom)    #, 1.0e-6)
[3804]1478
[6410]1479##
1480# @brief Apply friction to a surface (implicit).
1481# @param domain The domain to apply Manning friction to.
1482# @note Wrapper for C function manning_friction_c().
[4769]1483def manning_friction_implicit(domain):
[6410]1484    """Apply (Manning) friction to water momentum
[4769]1485    Wrapper for c version
[3804]1486    """
1487
[6410]1488    from shallow_water_ext import manning_friction as manning_friction_c
[3804]1489
1490    xmom = domain.quantities['xmomentum']
1491    ymom = domain.quantities['ymomentum']
1492
1493    w = domain.quantities['stage'].centroid_values
1494    z = domain.quantities['elevation'].centroid_values
1495
1496    uh = xmom.centroid_values
1497    vh = ymom.centroid_values
1498    eta = domain.quantities['friction'].centroid_values
1499
1500    xmom_update = xmom.semi_implicit_update
1501    ymom_update = ymom.semi_implicit_update
1502
[3928]1503    N = len(domain)
[3804]1504    eps = domain.minimum_allowed_height
1505    g = domain.g
1506
[4769]1507    manning_friction_c(g, eps, w, z, uh, vh, eta, xmom_update, ymom_update)
[3804]1508
[6410]1509##
1510# @brief Apply friction to a surface (explicit).
1511# @param domain The domain to apply Manning friction to.
1512# @note Wrapper for C function manning_friction_c().
[4769]1513def manning_friction_explicit(domain):
[6410]1514    """Apply (Manning) friction to water momentum
[4769]1515    Wrapper for c version
[3804]1516    """
1517
[6410]1518    from shallow_water_ext import manning_friction as manning_friction_c
[3804]1519
1520    xmom = domain.quantities['xmomentum']
1521    ymom = domain.quantities['ymomentum']
1522
1523    w = domain.quantities['stage'].centroid_values
1524    z = domain.quantities['elevation'].centroid_values
1525
1526    uh = xmom.centroid_values
1527    vh = ymom.centroid_values
1528    eta = domain.quantities['friction'].centroid_values
1529
1530    xmom_update = xmom.explicit_update
1531    ymom_update = ymom.explicit_update
1532
[3928]1533    N = len(domain)
[3804]1534    eps = domain.minimum_allowed_height
1535    g = domain.g
1536
[4769]1537    manning_friction_c(g, eps, w, z, uh, vh, eta, xmom_update, ymom_update)
[3804]1538
1539
[4769]1540# FIXME (Ole): This was implemented for use with one of the analytical solutions (Sampson?)
[6410]1541##
1542# @brief Apply linear friction to a surface.
1543# @param domain The domain to apply Manning friction to.
1544# @note Is this still used (30 Oct 2007)?
[3804]1545def linear_friction(domain):
1546    """Apply linear friction to water momentum
1547
1548    Assumes quantity: 'linear_friction' to be present
1549    """
1550
1551    from math import sqrt
1552
1553    w = domain.quantities['stage'].centroid_values
1554    z = domain.quantities['elevation'].centroid_values
1555    h = w-z
1556
1557    uh = domain.quantities['xmomentum'].centroid_values
1558    vh = domain.quantities['ymomentum'].centroid_values
1559    tau = domain.quantities['linear_friction'].centroid_values
1560
1561    xmom_update = domain.quantities['xmomentum'].semi_implicit_update
1562    ymom_update = domain.quantities['ymomentum'].semi_implicit_update
1563
[3928]1564    N = len(domain) # number_of_triangles
[3804]1565    eps = domain.minimum_allowed_height
1566    g = domain.g #Not necessary? Why was this added?
1567
1568    for k in range(N):
1569        if tau[k] >= eps:
1570            if h[k] >= eps:
1571                S = -tau[k]/h[k]
1572
1573                #Update momentum
1574                xmom_update[k] += S*uh[k]
1575                ymom_update[k] += S*vh[k]
1576
[6410]1577################################################################################
1578# Experimental auxiliary functions
1579################################################################################
[3804]1580
[6410]1581##
1582# @brief Check forcefield parameter.
1583# @param f Object to check.
1584# @note 'f' may be a callable object or a scalar value.
[3804]1585def check_forcefield(f):
[6410]1586    """Check that force object is as expected.
1587   
1588    Check that f is either:
[3804]1589    1: a callable object f(t,x,y), where x and y are vectors
1590       and that it returns an array or a list of same length
1591       as x and y
1592    2: a scalar
1593    """
1594
1595    if callable(f):
1596        N = 3
[6304]1597        x = num.ones(3, num.float)
1598        y = num.ones(3, num.float)
[3804]1599        try:
1600            q = f(1.0, x=x, y=y)
1601        except Exception, e:
1602            msg = 'Function %s could not be executed:\n%s' %(f, e)
[4769]1603            # FIXME: Reconsider this semantics
[3804]1604            raise msg
1605
1606        try:
[6304]1607            q = num.array(q, num.float)
[3804]1608        except:
1609            msg = 'Return value from vector function %s could ' %f
[6304]1610            msg += 'not be converted into a numeric array of floats.\n'
[3804]1611            msg += 'Specified function should return either list or array.'
[6410]1612            raise Exception, msg
[3804]1613
[4769]1614        # Is this really what we want?
[6410]1615        msg = 'Function %s must return vector' % str(f)
1616        assert hasattr(q, 'len'), msg
1617
1618        msg = ('Return vector from function %s must have same '
1619               'length as input vectors' % f)
[3804]1620        assert len(q) == N, msg
1621    else:
1622        try:
1623            f = float(f)
1624        except:
[6410]1625            msg = ('Force field %s must be either a vector function or a '
1626                   'scalar value (coercible to float).' % str(f))
1627            raise Exception, msg
1628
[3804]1629    return f
1630
1631
[6410]1632##
1633# Class to apply a wind stress to a domain.
[3804]1634class Wind_stress:
1635    """Apply wind stress to water momentum in terms of
1636    wind speed [m/s] and wind direction [degrees]
1637    """
1638
[6410]1639    ##
1640    # @brief Create an instance of Wind_stress.
1641    # @param *args
1642    # @param **kwargs
[3804]1643    def __init__(self, *args, **kwargs):
1644        """Initialise windfield from wind speed s [m/s]
1645        and wind direction phi [degrees]
1646
1647        Inputs v and phi can be either scalars or Python functions, e.g.
1648
1649        W = Wind_stress(10, 178)
1650
1651        #FIXME - 'normal' degrees are assumed for now, i.e. the
1652        vector (1,0) has zero degrees.
1653        We may need to convert from 'compass' degrees later on and also
1654        map from True north to grid north.
1655
1656        Arguments can also be Python functions of t,x,y as in
1657
1658        def speed(t,x,y):
1659            ...
1660            return s
1661
1662        def angle(t,x,y):
1663            ...
1664            return phi
1665
1666        where x and y are vectors.
1667
1668        and then pass the functions in
1669
1670        W = Wind_stress(speed, angle)
1671
1672        The instantiated object W can be appended to the list of
1673        forcing_terms as in
1674
1675        Alternatively, one vector valued function for (speed, angle)
1676        can be applied, providing both quantities simultaneously.
1677        As in
1678        W = Wind_stress(F), where returns (speed, angle) for each t.
1679
1680        domain.forcing_terms.append(W)
1681        """
1682
1683        from anuga.config import rho_a, rho_w, eta_w
1684
1685        if len(args) == 2:
1686            s = args[0]
1687            phi = args[1]
1688        elif len(args) == 1:
[4769]1689            # Assume vector function returning (s, phi)(t,x,y)
[3804]1690            vector_function = args[0]
1691            s = lambda t,x,y: vector_function(t,x=x,y=y)[0]
1692            phi = lambda t,x,y: vector_function(t,x=x,y=y)[1]
1693        else:
[4769]1694           # Assume info is in 2 keyword arguments
[3804]1695           if len(kwargs) == 2:
1696               s = kwargs['s']
1697               phi = kwargs['phi']
1698           else:
[6410]1699               raise Exception, 'Assumes two keyword arguments: s=..., phi=....'
[3804]1700
1701        self.speed = check_forcefield(s)
1702        self.phi = check_forcefield(phi)
1703
1704        self.const = eta_w*rho_a/rho_w
1705
[6410]1706    ##
1707    # @brief 'execute' this class instance.
1708    # @param domain
[3804]1709    def __call__(self, domain):
[6410]1710        """Evaluate windfield based on values found in domain"""
[3804]1711
1712        from math import pi, cos, sin, sqrt
1713
1714        xmom_update = domain.quantities['xmomentum'].explicit_update
1715        ymom_update = domain.quantities['ymomentum'].explicit_update
1716
[6410]1717        N = len(domain)    # number_of_triangles
[3804]1718        t = domain.time
1719
1720        if callable(self.speed):
1721            xc = domain.get_centroid_coordinates()
1722            s_vec = self.speed(t, xc[:,0], xc[:,1])
1723        else:
[4769]1724            # Assume s is a scalar
[3804]1725            try:
[6304]1726                s_vec = self.speed * num.ones(N, num.float)
[3804]1727            except:
1728                msg = 'Speed must be either callable or a scalar: %s' %self.s
1729                raise msg
1730
1731        if callable(self.phi):
1732            xc = domain.get_centroid_coordinates()
1733            phi_vec = self.phi(t, xc[:,0], xc[:,1])
1734        else:
[4769]1735            # Assume phi is a scalar
[3804]1736
1737            try:
[6304]1738                phi_vec = self.phi * num.ones(N, num.float)
[3804]1739            except:
1740                msg = 'Angle must be either callable or a scalar: %s' %self.phi
1741                raise msg
1742
1743        assign_windfield_values(xmom_update, ymom_update,
1744                                s_vec, phi_vec, self.const)
1745
1746
[6410]1747##
1748# @brief Assign wind field values
1749# @param xmom_update
1750# @param ymom_update
1751# @param s_vec
1752# @param phi_vec
1753# @param const
[3804]1754def assign_windfield_values(xmom_update, ymom_update,
1755                            s_vec, phi_vec, const):
1756    """Python version of assigning wind field to update vectors.
1757    A c version also exists (for speed)
1758    """
[6410]1759
[3804]1760    from math import pi, cos, sin, sqrt
1761
1762    N = len(s_vec)
1763    for k in range(N):
1764        s = s_vec[k]
1765        phi = phi_vec[k]
1766
[4769]1767        # Convert to radians
[3804]1768        phi = phi*pi/180
1769
[4769]1770        # Compute velocity vector (u, v)
[3804]1771        u = s*cos(phi)
1772        v = s*sin(phi)
1773
[4769]1774        # Compute wind stress
[3804]1775        S = const * sqrt(u**2 + v**2)
1776        xmom_update[k] += S*u
1777        ymom_update[k] += S*v
1778
1779
[6410]1780##
1781# @brief A class for a general explicit forcing term.
[5294]1782class General_forcing:
[5736]1783    """General explicit forcing term for update of quantity
[6410]1784
[5294]1785    This is used by Inflow and Rainfall for instance
1786
[6410]1787
[5294]1788    General_forcing(quantity_name, rate, center, radius, polygon)
1789
1790    domain:     ANUGA computational domain
[6410]1791    quantity_name: Name of quantity to update.
[5736]1792                   It must be a known conserved quantity.
[6410]1793
1794    rate [?/s]: Total rate of change over the specified area.
[5294]1795                This parameter can be either a constant or a
[6410]1796                function of time. Positive values indicate increases,
[5294]1797                negative values indicate decreases.
1798                Rate can be None at initialisation but must be specified
[5436]1799                before forcing term is applied (i.e. simulation has started).
[5294]1800
1801    center [m]: Coordinates at center of flow point
1802    radius [m]: Size of circular area
[6002]1803    polygon:    Arbitrary polygon
[5873]1804    default_rate: Rate to be used if rate fails (e.g. if model time exceeds its data)
1805                  Admissible types: None, constant number or function of t
[5294]1806
1807
1808    Either center, radius or polygon can be specified but not both.
1809    If neither are specified the entire domain gets updated.
[6410]1810    All coordinates to be specified in absolute UTM coordinates (x, y) assuming the zone of domain.
1811
[6002]1812    Inflow or Rainfall for examples of use
[5294]1813    """
1814
1815
1816    # FIXME (AnyOne) : Add various methods to allow spatial variations
1817
[6410]1818    ##
1819    # @brief Create an instance of this forcing term.
1820    # @param domain
1821    # @param quantity_name
1822    # @param rate
1823    # @param center
1824    # @param radius
1825    # @param polygon
1826    # @param default_rate
1827    # @param verbose
[5294]1828    def __init__(self,
1829                 domain,
1830                 quantity_name,
[5303]1831                 rate=0.0,
[6410]1832                 center=None,
1833                 radius=None,
[5294]1834                 polygon=None,
[5873]1835                 default_rate=None,
[5294]1836                 verbose=False):
[6410]1837
1838        from math import pi, cos, sin
1839
[5450]1840        if center is None:
[6410]1841            msg = 'I got radius but no center.'
[5450]1842            assert radius is None, msg
[6410]1843
[5450]1844        if radius is None:
[6410]1845            msg += 'I got center but no radius.'
[5450]1846            assert center is None, msg
[5294]1847
1848        self.domain = domain
1849        self.quantity_name = quantity_name
1850        self.rate = rate
1851        self.center = ensure_numeric(center)
1852        self.radius = radius
[6410]1853        self.polygon = polygon
[5294]1854        self.verbose = verbose
[6410]1855        self.value = 0.0    # Can be used to remember value at
1856                            # previous timestep in order to obtain rate
[5294]1857
[6006]1858        # Get boundary (in absolute coordinates)
1859        bounding_polygon = domain.get_boundary_polygon()
[5570]1860
[5294]1861        # Update area if applicable
[6410]1862        self.exchange_area = None
[5294]1863        if center is not None and radius is not None:
1864            assert len(center) == 2
1865            msg = 'Polygon cannot be specified when center and radius are'
1866            assert polygon is None, msg
1867
[5436]1868            self.exchange_area = radius**2*pi
[5570]1869
1870            # Check that circle center lies within the mesh.
[6410]1871            msg = 'Center %s specified for forcing term did not' % str(center)
[5570]1872            msg += 'fall within the domain boundary.'
1873            assert is_inside_polygon(center, bounding_polygon), msg
1874
1875            # Check that circle periphery lies within the mesh.
1876            N = 100
1877            periphery_points = []
1878            for i in range(N):
[6410]1879                theta = 2*pi*i/100
[5570]1880
1881                x = center[0] + radius*cos(theta)
1882                y = center[1] + radius*sin(theta)
1883
1884                periphery_points.append([x,y])
1885
1886            for point in periphery_points:
[6410]1887                msg = 'Point %s on periphery for forcing term' % str(point)
[5736]1888                msg += ' did not fall within the domain boundary.'
[5570]1889                assert is_inside_polygon(point, bounding_polygon), msg
1890
[5294]1891        if polygon is not None:
[5570]1892            # Check that polygon lies within the mesh.
1893            for point in self.polygon:
[6410]1894                msg = 'Point %s in polygon for forcing term' % str(point)
[5736]1895                msg += ' did not fall within the domain boundary.'
[5570]1896                assert is_inside_polygon(point, bounding_polygon), msg
[6410]1897
1898            # Compute area and check that it is greater than 0
[6006]1899            self.exchange_area = polygon_area(self.polygon)
[5294]1900
[6410]1901            msg = 'Polygon %s in forcing term' % str(self.polygon)
1902            msg += ' has area = %f' % self.exchange_area
1903            assert self.exchange_area > 0.0
1904
[5294]1905        # Pointer to update vector
[5736]1906        self.update = domain.quantities[self.quantity_name].explicit_update
[5294]1907
1908        # Determine indices in flow area
[6410]1909        N = len(domain)
[5294]1910        points = domain.get_centroid_coordinates(absolute=True)
1911
[5436]1912        # Calculate indices in exchange area for this forcing term
1913        self.exchange_indices = None
[5294]1914        if self.center is not None and self.radius is not None:
1915            # Inlet is circular
[6410]1916            inlet_region = 'center=%s, radius=%s' % (self.center, self.radius)
1917
[5436]1918            self.exchange_indices = []
[5294]1919            for k in range(N):
[6410]1920                x, y = points[k,:]    # Centroid
1921
[5736]1922                c = self.center
1923                if ((x-c[0])**2+(y-c[1])**2) < self.radius**2:
[5436]1924                    self.exchange_indices.append(k)
[6410]1925
1926        if self.polygon is not None:
[5294]1927            # Inlet is polygon
[6410]1928            inlet_region = 'polygon=%s, area=%f m^2' % (self.polygon,
1929                                                        self.exchange_area)
1930
[5436]1931            self.exchange_indices = inside_polygon(points, self.polygon)
[6410]1932
[5450]1933        if self.exchange_indices is not None:
1934            if len(self.exchange_indices) == 0:
[5736]1935                msg = 'No triangles have been identified in '
[6410]1936                msg += 'specified region: %s' % inlet_region
[5450]1937                raise Exception, msg
[6410]1938
[5873]1939        # Check and store default_rate
[6410]1940        msg = ('Keyword argument default_rate must be either None '
1941               'or a function of time.\nI got %s.' % str(default_rate))
1942        assert (default_rate is None or
1943                type(default_rate) in [IntType, FloatType] or
1944                callable(default_rate)), msg
1945
[5873]1946        if default_rate is not None:
1947            # If it is a constant, make it a function
1948            if not callable(default_rate):
1949                tmp = default_rate
1950                default_rate = lambda t: tmp
[5294]1951
[5873]1952            # Check that default_rate is a function of one argument
1953            try:
1954                default_rate(0.0)
1955            except:
1956                raise Exception, msg
1957
1958        self.default_rate = default_rate
[6410]1959        self.default_rate_invoked = False    # Flag
[5873]1960
[6410]1961    ##
1962    # @brief Execute this instance.
1963    # @param domain
[5294]1964    def __call__(self, domain):
[6410]1965        """Apply inflow function at time specified in domain, update stage"""
[5294]1966
1967        # Call virtual method allowing local modifications
[5873]1968        t = domain.get_time()
1969        try:
1970            rate = self.update_rate(t)
[5874]1971        except Modeltime_too_early, e:
1972            raise Modeltime_too_early, e
1973        except Modeltime_too_late, e:
1974            if self.default_rate is None:
[6410]1975                raise Exception, e    # Reraise exception
[5873]1976            else:
[5874]1977                # Pass control to default rate function
[5873]1978                rate = self.default_rate(t)
[6410]1979
[5874]1980                if self.default_rate_invoked is False:
1981                    # Issue warning the first time
[6410]1982                    msg = ('%s\n'
1983                           'Instead I will use the default rate: %s\n'
1984                           'Note: Further warnings will be supressed'
1985                           % (str(e), str(self.default_rate)))
[5874]1986                    warn(msg)
[6410]1987
[5874]1988                    # FIXME (Ole): Replace this crude flag with
1989                    # Python's ability to print warnings only once.
1990                    # See http://docs.python.org/lib/warning-filter.html
1991                    self.default_rate_invoked = True
[5873]1992
[5294]1993        if rate is None:
[6410]1994            msg = ('Attribute rate must be specified in General_forcing '
1995                   'or its descendants before attempting to call it')
[5294]1996            raise Exception, msg
1997
1998        # Now rate is a number
1999        if self.verbose is True:
[6410]2000            print 'Rate of %s at time = %.2f = %f' % (self.quantity_name,
2001                                                      domain.get_time(),
2002                                                      rate)
[5294]2003
[5436]2004        if self.exchange_indices is None:
[5294]2005            self.update[:] += rate
2006        else:
2007            # Brute force assignment of restricted rate
[5436]2008            for k in self.exchange_indices:
[5294]2009                self.update[k] += rate
2010
[6410]2011    ##
2012    # @brief Update the internal rate.
2013    # @param t A callable or scalar used to set the rate.
2014    # @return The new rate.
[5294]2015    def update_rate(self, t):
2016        """Virtual method allowing local modifications by writing an
2017        overriding version in descendant
2018        """
[6410]2019
[6304]2020        if callable(self.rate):
2021            rate = self.rate(t)
2022        else:
2023            rate = self.rate
[5294]2024
2025        return rate
2026
[6410]2027    ##
2028    # @brief Get values for the specified quantity.
2029    # @param quantity_name Name of the quantity of interest.
2030    # @return The value(s) of the quantity.
2031    # @note If 'quantity_name' is None, use self.quantity_name.
2032    def get_quantity_values(self, quantity_name=None):
2033        """Return values for specified quantity restricted to opening
[5294]2034
[6121]2035        Optionally a quantity name can be specified if values from another
2036        quantity is sought
[5294]2037        """
[6410]2038
[6121]2039        if quantity_name is None:
2040            quantity_name = self.quantity_name
[6410]2041
[6121]2042        q = self.domain.quantities[quantity_name]
[6410]2043        return q.get_values(location='centroids',
[6121]2044                            indices=self.exchange_indices)
[5294]2045
[6410]2046    ##
2047    # @brief Set value for the specified quantity.
2048    # @param val The value object used to set value.
2049    # @param quantity_name Name of the quantity of interest.
2050    # @note If 'quantity_name' is None, use self.quantity_name.
[6121]2051    def set_quantity_values(self, val, quantity_name=None):
[6410]2052        """Set values for specified quantity restricted to opening
2053
[6121]2054        Optionally a quantity name can be specified if values from another
[6410]2055        quantity is sought
[5294]2056        """
2057
[6121]2058        if quantity_name is None:
2059            quantity_name = self.quantity_name
[5294]2060
[6410]2061        q = self.domain.quantities[self.quantity_name]
2062        q.set_values(val,
2063                     location='centroids',
2064                     indices=self.exchange_indices)
[5294]2065
[5736]2066
[6410]2067##
2068# @brief A class for rainfall forcing function.
2069# @note Inherits from General_forcing.
[5294]2070class Rainfall(General_forcing):
[4530]2071    """Class Rainfall - general 'rain over entire domain' forcing term.
[6410]2072
[4530]2073    Used for implementing Rainfall over the entire domain.
[6410]2074
[6304]2075        Current Limited to only One Gauge..
[6410]2076
2077        Need to add Spatial Varying Capability
[6304]2078        (This module came from copying and amending the Inflow Code)
[6410]2079
[4530]2080    Rainfall(rain)
[5294]2081
[6410]2082    domain
2083    rain [mm/s]:  Total rain rate over the specified domain.
[6304]2084                  NOTE: Raingauge Data needs to reflect the time step.
2085                  IE: if Gauge is mm read at a time step, then the input
[4530]2086                  here is as mm/(timeStep) so 10mm in 5minutes becomes
2087                  10/(5x60) = 0.0333mm/s.
[6304]2088
[4530]2089                  This parameter can be either a constant or a
[6410]2090                  function of time. Positive values indicate inflow,
[4530]2091                  negative values indicate outflow.
2092                  (and be used for Infiltration - Write Seperate Module)
2093                  The specified flow will be divided by the area of
2094                  the inflow region and then applied to update the
[5294]2095                  stage quantity.
[5178]2096
2097    polygon: Specifies a polygon to restrict the rainfall.
[6410]2098
[4530]2099    Examples
2100    How to put them in a run File...
[6304]2101
[5736]2102    #------------------------------------------------------------------------
[4530]2103    # Setup specialised forcing terms
[5736]2104    #------------------------------------------------------------------------
[4530]2105    # This is the new element implemented by Ole and Rudy to allow direct
[6055]2106    # input of Rainfall in mm/s
[4438]2107
[6410]2108    catchmentrainfall = Rainfall(rain=file_function('Q100_2hr_Rain.tms'))
[4530]2109                        # Note need path to File in String.
2110                        # Else assumed in same directory
2111
2112    domain.forcing_terms.append(catchmentrainfall)
2113    """
2114
[6410]2115    ##
2116    # @brief Create an instance of the class.
2117    # @param domain Domain of interest.
2118    # @param rate Total rain rate over the specified domain (mm/s).
2119    # @param center
2120    # @param radius
2121    # @param polygon Polygon  to restrict rainfall.
2122    # @param default_rate
2123    # @param verbose True if this instance is to be verbose.
[4530]2124    def __init__(self,
[5294]2125                 domain,
[6304]2126                 rate=0.0,
[6410]2127                 center=None,
2128                 radius=None,
[5294]2129                 polygon=None,
[6410]2130                 default_rate=None,
[5294]2131                 verbose=False):
[4530]2132
[5294]2133        # Converting mm/s to m/s to apply in ANUGA)
2134        if callable(rate):
2135            rain = lambda t: rate(t)/1000.0
2136        else:
[5873]2137            rain = rate/1000.0
2138
[6410]2139        if default_rate is not None:
[5873]2140            if callable(default_rate):
2141                default_rain = lambda t: default_rate(t)/1000.0
2142            else:
2143                default_rain = default_rate/1000.0
2144        else:
2145            default_rain = None
2146
[6410]2147        # pass to parent class
[5294]2148        General_forcing.__init__(self,
2149                                 domain,
2150                                 'stage',
2151                                 rate=rain,
[6410]2152                                 center=center,
2153                                 radius=radius,
[5294]2154                                 polygon=polygon,
[5874]2155                                 default_rate=default_rain,
[5294]2156                                 verbose=verbose)
[5178]2157
[4530]2158
[6410]2159##
2160# @brief A class for inflow (rain and drain) forcing function.
2161# @note Inherits from General_forcing.
[5294]2162class Inflow(General_forcing):
[4438]2163    """Class Inflow - general 'rain and drain' forcing term.
[6410]2164
[4438]2165    Useful for implementing flows in and out of the domain.
[6410]2166
[5294]2167    Inflow(flow, center, radius, polygon)
2168
2169    domain
[6410]2170    rate [m^3/s]: Total flow rate over the specified area.
[4438]2171                  This parameter can be either a constant or a
[6410]2172                  function of time. Positive values indicate inflow,
[4438]2173                  negative values indicate outflow.
2174                  The specified flow will be divided by the area of
[6410]2175                  the inflow region and then applied to update stage.
[5294]2176    center [m]: Coordinates at center of flow point
2177    radius [m]: Size of circular area
2178    polygon:    Arbitrary polygon.
2179
2180    Either center, radius or polygon must be specified
[6410]2181
[4438]2182    Examples
2183
2184    # Constant drain at 0.003 m^3/s.
2185    # The outflow area is 0.07**2*pi=0.0154 m^2
[6410]2186    # This corresponds to a rate of change of 0.003/0.0154 = 0.2 m/s
2187    #
[4438]2188    Inflow((0.7, 0.4), 0.07, -0.003)
2189
2190
2191    # Tap turning up to a maximum inflow of 0.0142 m^3/s.
2192    # The inflow area is 0.03**2*pi = 0.00283 m^2
[6410]2193    # This corresponds to a rate of change of 0.0142/0.00283 = 5 m/s
[4438]2194    # over the specified area
2195    Inflow((0.5, 0.5), 0.03, lambda t: min(0.01*t, 0.0142))
[4530]2196
[5294]2197
[5736]2198    #------------------------------------------------------------------------
[4530]2199    # Setup specialised forcing terms
[5736]2200    #------------------------------------------------------------------------
[4530]2201    # This is the new element implemented by Ole to allow direct input
2202    # of Inflow in m^3/s
2203
2204    hydrograph = Inflow(center=(320, 300), radius=10,
[5506]2205                        rate=file_function('Q/QPMF_Rot_Sub13.tms'))
[4530]2206
2207    domain.forcing_terms.append(hydrograph)
[4438]2208    """
2209
[6410]2210    ##
2211    # @brief Create an instance of the class.
2212    # @param domain Domain of interest.
2213    # @param rate Total rain rate over the specified domain (mm/s).
2214    # @param center
2215    # @param radius
2216    # @param polygon Polygon  to restrict rainfall.
2217    # @param default_rate
2218    # @param verbose True if this instance is to be verbose.
[4438]2219    def __init__(self,
[5294]2220                 domain,
[6304]2221                 rate=0.0,
[6410]2222                 center=None,
2223                 radius=None,
[5294]2224                 polygon=None,
[5874]2225                 default_rate=None,
[6410]2226                 verbose=False):
[5294]2227        # Create object first to make area is available
2228        General_forcing.__init__(self,
2229                                 domain,
2230                                 'stage',
2231                                 rate=rate,
[6410]2232                                 center=center,
2233                                 radius=radius,
[5294]2234                                 polygon=polygon,
[5873]2235                                 default_rate=default_rate,
[5294]2236                                 verbose=verbose)
2237
[6410]2238    ##
2239    # @brief Update the instance rate.
2240    # @param t New rate object.
[5294]2241    def update_rate(self, t):
2242        """Virtual method allowing local modifications by writing an
2243        overriding version in descendant
2244
[6410]2245        This one converts m^3/s to m/s which can be added directly
[5736]2246        to 'stage' in ANUGA
[5294]2247        """
2248
[6304]2249        if callable(self.rate):
2250            _rate = self.rate(t)/self.exchange_area
2251        else:
2252            _rate = self.rate/self.exchange_area
[4438]2253
[5294]2254        return _rate
[4438]2255
2256
[6410]2257################################################################################
[4733]2258# Initialise module
[6410]2259################################################################################
[3804]2260
2261from anuga.utilities import compile
2262if compile.can_use_C_extension('shallow_water_ext.c'):
[6410]2263    # Underlying C implementations can be accessed
[3804]2264    from shallow_water_ext import rotate, assign_windfield_values
[4769]2265else:
[6410]2266    msg = 'C implementations could not be accessed by %s.\n ' % __file__
[4769]2267    msg += 'Make sure compile_all.py has been run as described in '
2268    msg += 'the ANUGA installation guide.'
2269    raise Exception, msg
[3804]2270
[4733]2271# Optimisation with psyco
[3804]2272from anuga.config import use_psyco
2273if use_psyco:
2274    try:
2275        import psyco
2276    except:
2277        import os
[5920]2278        if os.name == 'posix' and os.uname()[4] in ['x86_64', 'ia64']:
[3804]2279            pass
2280            #Psyco isn't supported on 64 bit systems, but it doesn't matter
2281        else:
[6410]2282            msg = ('WARNING: psyco (speedup) could not be imported, '
2283                   'you may want to consider installing it')
[3804]2284            print msg
2285    else:
2286        psyco.bind(Domain.distribute_to_vertices_and_edges)
2287        psyco.bind(Domain.compute_fluxes)
2288
[6410]2289
[3804]2290if __name__ == "__main__":
2291    pass
Note: See TracBrowser for help on using the repository browser.