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

Last change on this file since 7711 was 7711, checked in by hudson, 14 years ago

Refactored geometry classes to live in their own folder.

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