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

Last change on this file since 6441 was 6441, checked in by rwilson, 15 years ago

After changes to get_absolute, ensure_numeric, etc.

  • Property svn:keywords set to LastChangedDate LastChangedRevision LastChangedBy HeadURL Id
File size: 78.4 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: 6441 $)
71ModifiedBy:
[4005]72    $Author: rwilson $
73    $Date: 2009-03-03 21:26:22 +0000 (Tue, 03 Mar 2009) $
[3804]74"""
75
[4769]76# Subversion keywords:
[3804]77#
[4769]78# $LastChangedDate: 2009-03-03 21:26:22 +0000 (Tue, 03 Mar 2009) $
79# $LastChangedRevision: 6441 $
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?
[6441]1615        # info is "(func name, filename, defining line)"
1616        func_info = (f.func_name, f.func_code.co_filename,
1617                     f.func_code.co_firstlineno)
1618        msg = ('Function %s() must return vector (defined in %s, line %d)'
1619               % func_info)
[6410]1620        assert hasattr(q, 'len'), msg
1621
[6441]1622        msg = ('Return vector from function %s() must have same '
1623               'length as input vectors\nq=%s' % (f.func_name, str(q)))
[3804]1624        assert len(q) == N, msg
1625    else:
1626        try:
1627            f = float(f)
1628        except:
[6441]1629            msg = ('Force field %s must be a scalar value coercible to float.'
1630                   % str(f))
[6410]1631            raise Exception, msg
1632
[3804]1633    return f
1634
1635
[6410]1636##
1637# Class to apply a wind stress to a domain.
[3804]1638class Wind_stress:
1639    """Apply wind stress to water momentum in terms of
1640    wind speed [m/s] and wind direction [degrees]
1641    """
1642
[6410]1643    ##
1644    # @brief Create an instance of Wind_stress.
1645    # @param *args
1646    # @param **kwargs
[3804]1647    def __init__(self, *args, **kwargs):
1648        """Initialise windfield from wind speed s [m/s]
1649        and wind direction phi [degrees]
1650
1651        Inputs v and phi can be either scalars or Python functions, e.g.
1652
1653        W = Wind_stress(10, 178)
1654
1655        #FIXME - 'normal' degrees are assumed for now, i.e. the
1656        vector (1,0) has zero degrees.
1657        We may need to convert from 'compass' degrees later on and also
1658        map from True north to grid north.
1659
1660        Arguments can also be Python functions of t,x,y as in
1661
1662        def speed(t,x,y):
1663            ...
1664            return s
1665
1666        def angle(t,x,y):
1667            ...
1668            return phi
1669
1670        where x and y are vectors.
1671
1672        and then pass the functions in
1673
1674        W = Wind_stress(speed, angle)
1675
1676        The instantiated object W can be appended to the list of
1677        forcing_terms as in
1678
1679        Alternatively, one vector valued function for (speed, angle)
1680        can be applied, providing both quantities simultaneously.
1681        As in
1682        W = Wind_stress(F), where returns (speed, angle) for each t.
1683
1684        domain.forcing_terms.append(W)
1685        """
1686
1687        from anuga.config import rho_a, rho_w, eta_w
1688
1689        if len(args) == 2:
1690            s = args[0]
1691            phi = args[1]
1692        elif len(args) == 1:
[4769]1693            # Assume vector function returning (s, phi)(t,x,y)
[3804]1694            vector_function = args[0]
1695            s = lambda t,x,y: vector_function(t,x=x,y=y)[0]
1696            phi = lambda t,x,y: vector_function(t,x=x,y=y)[1]
1697        else:
[4769]1698           # Assume info is in 2 keyword arguments
[3804]1699           if len(kwargs) == 2:
1700               s = kwargs['s']
1701               phi = kwargs['phi']
1702           else:
[6410]1703               raise Exception, 'Assumes two keyword arguments: s=..., phi=....'
[3804]1704
1705        self.speed = check_forcefield(s)
1706        self.phi = check_forcefield(phi)
1707
1708        self.const = eta_w*rho_a/rho_w
1709
[6410]1710    ##
1711    # @brief 'execute' this class instance.
1712    # @param domain
[3804]1713    def __call__(self, domain):
[6410]1714        """Evaluate windfield based on values found in domain"""
[3804]1715
1716        from math import pi, cos, sin, sqrt
1717
1718        xmom_update = domain.quantities['xmomentum'].explicit_update
1719        ymom_update = domain.quantities['ymomentum'].explicit_update
1720
[6410]1721        N = len(domain)    # number_of_triangles
[3804]1722        t = domain.time
1723
1724        if callable(self.speed):
1725            xc = domain.get_centroid_coordinates()
1726            s_vec = self.speed(t, xc[:,0], xc[:,1])
1727        else:
[4769]1728            # Assume s is a scalar
[3804]1729            try:
[6304]1730                s_vec = self.speed * num.ones(N, num.float)
[3804]1731            except:
1732                msg = 'Speed must be either callable or a scalar: %s' %self.s
1733                raise msg
1734
1735        if callable(self.phi):
1736            xc = domain.get_centroid_coordinates()
1737            phi_vec = self.phi(t, xc[:,0], xc[:,1])
1738        else:
[4769]1739            # Assume phi is a scalar
[3804]1740
1741            try:
[6304]1742                phi_vec = self.phi * num.ones(N, num.float)
[3804]1743            except:
1744                msg = 'Angle must be either callable or a scalar: %s' %self.phi
1745                raise msg
1746
1747        assign_windfield_values(xmom_update, ymom_update,
1748                                s_vec, phi_vec, self.const)
1749
1750
[6410]1751##
1752# @brief Assign wind field values
1753# @param xmom_update
1754# @param ymom_update
1755# @param s_vec
1756# @param phi_vec
1757# @param const
[3804]1758def assign_windfield_values(xmom_update, ymom_update,
1759                            s_vec, phi_vec, const):
1760    """Python version of assigning wind field to update vectors.
1761    A c version also exists (for speed)
1762    """
[6410]1763
[3804]1764    from math import pi, cos, sin, sqrt
1765
1766    N = len(s_vec)
1767    for k in range(N):
1768        s = s_vec[k]
1769        phi = phi_vec[k]
1770
[4769]1771        # Convert to radians
[3804]1772        phi = phi*pi/180
1773
[4769]1774        # Compute velocity vector (u, v)
[3804]1775        u = s*cos(phi)
1776        v = s*sin(phi)
1777
[4769]1778        # Compute wind stress
[3804]1779        S = const * sqrt(u**2 + v**2)
1780        xmom_update[k] += S*u
1781        ymom_update[k] += S*v
1782
1783
[6410]1784##
1785# @brief A class for a general explicit forcing term.
[5294]1786class General_forcing:
[5736]1787    """General explicit forcing term for update of quantity
[6410]1788
[5294]1789    This is used by Inflow and Rainfall for instance
1790
[6410]1791
[5294]1792    General_forcing(quantity_name, rate, center, radius, polygon)
1793
1794    domain:     ANUGA computational domain
[6410]1795    quantity_name: Name of quantity to update.
[5736]1796                   It must be a known conserved quantity.
[6410]1797
1798    rate [?/s]: Total rate of change over the specified area.
[5294]1799                This parameter can be either a constant or a
[6410]1800                function of time. Positive values indicate increases,
[5294]1801                negative values indicate decreases.
1802                Rate can be None at initialisation but must be specified
[5436]1803                before forcing term is applied (i.e. simulation has started).
[5294]1804
1805    center [m]: Coordinates at center of flow point
1806    radius [m]: Size of circular area
[6002]1807    polygon:    Arbitrary polygon
[5873]1808    default_rate: Rate to be used if rate fails (e.g. if model time exceeds its data)
1809                  Admissible types: None, constant number or function of t
[5294]1810
1811
1812    Either center, radius or polygon can be specified but not both.
1813    If neither are specified the entire domain gets updated.
[6410]1814    All coordinates to be specified in absolute UTM coordinates (x, y) assuming the zone of domain.
1815
[6002]1816    Inflow or Rainfall for examples of use
[5294]1817    """
1818
1819
1820    # FIXME (AnyOne) : Add various methods to allow spatial variations
1821
[6410]1822    ##
1823    # @brief Create an instance of this forcing term.
1824    # @param domain
1825    # @param quantity_name
1826    # @param rate
1827    # @param center
1828    # @param radius
1829    # @param polygon
1830    # @param default_rate
1831    # @param verbose
[5294]1832    def __init__(self,
1833                 domain,
1834                 quantity_name,
[5303]1835                 rate=0.0,
[6410]1836                 center=None,
1837                 radius=None,
[5294]1838                 polygon=None,
[5873]1839                 default_rate=None,
[5294]1840                 verbose=False):
[6410]1841
1842        from math import pi, cos, sin
1843
[5450]1844        if center is None:
[6410]1845            msg = 'I got radius but no center.'
[5450]1846            assert radius is None, msg
[6410]1847
[5450]1848        if radius is None:
[6410]1849            msg += 'I got center but no radius.'
[5450]1850            assert center is None, msg
[5294]1851
1852        self.domain = domain
1853        self.quantity_name = quantity_name
1854        self.rate = rate
1855        self.center = ensure_numeric(center)
1856        self.radius = radius
[6410]1857        self.polygon = polygon
[5294]1858        self.verbose = verbose
[6410]1859        self.value = 0.0    # Can be used to remember value at
1860                            # previous timestep in order to obtain rate
[5294]1861
[6006]1862        # Get boundary (in absolute coordinates)
1863        bounding_polygon = domain.get_boundary_polygon()
[5570]1864
[5294]1865        # Update area if applicable
[6410]1866        self.exchange_area = None
[5294]1867        if center is not None and radius is not None:
1868            assert len(center) == 2
1869            msg = 'Polygon cannot be specified when center and radius are'
1870            assert polygon is None, msg
1871
[5436]1872            self.exchange_area = radius**2*pi
[5570]1873
1874            # Check that circle center lies within the mesh.
[6410]1875            msg = 'Center %s specified for forcing term did not' % str(center)
[5570]1876            msg += 'fall within the domain boundary.'
1877            assert is_inside_polygon(center, bounding_polygon), msg
1878
1879            # Check that circle periphery lies within the mesh.
1880            N = 100
1881            periphery_points = []
1882            for i in range(N):
[6410]1883                theta = 2*pi*i/100
[5570]1884
1885                x = center[0] + radius*cos(theta)
1886                y = center[1] + radius*sin(theta)
1887
1888                periphery_points.append([x,y])
1889
1890            for point in periphery_points:
[6410]1891                msg = 'Point %s on periphery for forcing term' % str(point)
[5736]1892                msg += ' did not fall within the domain boundary.'
[5570]1893                assert is_inside_polygon(point, bounding_polygon), msg
1894
[5294]1895        if polygon is not None:
[5570]1896            # Check that polygon lies within the mesh.
1897            for point in self.polygon:
[6410]1898                msg = 'Point %s in polygon for forcing term' % str(point)
[5736]1899                msg += ' did not fall within the domain boundary.'
[5570]1900                assert is_inside_polygon(point, bounding_polygon), msg
[6410]1901
1902            # Compute area and check that it is greater than 0
[6006]1903            self.exchange_area = polygon_area(self.polygon)
[5294]1904
[6410]1905            msg = 'Polygon %s in forcing term' % str(self.polygon)
1906            msg += ' has area = %f' % self.exchange_area
1907            assert self.exchange_area > 0.0
1908
[5294]1909        # Pointer to update vector
[5736]1910        self.update = domain.quantities[self.quantity_name].explicit_update
[5294]1911
1912        # Determine indices in flow area
[6410]1913        N = len(domain)
[5294]1914        points = domain.get_centroid_coordinates(absolute=True)
1915
[5436]1916        # Calculate indices in exchange area for this forcing term
1917        self.exchange_indices = None
[5294]1918        if self.center is not None and self.radius is not None:
1919            # Inlet is circular
[6410]1920            inlet_region = 'center=%s, radius=%s' % (self.center, self.radius)
1921
[5436]1922            self.exchange_indices = []
[5294]1923            for k in range(N):
[6410]1924                x, y = points[k,:]    # Centroid
1925
[5736]1926                c = self.center
1927                if ((x-c[0])**2+(y-c[1])**2) < self.radius**2:
[5436]1928                    self.exchange_indices.append(k)
[6410]1929
1930        if self.polygon is not None:
[5294]1931            # Inlet is polygon
[6441]1932            inlet_region = ('polygon=\n%s, area=%f m^2' %
1933                            (self.polygon, self.exchange_area))
[6410]1934
[5436]1935            self.exchange_indices = inside_polygon(points, self.polygon)
[6410]1936
[5450]1937        if self.exchange_indices is not None:
1938            if len(self.exchange_indices) == 0:
[5736]1939                msg = 'No triangles have been identified in '
[6410]1940                msg += 'specified region: %s' % inlet_region
[5450]1941                raise Exception, msg
[6410]1942
[5873]1943        # Check and store default_rate
[6410]1944        msg = ('Keyword argument default_rate must be either None '
1945               'or a function of time.\nI got %s.' % str(default_rate))
1946        assert (default_rate is None or
1947                type(default_rate) in [IntType, FloatType] or
1948                callable(default_rate)), msg
1949
[5873]1950        if default_rate is not None:
1951            # If it is a constant, make it a function
1952            if not callable(default_rate):
1953                tmp = default_rate
1954                default_rate = lambda t: tmp
[5294]1955
[5873]1956            # Check that default_rate is a function of one argument
1957            try:
1958                default_rate(0.0)
1959            except:
1960                raise Exception, msg
1961
1962        self.default_rate = default_rate
[6410]1963        self.default_rate_invoked = False    # Flag
[5873]1964
[6410]1965    ##
1966    # @brief Execute this instance.
1967    # @param domain
[5294]1968    def __call__(self, domain):
[6410]1969        """Apply inflow function at time specified in domain, update stage"""
[5294]1970
1971        # Call virtual method allowing local modifications
[5873]1972        t = domain.get_time()
1973        try:
1974            rate = self.update_rate(t)
[5874]1975        except Modeltime_too_early, e:
1976            raise Modeltime_too_early, e
1977        except Modeltime_too_late, e:
1978            if self.default_rate is None:
[6410]1979                raise Exception, e    # Reraise exception
[5873]1980            else:
[5874]1981                # Pass control to default rate function
[5873]1982                rate = self.default_rate(t)
[6410]1983
[5874]1984                if self.default_rate_invoked is False:
1985                    # Issue warning the first time
[6410]1986                    msg = ('%s\n'
1987                           'Instead I will use the default rate: %s\n'
1988                           'Note: Further warnings will be supressed'
1989                           % (str(e), str(self.default_rate)))
[5874]1990                    warn(msg)
[6410]1991
[5874]1992                    # FIXME (Ole): Replace this crude flag with
1993                    # Python's ability to print warnings only once.
1994                    # See http://docs.python.org/lib/warning-filter.html
1995                    self.default_rate_invoked = True
[5873]1996
[5294]1997        if rate is None:
[6410]1998            msg = ('Attribute rate must be specified in General_forcing '
1999                   'or its descendants before attempting to call it')
[5294]2000            raise Exception, msg
2001
2002        # Now rate is a number
2003        if self.verbose is True:
[6410]2004            print 'Rate of %s at time = %.2f = %f' % (self.quantity_name,
2005                                                      domain.get_time(),
2006                                                      rate)
[5294]2007
[5436]2008        if self.exchange_indices is None:
[5294]2009            self.update[:] += rate
2010        else:
2011            # Brute force assignment of restricted rate
[5436]2012            for k in self.exchange_indices:
[5294]2013                self.update[k] += rate
2014
[6410]2015    ##
2016    # @brief Update the internal rate.
2017    # @param t A callable or scalar used to set the rate.
2018    # @return The new rate.
[5294]2019    def update_rate(self, t):
2020        """Virtual method allowing local modifications by writing an
2021        overriding version in descendant
2022        """
[6410]2023
[6304]2024        if callable(self.rate):
2025            rate = self.rate(t)
2026        else:
2027            rate = self.rate
[5294]2028
2029        return rate
2030
[6410]2031    ##
2032    # @brief Get values for the specified quantity.
2033    # @param quantity_name Name of the quantity of interest.
2034    # @return The value(s) of the quantity.
2035    # @note If 'quantity_name' is None, use self.quantity_name.
2036    def get_quantity_values(self, quantity_name=None):
2037        """Return values for specified quantity restricted to opening
[5294]2038
[6121]2039        Optionally a quantity name can be specified if values from another
2040        quantity is sought
[5294]2041        """
[6410]2042
[6121]2043        if quantity_name is None:
2044            quantity_name = self.quantity_name
[6410]2045
[6121]2046        q = self.domain.quantities[quantity_name]
[6410]2047        return q.get_values(location='centroids',
[6121]2048                            indices=self.exchange_indices)
[5294]2049
[6410]2050    ##
2051    # @brief Set value for the specified quantity.
2052    # @param val The value object used to set value.
2053    # @param quantity_name Name of the quantity of interest.
2054    # @note If 'quantity_name' is None, use self.quantity_name.
[6121]2055    def set_quantity_values(self, val, quantity_name=None):
[6410]2056        """Set values for specified quantity restricted to opening
2057
[6121]2058        Optionally a quantity name can be specified if values from another
[6410]2059        quantity is sought
[5294]2060        """
2061
[6121]2062        if quantity_name is None:
2063            quantity_name = self.quantity_name
[5294]2064
[6410]2065        q = self.domain.quantities[self.quantity_name]
2066        q.set_values(val,
2067                     location='centroids',
2068                     indices=self.exchange_indices)
[5294]2069
[5736]2070
[6410]2071##
2072# @brief A class for rainfall forcing function.
2073# @note Inherits from General_forcing.
[5294]2074class Rainfall(General_forcing):
[4530]2075    """Class Rainfall - general 'rain over entire domain' forcing term.
[6410]2076
[4530]2077    Used for implementing Rainfall over the entire domain.
[6410]2078
[6304]2079        Current Limited to only One Gauge..
[6410]2080
2081        Need to add Spatial Varying Capability
[6304]2082        (This module came from copying and amending the Inflow Code)
[6410]2083
[4530]2084    Rainfall(rain)
[5294]2085
[6410]2086    domain
2087    rain [mm/s]:  Total rain rate over the specified domain.
[6304]2088                  NOTE: Raingauge Data needs to reflect the time step.
2089                  IE: if Gauge is mm read at a time step, then the input
[4530]2090                  here is as mm/(timeStep) so 10mm in 5minutes becomes
2091                  10/(5x60) = 0.0333mm/s.
[6304]2092
[4530]2093                  This parameter can be either a constant or a
[6410]2094                  function of time. Positive values indicate inflow,
[4530]2095                  negative values indicate outflow.
2096                  (and be used for Infiltration - Write Seperate Module)
2097                  The specified flow will be divided by the area of
2098                  the inflow region and then applied to update the
[5294]2099                  stage quantity.
[5178]2100
2101    polygon: Specifies a polygon to restrict the rainfall.
[6410]2102
[4530]2103    Examples
2104    How to put them in a run File...
[6304]2105
[5736]2106    #------------------------------------------------------------------------
[4530]2107    # Setup specialised forcing terms
[5736]2108    #------------------------------------------------------------------------
[4530]2109    # This is the new element implemented by Ole and Rudy to allow direct
[6055]2110    # input of Rainfall in mm/s
[4438]2111
[6410]2112    catchmentrainfall = Rainfall(rain=file_function('Q100_2hr_Rain.tms'))
[4530]2113                        # Note need path to File in String.
2114                        # Else assumed in same directory
2115
2116    domain.forcing_terms.append(catchmentrainfall)
2117    """
2118
[6410]2119    ##
2120    # @brief Create an instance of the class.
2121    # @param domain Domain of interest.
2122    # @param rate Total rain rate over the specified domain (mm/s).
2123    # @param center
2124    # @param radius
2125    # @param polygon Polygon  to restrict rainfall.
2126    # @param default_rate
2127    # @param verbose True if this instance is to be verbose.
[4530]2128    def __init__(self,
[5294]2129                 domain,
[6304]2130                 rate=0.0,
[6410]2131                 center=None,
2132                 radius=None,
[5294]2133                 polygon=None,
[6410]2134                 default_rate=None,
[5294]2135                 verbose=False):
[4530]2136
[5294]2137        # Converting mm/s to m/s to apply in ANUGA)
2138        if callable(rate):
2139            rain = lambda t: rate(t)/1000.0
2140        else:
[5873]2141            rain = rate/1000.0
2142
[6410]2143        if default_rate is not None:
[5873]2144            if callable(default_rate):
2145                default_rain = lambda t: default_rate(t)/1000.0
2146            else:
2147                default_rain = default_rate/1000.0
2148        else:
2149            default_rain = None
2150
[6410]2151        # pass to parent class
[5294]2152        General_forcing.__init__(self,
2153                                 domain,
2154                                 'stage',
2155                                 rate=rain,
[6410]2156                                 center=center,
2157                                 radius=radius,
[5294]2158                                 polygon=polygon,
[5874]2159                                 default_rate=default_rain,
[5294]2160                                 verbose=verbose)
[5178]2161
[4530]2162
[6410]2163##
2164# @brief A class for inflow (rain and drain) forcing function.
2165# @note Inherits from General_forcing.
[5294]2166class Inflow(General_forcing):
[4438]2167    """Class Inflow - general 'rain and drain' forcing term.
[6410]2168
[4438]2169    Useful for implementing flows in and out of the domain.
[6410]2170
[5294]2171    Inflow(flow, center, radius, polygon)
2172
2173    domain
[6410]2174    rate [m^3/s]: Total flow rate over the specified area.
[4438]2175                  This parameter can be either a constant or a
[6410]2176                  function of time. Positive values indicate inflow,
[4438]2177                  negative values indicate outflow.
2178                  The specified flow will be divided by the area of
[6410]2179                  the inflow region and then applied to update stage.
[5294]2180    center [m]: Coordinates at center of flow point
2181    radius [m]: Size of circular area
2182    polygon:    Arbitrary polygon.
2183
2184    Either center, radius or polygon must be specified
[6410]2185
[4438]2186    Examples
2187
2188    # Constant drain at 0.003 m^3/s.
2189    # The outflow area is 0.07**2*pi=0.0154 m^2
[6410]2190    # This corresponds to a rate of change of 0.003/0.0154 = 0.2 m/s
2191    #
[4438]2192    Inflow((0.7, 0.4), 0.07, -0.003)
2193
2194
2195    # Tap turning up to a maximum inflow of 0.0142 m^3/s.
2196    # The inflow area is 0.03**2*pi = 0.00283 m^2
[6410]2197    # This corresponds to a rate of change of 0.0142/0.00283 = 5 m/s
[4438]2198    # over the specified area
2199    Inflow((0.5, 0.5), 0.03, lambda t: min(0.01*t, 0.0142))
[4530]2200
[5294]2201
[5736]2202    #------------------------------------------------------------------------
[4530]2203    # Setup specialised forcing terms
[5736]2204    #------------------------------------------------------------------------
[4530]2205    # This is the new element implemented by Ole to allow direct input
2206    # of Inflow in m^3/s
2207
2208    hydrograph = Inflow(center=(320, 300), radius=10,
[5506]2209                        rate=file_function('Q/QPMF_Rot_Sub13.tms'))
[4530]2210
2211    domain.forcing_terms.append(hydrograph)
[4438]2212    """
2213
[6410]2214    ##
2215    # @brief Create an instance of the class.
2216    # @param domain Domain of interest.
2217    # @param rate Total rain rate over the specified domain (mm/s).
2218    # @param center
2219    # @param radius
2220    # @param polygon Polygon  to restrict rainfall.
2221    # @param default_rate
2222    # @param verbose True if this instance is to be verbose.
[4438]2223    def __init__(self,
[5294]2224                 domain,
[6304]2225                 rate=0.0,
[6410]2226                 center=None,
2227                 radius=None,
[5294]2228                 polygon=None,
[5874]2229                 default_rate=None,
[6410]2230                 verbose=False):
[5294]2231        # Create object first to make area is available
2232        General_forcing.__init__(self,
2233                                 domain,
2234                                 'stage',
2235                                 rate=rate,
[6410]2236                                 center=center,
2237                                 radius=radius,
[5294]2238                                 polygon=polygon,
[5873]2239                                 default_rate=default_rate,
[5294]2240                                 verbose=verbose)
2241
[6410]2242    ##
2243    # @brief Update the instance rate.
2244    # @param t New rate object.
[5294]2245    def update_rate(self, t):
2246        """Virtual method allowing local modifications by writing an
2247        overriding version in descendant
2248
[6410]2249        This one converts m^3/s to m/s which can be added directly
[5736]2250        to 'stage' in ANUGA
[5294]2251        """
2252
[6304]2253        if callable(self.rate):
2254            _rate = self.rate(t)/self.exchange_area
2255        else:
2256            _rate = self.rate/self.exchange_area
[4438]2257
[5294]2258        return _rate
[4438]2259
2260
[6410]2261################################################################################
[4733]2262# Initialise module
[6410]2263################################################################################
[3804]2264
2265from anuga.utilities import compile
2266if compile.can_use_C_extension('shallow_water_ext.c'):
[6410]2267    # Underlying C implementations can be accessed
[3804]2268    from shallow_water_ext import rotate, assign_windfield_values
[4769]2269else:
[6410]2270    msg = 'C implementations could not be accessed by %s.\n ' % __file__
[4769]2271    msg += 'Make sure compile_all.py has been run as described in '
2272    msg += 'the ANUGA installation guide.'
2273    raise Exception, msg
[3804]2274
[4733]2275# Optimisation with psyco
[3804]2276from anuga.config import use_psyco
2277if use_psyco:
2278    try:
2279        import psyco
2280    except:
2281        import os
[5920]2282        if os.name == 'posix' and os.uname()[4] in ['x86_64', 'ia64']:
[3804]2283            pass
2284            #Psyco isn't supported on 64 bit systems, but it doesn't matter
2285        else:
[6410]2286            msg = ('WARNING: psyco (speedup) could not be imported, '
2287                   'you may want to consider installing it')
[3804]2288            print msg
2289    else:
2290        psyco.bind(Domain.distribute_to_vertices_and_edges)
2291        psyco.bind(Domain.compute_fluxes)
2292
[6410]2293
[3804]2294if __name__ == "__main__":
2295    pass
Note: See TracBrowser for help on using the repository browser.