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

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

Back-merge from Numeric to numpy.

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