source: anuga_core/source/anuga/shallow_water/shallow_water_domain.py @ 7312

Last change on this file since 7312 was 7312, checked in by ole, 15 years ago

Clarified uniquely stored vertices

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