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

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

numpy changes.

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