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

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

Back-merge from Numeric trunk to numpy branch.

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