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

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

After changes to get_absolute, ensure_numeric, etc.

  • Property svn:keywords set to LastChangedDate LastChangedRevision LastChangedBy HeadURL Id
File size: 78.4 KB
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: 6441 $)
71ModifiedBy:
72    $Author: rwilson $
73    $Date: 2009-03-03 21:26:22 +0000 (Tue, 03 Mar 2009) $
74"""
75
76# Subversion keywords:
77#
78# $LastChangedDate: 2009-03-03 21:26:22 +0000 (Tue, 03 Mar 2009) $
79# $LastChangedRevision: 6441 $
80# $LastChangedBy: rwilson $
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)    # super integrity check
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# End of class Shallow Water Domain
817################################################################################
818
819#-----------------
820# Flux computation
821#-----------------
822
823## @brief Compute fluxes and timestep suitable for all volumes in domain.
824# @param domain The domain to calculate fluxes for.
825def compute_fluxes(domain):
826    """Compute fluxes and timestep suitable for all volumes in domain.
827
828    Compute total flux for each conserved quantity using "flux_function"
829
830    Fluxes across each edge are scaled by edgelengths and summed up
831    Resulting flux is then scaled by area and stored in
832    explicit_update for each of the three conserved quantities
833    stage, xmomentum and ymomentum
834
835    The maximal allowable speed computed by the flux_function for each volume
836    is converted to a timestep that must not be exceeded. The minimum of
837    those is computed as the next overall timestep.
838
839    Post conditions:
840      domain.explicit_update is reset to computed flux values
841      domain.timestep is set to the largest step satisfying all volumes.
842
843    This wrapper calls the underlying C version of compute fluxes
844    """
845
846    import sys
847    from shallow_water_ext import compute_fluxes_ext_central \
848                                  as compute_fluxes_ext
849
850    N = len(domain)    # number_of_triangles
851
852    # Shortcuts
853    Stage = domain.quantities['stage']
854    Xmom = domain.quantities['xmomentum']
855    Ymom = domain.quantities['ymomentum']
856    Bed = domain.quantities['elevation']
857
858    timestep = float(sys.maxint)
859
860    flux_timestep = compute_fluxes_ext(timestep,
861                                       domain.epsilon,
862                                       domain.H0,
863                                       domain.g,
864                                       domain.neighbours,
865                                       domain.neighbour_edges,
866                                       domain.normals,
867                                       domain.edgelengths,
868                                       domain.radii,
869                                       domain.areas,
870                                       domain.tri_full_flag,
871                                       Stage.edge_values,
872                                       Xmom.edge_values,
873                                       Ymom.edge_values,
874                                       Bed.edge_values,
875                                       Stage.boundary_values,
876                                       Xmom.boundary_values,
877                                       Ymom.boundary_values,
878                                       Stage.explicit_update,
879                                       Xmom.explicit_update,
880                                       Ymom.explicit_update,
881                                       domain.already_computed_flux,
882                                       domain.max_speed,
883                                       int(domain.optimise_dry_cells))
884
885    domain.flux_timestep = flux_timestep
886
887################################################################################
888# Module functions for gradient limiting
889################################################################################
890
891##
892# @brief Wrapper for C version of extrapolate_second_order_sw.
893# @param domain The domain to operate on.
894# @note MH090605 The following method belongs to the shallow_water domain class
895#       see comments in the corresponding method in shallow_water_ext.c
896def extrapolate_second_order_sw(domain):
897    """Wrapper calling C version of extrapolate_second_order_sw"""
898
899    import sys
900    from shallow_water_ext import extrapolate_second_order_sw as extrapol2
901
902    N = len(domain) # number_of_triangles
903
904    # Shortcuts
905    Stage = domain.quantities['stage']
906    Xmom = domain.quantities['xmomentum']
907    Ymom = domain.quantities['ymomentum']
908    Elevation = domain.quantities['elevation']
909
910    extrapol2(domain,
911              domain.surrogate_neighbours,
912              domain.number_of_boundaries,
913              domain.centroid_coordinates,
914              Stage.centroid_values,
915              Xmom.centroid_values,
916              Ymom.centroid_values,
917              Elevation.centroid_values,
918              domain.vertex_coordinates,
919              Stage.vertex_values,
920              Xmom.vertex_values,
921              Ymom.vertex_values,
922              Elevation.vertex_values,
923              int(domain.optimise_dry_cells))
924
925##
926# @brief Distribution from centroids to vertices specific to the SWW eqn.
927# @param domain The domain to operate on.
928def distribute_using_vertex_limiter(domain):
929    """Distribution from centroids to vertices specific to the SWW equation.
930
931    It will ensure that h (w-z) is always non-negative even in the
932    presence of steep bed-slopes by taking a weighted average between shallow
933    and deep cases.
934
935    In addition, all conserved quantities get distributed as per either a
936    constant (order==1) or a piecewise linear function (order==2).
937
938    FIXME: more explanation about removal of artificial variability etc
939
940    Precondition:
941      All quantities defined at centroids and bed elevation defined at
942      vertices.
943
944    Postcondition
945      Conserved quantities defined at vertices
946    """
947
948    # Remove very thin layers of water
949    protect_against_infinitesimal_and_negative_heights(domain)
950
951    # Extrapolate all conserved quantities
952    if domain.optimised_gradient_limiter:
953        # MH090605 if second order,
954        # perform the extrapolation and limiting on
955        # all of the conserved quantities
956
957        if (domain._order_ == 1):
958            for name in domain.conserved_quantities:
959                Q = domain.quantities[name]
960                Q.extrapolate_first_order()
961        elif domain._order_ == 2:
962            domain.extrapolate_second_order_sw()
963        else:
964            raise 'Unknown order'
965    else:
966        # Old code:
967        for name in domain.conserved_quantities:
968            Q = domain.quantities[name]
969
970            if domain._order_ == 1:
971                Q.extrapolate_first_order()
972            elif domain._order_ == 2:
973                Q.extrapolate_second_order_and_limit_by_vertex()
974            else:
975                raise 'Unknown order'
976
977    # Take bed elevation into account when water heights are small
978    balance_deep_and_shallow(domain)
979
980    # Compute edge values by interpolation
981    for name in domain.conserved_quantities:
982        Q = domain.quantities[name]
983        Q.interpolate_from_vertices_to_edges()
984
985##
986# @brief Distribution from centroids to edges specific to the SWW eqn.
987# @param domain The domain to operate on.
988def distribute_using_edge_limiter(domain):
989    """Distribution from centroids to edges specific to the SWW eqn.
990
991    It will ensure that h (w-z) is always non-negative even in the
992    presence of steep bed-slopes by taking a weighted average between shallow
993    and deep cases.
994
995    In addition, all conserved quantities get distributed as per either a
996    constant (order==1) or a piecewise linear function (order==2).
997
998
999    Precondition:
1000      All quantities defined at centroids and bed elevation defined at
1001      vertices.
1002
1003    Postcondition
1004      Conserved quantities defined at vertices
1005    """
1006
1007    # Remove very thin layers of water
1008    protect_against_infinitesimal_and_negative_heights(domain)
1009
1010    for name in domain.conserved_quantities:
1011        Q = domain.quantities[name]
1012        if domain._order_ == 1:
1013            Q.extrapolate_first_order()
1014        elif domain._order_ == 2:
1015            Q.extrapolate_second_order_and_limit_by_edge()
1016        else:
1017            raise 'Unknown order'
1018
1019    balance_deep_and_shallow(domain)
1020
1021    # Compute edge values by interpolation
1022    for name in domain.conserved_quantities:
1023        Q = domain.quantities[name]
1024        Q.interpolate_from_vertices_to_edges()
1025
1026##
1027# @brief  Protect against infinitesimal heights and associated high velocities.
1028# @param domain The domain to operate on.
1029def protect_against_infinitesimal_and_negative_heights(domain):
1030    """Protect against infinitesimal heights and associated high velocities"""
1031
1032    from shallow_water_ext import protect
1033
1034    # Shortcuts
1035    wc = domain.quantities['stage'].centroid_values
1036    zc = domain.quantities['elevation'].centroid_values
1037    xmomc = domain.quantities['xmomentum'].centroid_values
1038    ymomc = domain.quantities['ymomentum'].centroid_values
1039
1040    protect(domain.minimum_allowed_height, domain.maximum_allowed_speed,
1041            domain.epsilon, wc, zc, xmomc, ymomc)
1042
1043##
1044# @brief Wrapper for C function balance_deep_and_shallow_c().
1045# @param domain The domain to operate on.
1046def balance_deep_and_shallow(domain):
1047    """Compute linear combination between stage as computed by
1048    gradient-limiters limiting using w, and stage computed by
1049    gradient-limiters limiting using h (h-limiter).
1050    The former takes precedence when heights are large compared to the
1051    bed slope while the latter takes precedence when heights are
1052    relatively small.  Anything in between is computed as a balanced
1053    linear combination in order to avoid numerical disturbances which
1054    would otherwise appear as a result of hard switching between
1055    modes.
1056
1057    Wrapper for C implementation
1058    """
1059
1060    from shallow_water_ext import balance_deep_and_shallow \
1061                                  as balance_deep_and_shallow_c
1062
1063    # Shortcuts
1064    wc = domain.quantities['stage'].centroid_values
1065    zc = domain.quantities['elevation'].centroid_values
1066    wv = domain.quantities['stage'].vertex_values
1067    zv = domain.quantities['elevation'].vertex_values
1068
1069    # Momentums at centroids
1070    xmomc = domain.quantities['xmomentum'].centroid_values
1071    ymomc = domain.quantities['ymomentum'].centroid_values
1072
1073    # Momentums at vertices
1074    xmomv = domain.quantities['xmomentum'].vertex_values
1075    ymomv = domain.quantities['ymomentum'].vertex_values
1076
1077    balance_deep_and_shallow_c(domain,
1078                               wc, zc, wv, zv, wc,
1079                               xmomc, ymomc, xmomv, ymomv)
1080
1081
1082################################################################################
1083# Boundary conditions - specific to the shallow water wave equation
1084################################################################################
1085
1086##
1087# @brief Class for a reflective boundary.
1088# @note Inherits from Boundary.
1089class Reflective_boundary(Boundary):
1090    """Reflective boundary returns same conserved quantities as
1091    those present in its neighbour volume but reflected.
1092
1093    This class is specific to the shallow water equation as it
1094    works with the momentum quantities assumed to be the second
1095    and third conserved quantities.
1096    """
1097
1098    ##
1099    # @brief Instantiate a Reflective_boundary.
1100    # @param domain
1101    def __init__(self, domain=None):
1102        Boundary.__init__(self)
1103
1104        if domain is None:
1105            msg = 'Domain must be specified for reflective boundary'
1106            raise Exception, msg
1107
1108        # Handy shorthands
1109        self.stage = domain.quantities['stage'].edge_values
1110        self.xmom = domain.quantities['xmomentum'].edge_values
1111        self.ymom = domain.quantities['ymomentum'].edge_values
1112        self.normals = domain.normals
1113
1114        self.conserved_quantities = num.zeros(3, num.float)
1115
1116    ##
1117    # @brief Return a representation of this instance.
1118    def __repr__(self):
1119        return 'Reflective_boundary'
1120
1121    ##
1122    # @brief Calculate reflections (reverse outward momentum).
1123    # @param vol_id
1124    # @param edge_id
1125    def evaluate(self, vol_id, edge_id):
1126        """Reflective boundaries reverses the outward momentum
1127        of the volume they serve.
1128        """
1129
1130        q = self.conserved_quantities
1131        q[0] = self.stage[vol_id, edge_id]
1132        q[1] = self.xmom[vol_id, edge_id]
1133        q[2] = self.ymom[vol_id, edge_id]
1134
1135        normal = self.normals[vol_id, 2*edge_id:2*edge_id+2]
1136
1137        r = rotate(q, normal, direction = 1)
1138        r[1] = -r[1]
1139        q = rotate(r, normal, direction = -1)
1140
1141        return q
1142
1143
1144##
1145# @brief Class for a transmissive boundary.
1146# @note Inherits from Boundary.
1147class Transmissive_momentum_set_stage_boundary(Boundary):
1148    """Returns same momentum conserved quantities as
1149    those present in its neighbour volume.
1150    Sets stage by specifying a function f of time which may either be a
1151    vector function or a scalar function
1152
1153    Example:
1154
1155    def waveform(t):
1156        return sea_level + normalized_amplitude/cosh(t-25)**2
1157
1158    Bts = Transmissive_momentum_set_stage_boundary(domain, waveform)
1159
1160    Underlying domain must be specified when boundary is instantiated
1161    """
1162
1163    ##
1164    # @brief Instantiate a Reflective_boundary.
1165    # @param domain
1166    # @param function
1167    def __init__(self, domain=None, function=None):
1168        Boundary.__init__(self)
1169
1170        if domain is None:
1171            msg = 'Domain must be specified for this type boundary'
1172            raise Exception, msg
1173
1174        if function is None:
1175            msg = 'Function must be specified for this type boundary'
1176            raise Exception, msg
1177
1178        self.domain = domain
1179        self.function = function
1180
1181    ##
1182    # @brief Return a representation of this instance.
1183    def __repr__(self):
1184        return 'Transmissive_momentum_set_stage_boundary(%s)' %self.domain
1185
1186    ##
1187    # @brief Calculate transmissive results.
1188    # @param vol_id
1189    # @param edge_id
1190    def evaluate(self, vol_id, edge_id):
1191        """Transmissive momentum set stage boundaries return the edge momentum
1192        values of the volume they serve.
1193        """
1194
1195        q = self.domain.get_conserved_quantities(vol_id, edge = edge_id)
1196
1197        t = self.domain.get_time()
1198
1199        if hasattr(self.function, 'time'):
1200            # Roll boundary over if time exceeds
1201            while t > self.function.time[-1]:
1202                msg = 'WARNING: domain time %.2f has exceeded' % t
1203                msg += 'time provided in '
1204                msg += 'transmissive_momentum_set_stage boundary object.\n'
1205                msg += 'I will continue, reusing the object from t==0'
1206                print msg
1207                t -= self.function.time[-1]
1208
1209        value = self.function(t)
1210        try:
1211            x = float(value)
1212        except:
1213            x = float(value[0])
1214
1215        q[0] = x
1216        return q
1217
1218        # FIXME: Consider this (taken from File_boundary) to allow
1219        # spatial variation
1220        # if vol_id is not None and edge_id is not None:
1221        #     i = self.boundary_indices[ vol_id, edge_id ]
1222        #     return self.F(t, point_id = i)
1223        # else:
1224        #     return self.F(t)
1225
1226
1227##
1228# @brief Deprecated boundary class.
1229class Transmissive_Momentum_Set_Stage_boundary(Transmissive_momentum_set_stage_boundary):
1230    pass
1231
1232
1233##
1234# @brief A transmissive boundary, momentum set to zero.
1235# @note Inherits from Bouondary.
1236class Transmissive_stage_zero_momentum_boundary(Boundary):
1237    """Return same stage as those present in its neighbour volume.
1238    Set momentum to zero.
1239
1240    Underlying domain must be specified when boundary is instantiated
1241    """
1242
1243    ##
1244    # @brief Instantiate a Transmissive (zero momentum) boundary.
1245    # @param domain
1246    def __init__(self, domain=None):
1247        Boundary.__init__(self)
1248
1249        if domain is None:
1250            msg = ('Domain must be specified for '
1251                   'Transmissive_stage_zero_momentum boundary')
1252            raise Exception, msg
1253
1254        self.domain = domain
1255
1256    ##
1257    # @brief Return a representation of this instance.
1258    def __repr__(self):
1259        return 'Transmissive_stage_zero_momentum_boundary(%s)' % self.domain
1260
1261    ##
1262    # @brief Calculate transmissive (zero momentum) results.
1263    # @param vol_id
1264    # @param edge_id
1265    def evaluate(self, vol_id, edge_id):
1266        """Transmissive boundaries return the edge values
1267        of the volume they serve.
1268        """
1269
1270        q = self.domain.get_conserved_quantities(vol_id, edge=edge_id)
1271
1272        q[1] = q[2] = 0.0
1273        return q
1274
1275
1276##
1277# @brief Class for a Dirichlet discharge boundary.
1278# @note Inherits from Boundary.
1279class Dirichlet_discharge_boundary(Boundary):
1280    """
1281    Sets stage (stage0)
1282    Sets momentum (wh0) in the inward normal direction.
1283
1284    Underlying domain must be specified when boundary is instantiated
1285    """
1286
1287    ##
1288    # @brief Instantiate a Dirichlet discharge boundary.
1289    # @param domain
1290    # @param stage0
1291    # @param wh0
1292    def __init__(self, domain=None, stage0=None, wh0=None):
1293        Boundary.__init__(self)
1294
1295        if domain is None:
1296            msg = 'Domain must be specified for this type of boundary'
1297            raise Exception, msg
1298
1299        if stage0 is None:
1300            raise Exception, 'Stage must be specified for this type of boundary'
1301
1302        if wh0 is None:
1303            wh0 = 0.0
1304
1305        self.domain = domain
1306        self.stage0 = stage0
1307        self.wh0 = wh0
1308
1309    ##
1310    # @brief Return a representation of this instance.
1311    def __repr__(self):
1312        return 'Dirichlet_Discharge_boundary(%s)' % self.domain
1313
1314    ##
1315    # @brief Calculate Dirichlet discharge boundary results.
1316    # @param vol_id
1317    # @param edge_id
1318    def evaluate(self, vol_id, edge_id):
1319        """Set discharge in the (inward) normal direction"""
1320
1321        normal = self.domain.get_normal(vol_id,edge_id)
1322        q = [self.stage0, -self.wh0*normal[0], -self.wh0*normal[1]]
1323        return q
1324
1325        # FIXME: Consider this (taken from File_boundary) to allow
1326        # spatial variation
1327        # if vol_id is not None and edge_id is not None:
1328        #     i = self.boundary_indices[ vol_id, edge_id ]
1329        #     return self.F(t, point_id = i)
1330        # else:
1331        #     return self.F(t)
1332
1333
1334# Backward compatibility
1335# FIXME(Ole): Deprecate
1336##
1337# @brief Deprecated
1338class Dirichlet_Discharge_boundary(Dirichlet_discharge_boundary):
1339    pass
1340
1341
1342##
1343# @brief A class for a 'field' boundary.
1344# @note Inherits from Boundary.
1345class Field_boundary(Boundary):
1346    """Set boundary from given field represented in an sww file containing
1347    values for stage, xmomentum and ymomentum.
1348
1349    Optionally, the user can specify mean_stage to offset the stage provided
1350    in the sww file.
1351
1352    This function is a thin wrapper around the generic File_boundary. The
1353    difference between the file_boundary and field_boundary is only that the
1354    field_boundary will allow you to change the level of the stage height when
1355    you read in the boundary condition. This is very useful when running
1356    different tide heights in the same area as you need only to convert one
1357    boundary condition to a SWW file, ideally for tide height of 0 m
1358    (saving disk space). Then you can use field_boundary to read this SWW file
1359    and change the stage height (tide) on the fly depending on the scenario.
1360    """
1361
1362    ##
1363    # @brief Construct an instance of a 'field' boundary.
1364    # @param filename Name of SWW file containing stage, x and ymomentum.
1365    # @param domain Shallow water domain for which the boundary applies.
1366    # @param mean_stage Mean water level added to stage derived from SWW file.
1367    # @param time_thinning Time step thinning factor.
1368    # @param time_limit
1369    # @param boundary_polygon
1370    # @param default_boundary None or an instance of Boundary.
1371    # @param use_cache True if caching is to be used.
1372    # @param verbose True if this method is to be verbose.
1373    def __init__(self,
1374                 filename,
1375                 domain,
1376                 mean_stage=0.0,
1377                 time_thinning=1,
1378                 time_limit=None,
1379                 boundary_polygon=None,
1380                 default_boundary=None,
1381                 use_cache=False,
1382                 verbose=False):
1383        """Constructor
1384
1385        filename: Name of sww file
1386        domain: pointer to shallow water domain for which the boundary applies
1387        mean_stage: The mean water level which will be added to stage derived
1388                    from the sww file
1389        time_thinning: Will set how many time steps from the sww file read in
1390                       will be interpolated to the boundary. For example if
1391                       the sww file has 1 second time steps and is 24 hours
1392                       in length it has 86400 time steps. If you set
1393                       time_thinning to 1 it will read all these steps.
1394                       If you set it to 100 it will read every 100th step eg
1395                       only 864 step. This parameter is very useful to increase
1396                       the speed of a model run that you are setting up
1397                       and testing.
1398
1399        default_boundary: Must be either None or an instance of a
1400                          class descending from class Boundary.
1401                          This will be used in case model time exceeds
1402                          that available in the underlying data.
1403
1404        use_cache:
1405        verbose:
1406        """
1407
1408        # Create generic file_boundary object
1409        self.file_boundary = File_boundary(filename,
1410                                           domain,
1411                                           time_thinning=time_thinning,
1412                                           time_limit=time_limit,
1413                                           boundary_polygon=boundary_polygon,
1414                                           default_boundary=default_boundary,
1415                                           use_cache=use_cache,
1416                                           verbose=verbose)
1417
1418        # Record information from File_boundary
1419        self.F = self.file_boundary.F
1420        self.domain = self.file_boundary.domain
1421
1422        # Record mean stage
1423        self.mean_stage = mean_stage
1424
1425    ##
1426    # @note Generate a string representation of this instance.
1427    def __repr__(self):
1428        return 'Field boundary'
1429
1430    ##
1431    # @brief Calculate 'field' boundary results.
1432    # @param vol_id
1433    # @param edge_id
1434    def evaluate(self, vol_id=None, edge_id=None):
1435        """Return linearly interpolated values based on domain.time
1436
1437        vol_id and edge_id are ignored
1438        """
1439
1440        # Evaluate file boundary
1441        q = self.file_boundary.evaluate(vol_id, edge_id)
1442
1443        # Adjust stage
1444        for j, name in enumerate(self.domain.conserved_quantities):
1445            if name == 'stage':
1446                q[j] += self.mean_stage
1447        return q
1448
1449
1450################################################################################
1451# Standard forcing terms
1452################################################################################
1453
1454##
1455# @brief Apply gravitational pull in the presence of bed slope.
1456# @param domain The domain to apply gravity to.
1457# @note Wrapper for C function gravity_c().
1458def gravity(domain):
1459    """Apply gravitational pull in the presence of bed slope
1460    Wrapper calls underlying C implementation
1461    """
1462
1463    from shallow_water_ext import gravity as gravity_c
1464
1465    xmom = domain.quantities['xmomentum'].explicit_update
1466    ymom = domain.quantities['ymomentum'].explicit_update
1467
1468    stage = domain.quantities['stage']
1469    elevation = domain.quantities['elevation']
1470
1471    h = stage.centroid_values - elevation.centroid_values
1472    z = elevation.vertex_values
1473
1474    x = domain.get_vertex_coordinates()
1475    g = domain.g
1476
1477    gravity_c(g, h, z, x, xmom, ymom)    #, 1.0e-6)
1478
1479##
1480# @brief Apply friction to a surface (implicit).
1481# @param domain The domain to apply Manning friction to.
1482# @note Wrapper for C function manning_friction_c().
1483def manning_friction_implicit(domain):
1484    """Apply (Manning) friction to water momentum
1485    Wrapper for c version
1486    """
1487
1488    from shallow_water_ext import manning_friction as manning_friction_c
1489
1490    xmom = domain.quantities['xmomentum']
1491    ymom = domain.quantities['ymomentum']
1492
1493    w = domain.quantities['stage'].centroid_values
1494    z = domain.quantities['elevation'].centroid_values
1495
1496    uh = xmom.centroid_values
1497    vh = ymom.centroid_values
1498    eta = domain.quantities['friction'].centroid_values
1499
1500    xmom_update = xmom.semi_implicit_update
1501    ymom_update = ymom.semi_implicit_update
1502
1503    N = len(domain)
1504    eps = domain.minimum_allowed_height
1505    g = domain.g
1506
1507    manning_friction_c(g, eps, w, z, uh, vh, eta, xmom_update, ymom_update)
1508
1509##
1510# @brief Apply friction to a surface (explicit).
1511# @param domain The domain to apply Manning friction to.
1512# @note Wrapper for C function manning_friction_c().
1513def manning_friction_explicit(domain):
1514    """Apply (Manning) friction to water momentum
1515    Wrapper for c version
1516    """
1517
1518    from shallow_water_ext import manning_friction as manning_friction_c
1519
1520    xmom = domain.quantities['xmomentum']
1521    ymom = domain.quantities['ymomentum']
1522
1523    w = domain.quantities['stage'].centroid_values
1524    z = domain.quantities['elevation'].centroid_values
1525
1526    uh = xmom.centroid_values
1527    vh = ymom.centroid_values
1528    eta = domain.quantities['friction'].centroid_values
1529
1530    xmom_update = xmom.explicit_update
1531    ymom_update = ymom.explicit_update
1532
1533    N = len(domain)
1534    eps = domain.minimum_allowed_height
1535    g = domain.g
1536
1537    manning_friction_c(g, eps, w, z, uh, vh, eta, xmom_update, ymom_update)
1538
1539
1540# FIXME (Ole): This was implemented for use with one of the analytical solutions (Sampson?)
1541##
1542# @brief Apply linear friction to a surface.
1543# @param domain The domain to apply Manning friction to.
1544# @note Is this still used (30 Oct 2007)?
1545def linear_friction(domain):
1546    """Apply linear friction to water momentum
1547
1548    Assumes quantity: 'linear_friction' to be present
1549    """
1550
1551    from math import sqrt
1552
1553    w = domain.quantities['stage'].centroid_values
1554    z = domain.quantities['elevation'].centroid_values
1555    h = w-z
1556
1557    uh = domain.quantities['xmomentum'].centroid_values
1558    vh = domain.quantities['ymomentum'].centroid_values
1559    tau = domain.quantities['linear_friction'].centroid_values
1560
1561    xmom_update = domain.quantities['xmomentum'].semi_implicit_update
1562    ymom_update = domain.quantities['ymomentum'].semi_implicit_update
1563
1564    N = len(domain) # number_of_triangles
1565    eps = domain.minimum_allowed_height
1566    g = domain.g #Not necessary? Why was this added?
1567
1568    for k in range(N):
1569        if tau[k] >= eps:
1570            if h[k] >= eps:
1571                S = -tau[k]/h[k]
1572
1573                #Update momentum
1574                xmom_update[k] += S*uh[k]
1575                ymom_update[k] += S*vh[k]
1576
1577################################################################################
1578# Experimental auxiliary functions
1579################################################################################
1580
1581##
1582# @brief Check forcefield parameter.
1583# @param f Object to check.
1584# @note 'f' may be a callable object or a scalar value.
1585def check_forcefield(f):
1586    """Check that force object is as expected.
1587   
1588    Check that f is either:
1589    1: a callable object f(t,x,y), where x and y are vectors
1590       and that it returns an array or a list of same length
1591       as x and y
1592    2: a scalar
1593    """
1594
1595    if callable(f):
1596        N = 3
1597        x = num.ones(3, num.float)
1598        y = num.ones(3, num.float)
1599        try:
1600            q = f(1.0, x=x, y=y)
1601        except Exception, e:
1602            msg = 'Function %s could not be executed:\n%s' %(f, e)
1603            # FIXME: Reconsider this semantics
1604            raise msg
1605
1606        try:
1607            q = num.array(q, num.float)
1608        except:
1609            msg = 'Return value from vector function %s could ' %f
1610            msg += 'not be converted into a numeric array of floats.\n'
1611            msg += 'Specified function should return either list or array.'
1612            raise Exception, msg
1613
1614        # Is this really what we want?
1615        # info is "(func name, filename, defining line)"
1616        func_info = (f.func_name, f.func_code.co_filename,
1617                     f.func_code.co_firstlineno)
1618        msg = ('Function %s() must return vector (defined in %s, line %d)'
1619               % func_info)
1620        assert hasattr(q, 'len'), msg
1621
1622        msg = ('Return vector from function %s() must have same '
1623               'length as input vectors\nq=%s' % (f.func_name, str(q)))
1624        assert len(q) == N, msg
1625    else:
1626        try:
1627            f = float(f)
1628        except:
1629            msg = ('Force field %s must be a scalar value coercible to float.'
1630                   % str(f))
1631            raise Exception, msg
1632
1633    return f
1634
1635
1636##
1637# Class to apply a wind stress to a domain.
1638class Wind_stress:
1639    """Apply wind stress to water momentum in terms of
1640    wind speed [m/s] and wind direction [degrees]
1641    """
1642
1643    ##
1644    # @brief Create an instance of Wind_stress.
1645    # @param *args
1646    # @param **kwargs
1647    def __init__(self, *args, **kwargs):
1648        """Initialise windfield from wind speed s [m/s]
1649        and wind direction phi [degrees]
1650
1651        Inputs v and phi can be either scalars or Python functions, e.g.
1652
1653        W = Wind_stress(10, 178)
1654
1655        #FIXME - 'normal' degrees are assumed for now, i.e. the
1656        vector (1,0) has zero degrees.
1657        We may need to convert from 'compass' degrees later on and also
1658        map from True north to grid north.
1659
1660        Arguments can also be Python functions of t,x,y as in
1661
1662        def speed(t,x,y):
1663            ...
1664            return s
1665
1666        def angle(t,x,y):
1667            ...
1668            return phi
1669
1670        where x and y are vectors.
1671
1672        and then pass the functions in
1673
1674        W = Wind_stress(speed, angle)
1675
1676        The instantiated object W can be appended to the list of
1677        forcing_terms as in
1678
1679        Alternatively, one vector valued function for (speed, angle)
1680        can be applied, providing both quantities simultaneously.
1681        As in
1682        W = Wind_stress(F), where returns (speed, angle) for each t.
1683
1684        domain.forcing_terms.append(W)
1685        """
1686
1687        from anuga.config import rho_a, rho_w, eta_w
1688
1689        if len(args) == 2:
1690            s = args[0]
1691            phi = args[1]
1692        elif len(args) == 1:
1693            # Assume vector function returning (s, phi)(t,x,y)
1694            vector_function = args[0]
1695            s = lambda t,x,y: vector_function(t,x=x,y=y)[0]
1696            phi = lambda t,x,y: vector_function(t,x=x,y=y)[1]
1697        else:
1698           # Assume info is in 2 keyword arguments
1699           if len(kwargs) == 2:
1700               s = kwargs['s']
1701               phi = kwargs['phi']
1702           else:
1703               raise Exception, 'Assumes two keyword arguments: s=..., phi=....'
1704
1705        self.speed = check_forcefield(s)
1706        self.phi = check_forcefield(phi)
1707
1708        self.const = eta_w*rho_a/rho_w
1709
1710    ##
1711    # @brief 'execute' this class instance.
1712    # @param domain
1713    def __call__(self, domain):
1714        """Evaluate windfield based on values found in domain"""
1715
1716        from math import pi, cos, sin, sqrt
1717
1718        xmom_update = domain.quantities['xmomentum'].explicit_update
1719        ymom_update = domain.quantities['ymomentum'].explicit_update
1720
1721        N = len(domain)    # number_of_triangles
1722        t = domain.time
1723
1724        if callable(self.speed):
1725            xc = domain.get_centroid_coordinates()
1726            s_vec = self.speed(t, xc[:,0], xc[:,1])
1727        else:
1728            # Assume s is a scalar
1729            try:
1730                s_vec = self.speed * num.ones(N, num.float)
1731            except:
1732                msg = 'Speed must be either callable or a scalar: %s' %self.s
1733                raise msg
1734
1735        if callable(self.phi):
1736            xc = domain.get_centroid_coordinates()
1737            phi_vec = self.phi(t, xc[:,0], xc[:,1])
1738        else:
1739            # Assume phi is a scalar
1740
1741            try:
1742                phi_vec = self.phi * num.ones(N, num.float)
1743            except:
1744                msg = 'Angle must be either callable or a scalar: %s' %self.phi
1745                raise msg
1746
1747        assign_windfield_values(xmom_update, ymom_update,
1748                                s_vec, phi_vec, self.const)
1749
1750
1751##
1752# @brief Assign wind field values
1753# @param xmom_update
1754# @param ymom_update
1755# @param s_vec
1756# @param phi_vec
1757# @param const
1758def assign_windfield_values(xmom_update, ymom_update,
1759                            s_vec, phi_vec, const):
1760    """Python version of assigning wind field to update vectors.
1761    A c version also exists (for speed)
1762    """
1763
1764    from math import pi, cos, sin, sqrt
1765
1766    N = len(s_vec)
1767    for k in range(N):
1768        s = s_vec[k]
1769        phi = phi_vec[k]
1770
1771        # Convert to radians
1772        phi = phi*pi/180
1773
1774        # Compute velocity vector (u, v)
1775        u = s*cos(phi)
1776        v = s*sin(phi)
1777
1778        # Compute wind stress
1779        S = const * sqrt(u**2 + v**2)
1780        xmom_update[k] += S*u
1781        ymom_update[k] += S*v
1782
1783
1784##
1785# @brief A class for a general explicit forcing term.
1786class General_forcing:
1787    """General explicit forcing term for update of quantity
1788
1789    This is used by Inflow and Rainfall for instance
1790
1791
1792    General_forcing(quantity_name, rate, center, radius, polygon)
1793
1794    domain:     ANUGA computational domain
1795    quantity_name: Name of quantity to update.
1796                   It must be a known conserved quantity.
1797
1798    rate [?/s]: Total rate of change over the specified area.
1799                This parameter can be either a constant or a
1800                function of time. Positive values indicate increases,
1801                negative values indicate decreases.
1802                Rate can be None at initialisation but must be specified
1803                before forcing term is applied (i.e. simulation has started).
1804
1805    center [m]: Coordinates at center of flow point
1806    radius [m]: Size of circular area
1807    polygon:    Arbitrary polygon
1808    default_rate: Rate to be used if rate fails (e.g. if model time exceeds its data)
1809                  Admissible types: None, constant number or function of t
1810
1811
1812    Either center, radius or polygon can be specified but not both.
1813    If neither are specified the entire domain gets updated.
1814    All coordinates to be specified in absolute UTM coordinates (x, y) assuming the zone of domain.
1815
1816    Inflow or Rainfall for examples of use
1817    """
1818
1819
1820    # FIXME (AnyOne) : Add various methods to allow spatial variations
1821
1822    ##
1823    # @brief Create an instance of this forcing term.
1824    # @param domain
1825    # @param quantity_name
1826    # @param rate
1827    # @param center
1828    # @param radius
1829    # @param polygon
1830    # @param default_rate
1831    # @param verbose
1832    def __init__(self,
1833                 domain,
1834                 quantity_name,
1835                 rate=0.0,
1836                 center=None,
1837                 radius=None,
1838                 polygon=None,
1839                 default_rate=None,
1840                 verbose=False):
1841
1842        from math import pi, cos, sin
1843
1844        if center is None:
1845            msg = 'I got radius but no center.'
1846            assert radius is None, msg
1847
1848        if radius is None:
1849            msg += 'I got center but no radius.'
1850            assert center is None, msg
1851
1852        self.domain = domain
1853        self.quantity_name = quantity_name
1854        self.rate = rate
1855        self.center = ensure_numeric(center)
1856        self.radius = radius
1857        self.polygon = polygon
1858        self.verbose = verbose
1859        self.value = 0.0    # Can be used to remember value at
1860                            # previous timestep in order to obtain rate
1861
1862        # Get boundary (in absolute coordinates)
1863        bounding_polygon = domain.get_boundary_polygon()
1864
1865        # Update area if applicable
1866        self.exchange_area = None
1867        if center is not None and radius is not None:
1868            assert len(center) == 2
1869            msg = 'Polygon cannot be specified when center and radius are'
1870            assert polygon is None, msg
1871
1872            self.exchange_area = radius**2*pi
1873
1874            # Check that circle center lies within the mesh.
1875            msg = 'Center %s specified for forcing term did not' % str(center)
1876            msg += 'fall within the domain boundary.'
1877            assert is_inside_polygon(center, bounding_polygon), msg
1878
1879            # Check that circle periphery lies within the mesh.
1880            N = 100
1881            periphery_points = []
1882            for i in range(N):
1883                theta = 2*pi*i/100
1884
1885                x = center[0] + radius*cos(theta)
1886                y = center[1] + radius*sin(theta)
1887
1888                periphery_points.append([x,y])
1889
1890            for point in periphery_points:
1891                msg = 'Point %s on periphery for forcing term' % str(point)
1892                msg += ' did not fall within the domain boundary.'
1893                assert is_inside_polygon(point, bounding_polygon), msg
1894
1895        if polygon is not None:
1896            # Check that polygon lies within the mesh.
1897            for point in self.polygon:
1898                msg = 'Point %s in polygon for forcing term' % str(point)
1899                msg += ' did not fall within the domain boundary.'
1900                assert is_inside_polygon(point, bounding_polygon), msg
1901
1902            # Compute area and check that it is greater than 0
1903            self.exchange_area = polygon_area(self.polygon)
1904
1905            msg = 'Polygon %s in forcing term' % str(self.polygon)
1906            msg += ' has area = %f' % self.exchange_area
1907            assert self.exchange_area > 0.0
1908
1909        # Pointer to update vector
1910        self.update = domain.quantities[self.quantity_name].explicit_update
1911
1912        # Determine indices in flow area
1913        N = len(domain)
1914        points = domain.get_centroid_coordinates(absolute=True)
1915
1916        # Calculate indices in exchange area for this forcing term
1917        self.exchange_indices = None
1918        if self.center is not None and self.radius is not None:
1919            # Inlet is circular
1920            inlet_region = 'center=%s, radius=%s' % (self.center, self.radius)
1921
1922            self.exchange_indices = []
1923            for k in range(N):
1924                x, y = points[k,:]    # Centroid
1925
1926                c = self.center
1927                if ((x-c[0])**2+(y-c[1])**2) < self.radius**2:
1928                    self.exchange_indices.append(k)
1929
1930        if self.polygon is not None:
1931            # Inlet is polygon
1932            inlet_region = ('polygon=\n%s, area=%f m^2' %
1933                            (self.polygon, self.exchange_area))
1934
1935            self.exchange_indices = inside_polygon(points, self.polygon)
1936
1937        if self.exchange_indices is not None:
1938            if len(self.exchange_indices) == 0:
1939                msg = 'No triangles have been identified in '
1940                msg += 'specified region: %s' % inlet_region
1941                raise Exception, msg
1942
1943        # Check and store default_rate
1944        msg = ('Keyword argument default_rate must be either None '
1945               'or a function of time.\nI got %s.' % str(default_rate))
1946        assert (default_rate is None or
1947                type(default_rate) in [IntType, FloatType] or
1948                callable(default_rate)), msg
1949
1950        if default_rate is not None:
1951            # If it is a constant, make it a function
1952            if not callable(default_rate):
1953                tmp = default_rate
1954                default_rate = lambda t: tmp
1955
1956            # Check that default_rate is a function of one argument
1957            try:
1958                default_rate(0.0)
1959            except:
1960                raise Exception, msg
1961
1962        self.default_rate = default_rate
1963        self.default_rate_invoked = False    # Flag
1964
1965    ##
1966    # @brief Execute this instance.
1967    # @param domain
1968    def __call__(self, domain):
1969        """Apply inflow function at time specified in domain, update stage"""
1970
1971        # Call virtual method allowing local modifications
1972        t = domain.get_time()
1973        try:
1974            rate = self.update_rate(t)
1975        except Modeltime_too_early, e:
1976            raise Modeltime_too_early, e
1977        except Modeltime_too_late, e:
1978            if self.default_rate is None:
1979                raise Exception, e    # Reraise exception
1980            else:
1981                # Pass control to default rate function
1982                rate = self.default_rate(t)
1983
1984                if self.default_rate_invoked is False:
1985                    # Issue warning the first time
1986                    msg = ('%s\n'
1987                           'Instead I will use the default rate: %s\n'
1988                           'Note: Further warnings will be supressed'
1989                           % (str(e), str(self.default_rate)))
1990                    warn(msg)
1991
1992                    # FIXME (Ole): Replace this crude flag with
1993                    # Python's ability to print warnings only once.
1994                    # See http://docs.python.org/lib/warning-filter.html
1995                    self.default_rate_invoked = True
1996
1997        if rate is None:
1998            msg = ('Attribute rate must be specified in General_forcing '
1999                   'or its descendants before attempting to call it')
2000            raise Exception, msg
2001
2002        # Now rate is a number
2003        if self.verbose is True:
2004            print 'Rate of %s at time = %.2f = %f' % (self.quantity_name,
2005                                                      domain.get_time(),
2006                                                      rate)
2007
2008        if self.exchange_indices is None:
2009            self.update[:] += rate
2010        else:
2011            # Brute force assignment of restricted rate
2012            for k in self.exchange_indices:
2013                self.update[k] += rate
2014
2015    ##
2016    # @brief Update the internal rate.
2017    # @param t A callable or scalar used to set the rate.
2018    # @return The new rate.
2019    def update_rate(self, t):
2020        """Virtual method allowing local modifications by writing an
2021        overriding version in descendant
2022        """
2023
2024        if callable(self.rate):
2025            rate = self.rate(t)
2026        else:
2027            rate = self.rate
2028
2029        return rate
2030
2031    ##
2032    # @brief Get values for the specified quantity.
2033    # @param quantity_name Name of the quantity of interest.
2034    # @return The value(s) of the quantity.
2035    # @note If 'quantity_name' is None, use self.quantity_name.
2036    def get_quantity_values(self, quantity_name=None):
2037        """Return values for specified quantity restricted to opening
2038
2039        Optionally a quantity name can be specified if values from another
2040        quantity is sought
2041        """
2042
2043        if quantity_name is None:
2044            quantity_name = self.quantity_name
2045
2046        q = self.domain.quantities[quantity_name]
2047        return q.get_values(location='centroids',
2048                            indices=self.exchange_indices)
2049
2050    ##
2051    # @brief Set value for the specified quantity.
2052    # @param val The value object used to set value.
2053    # @param quantity_name Name of the quantity of interest.
2054    # @note If 'quantity_name' is None, use self.quantity_name.
2055    def set_quantity_values(self, val, quantity_name=None):
2056        """Set values for specified quantity restricted to opening
2057
2058        Optionally a quantity name can be specified if values from another
2059        quantity is sought
2060        """
2061
2062        if quantity_name is None:
2063            quantity_name = self.quantity_name
2064
2065        q = self.domain.quantities[self.quantity_name]
2066        q.set_values(val,
2067                     location='centroids',
2068                     indices=self.exchange_indices)
2069
2070
2071##
2072# @brief A class for rainfall forcing function.
2073# @note Inherits from General_forcing.
2074class Rainfall(General_forcing):
2075    """Class Rainfall - general 'rain over entire domain' forcing term.
2076
2077    Used for implementing Rainfall over the entire domain.
2078
2079        Current Limited to only One Gauge..
2080
2081        Need to add Spatial Varying Capability
2082        (This module came from copying and amending the Inflow Code)
2083
2084    Rainfall(rain)
2085
2086    domain
2087    rain [mm/s]:  Total rain rate over the specified domain.
2088                  NOTE: Raingauge Data needs to reflect the time step.
2089                  IE: if Gauge is mm read at a time step, then the input
2090                  here is as mm/(timeStep) so 10mm in 5minutes becomes
2091                  10/(5x60) = 0.0333mm/s.
2092
2093                  This parameter can be either a constant or a
2094                  function of time. Positive values indicate inflow,
2095                  negative values indicate outflow.
2096                  (and be used for Infiltration - Write Seperate Module)
2097                  The specified flow will be divided by the area of
2098                  the inflow region and then applied to update the
2099                  stage quantity.
2100
2101    polygon: Specifies a polygon to restrict the rainfall.
2102
2103    Examples
2104    How to put them in a run File...
2105
2106    #------------------------------------------------------------------------
2107    # Setup specialised forcing terms
2108    #------------------------------------------------------------------------
2109    # This is the new element implemented by Ole and Rudy to allow direct
2110    # input of Rainfall in mm/s
2111
2112    catchmentrainfall = Rainfall(rain=file_function('Q100_2hr_Rain.tms'))
2113                        # Note need path to File in String.
2114                        # Else assumed in same directory
2115
2116    domain.forcing_terms.append(catchmentrainfall)
2117    """
2118
2119    ##
2120    # @brief Create an instance of the class.
2121    # @param domain Domain of interest.
2122    # @param rate Total rain rate over the specified domain (mm/s).
2123    # @param center
2124    # @param radius
2125    # @param polygon Polygon  to restrict rainfall.
2126    # @param default_rate
2127    # @param verbose True if this instance is to be verbose.
2128    def __init__(self,
2129                 domain,
2130                 rate=0.0,
2131                 center=None,
2132                 radius=None,
2133                 polygon=None,
2134                 default_rate=None,
2135                 verbose=False):
2136
2137        # Converting mm/s to m/s to apply in ANUGA)
2138        if callable(rate):
2139            rain = lambda t: rate(t)/1000.0
2140        else:
2141            rain = rate/1000.0
2142
2143        if default_rate is not None:
2144            if callable(default_rate):
2145                default_rain = lambda t: default_rate(t)/1000.0
2146            else:
2147                default_rain = default_rate/1000.0
2148        else:
2149            default_rain = None
2150
2151        # pass to parent class
2152        General_forcing.__init__(self,
2153                                 domain,
2154                                 'stage',
2155                                 rate=rain,
2156                                 center=center,
2157                                 radius=radius,
2158                                 polygon=polygon,
2159                                 default_rate=default_rain,
2160                                 verbose=verbose)
2161
2162
2163##
2164# @brief A class for inflow (rain and drain) forcing function.
2165# @note Inherits from General_forcing.
2166class Inflow(General_forcing):
2167    """Class Inflow - general 'rain and drain' forcing term.
2168
2169    Useful for implementing flows in and out of the domain.
2170
2171    Inflow(flow, center, radius, polygon)
2172
2173    domain
2174    rate [m^3/s]: Total flow rate over the specified area.
2175                  This parameter can be either a constant or a
2176                  function of time. Positive values indicate inflow,
2177                  negative values indicate outflow.
2178                  The specified flow will be divided by the area of
2179                  the inflow region and then applied to update stage.
2180    center [m]: Coordinates at center of flow point
2181    radius [m]: Size of circular area
2182    polygon:    Arbitrary polygon.
2183
2184    Either center, radius or polygon must be specified
2185
2186    Examples
2187
2188    # Constant drain at 0.003 m^3/s.
2189    # The outflow area is 0.07**2*pi=0.0154 m^2
2190    # This corresponds to a rate of change of 0.003/0.0154 = 0.2 m/s
2191    #
2192    Inflow((0.7, 0.4), 0.07, -0.003)
2193
2194
2195    # Tap turning up to a maximum inflow of 0.0142 m^3/s.
2196    # The inflow area is 0.03**2*pi = 0.00283 m^2
2197    # This corresponds to a rate of change of 0.0142/0.00283 = 5 m/s
2198    # over the specified area
2199    Inflow((0.5, 0.5), 0.03, lambda t: min(0.01*t, 0.0142))
2200
2201
2202    #------------------------------------------------------------------------
2203    # Setup specialised forcing terms
2204    #------------------------------------------------------------------------
2205    # This is the new element implemented by Ole to allow direct input
2206    # of Inflow in m^3/s
2207
2208    hydrograph = Inflow(center=(320, 300), radius=10,
2209                        rate=file_function('Q/QPMF_Rot_Sub13.tms'))
2210
2211    domain.forcing_terms.append(hydrograph)
2212    """
2213
2214    ##
2215    # @brief Create an instance of the class.
2216    # @param domain Domain of interest.
2217    # @param rate Total rain rate over the specified domain (mm/s).
2218    # @param center
2219    # @param radius
2220    # @param polygon Polygon  to restrict rainfall.
2221    # @param default_rate
2222    # @param verbose True if this instance is to be verbose.
2223    def __init__(self,
2224                 domain,
2225                 rate=0.0,
2226                 center=None,
2227                 radius=None,
2228                 polygon=None,
2229                 default_rate=None,
2230                 verbose=False):
2231        # Create object first to make area is available
2232        General_forcing.__init__(self,
2233                                 domain,
2234                                 'stage',
2235                                 rate=rate,
2236                                 center=center,
2237                                 radius=radius,
2238                                 polygon=polygon,
2239                                 default_rate=default_rate,
2240                                 verbose=verbose)
2241
2242    ##
2243    # @brief Update the instance rate.
2244    # @param t New rate object.
2245    def update_rate(self, t):
2246        """Virtual method allowing local modifications by writing an
2247        overriding version in descendant
2248
2249        This one converts m^3/s to m/s which can be added directly
2250        to 'stage' in ANUGA
2251        """
2252
2253        if callable(self.rate):
2254            _rate = self.rate(t)/self.exchange_area
2255        else:
2256            _rate = self.rate/self.exchange_area
2257
2258        return _rate
2259
2260
2261################################################################################
2262# Initialise module
2263################################################################################
2264
2265from anuga.utilities import compile
2266if compile.can_use_C_extension('shallow_water_ext.c'):
2267    # Underlying C implementations can be accessed
2268    from shallow_water_ext import rotate, assign_windfield_values
2269else:
2270    msg = 'C implementations could not be accessed by %s.\n ' % __file__
2271    msg += 'Make sure compile_all.py has been run as described in '
2272    msg += 'the ANUGA installation guide.'
2273    raise Exception, msg
2274
2275# Optimisation with psyco
2276from anuga.config import use_psyco
2277if use_psyco:
2278    try:
2279        import psyco
2280    except:
2281        import os
2282        if os.name == 'posix' and os.uname()[4] in ['x86_64', 'ia64']:
2283            pass
2284            #Psyco isn't supported on 64 bit systems, but it doesn't matter
2285        else:
2286            msg = ('WARNING: psyco (speedup) could not be imported, '
2287                   'you may want to consider installing it')
2288            print msg
2289    else:
2290        psyco.bind(Domain.distribute_to_vertices_and_edges)
2291        psyco.bind(Domain.compute_fluxes)
2292
2293
2294if __name__ == "__main__":
2295    pass
Note: See TracBrowser for help on using the repository browser.