source: trunk/anuga_core/source/anuga/shallow_water/shallow_water_domain.py @ 7810

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

Added aggressive psyco optimisation, fixed benchmark app to work with new API.

  • Property svn:keywords set to LastChangedDate LastChangedRevision LastChangedBy HeadURL Id
File size: 49.7 KB
RevLine 
[7276]1"""
2Finite-volume computations of the shallow water wave equation.
[4004]3
[4005]4Title: ANGUA shallow_water_domain - 2D triangular domains for finite-volume
5       computations of the shallow water wave equation.
[3804]6
7
[5186]8Author: Ole Nielsen, Ole.Nielsen@ga.gov.au
9        Stephen Roberts, Stephen.Roberts@anu.edu.au
10        Duncan Gray, Duncan.Gray@ga.gov.au
[4004]11
12CreationDate: 2004
13
14Description:
[4005]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
[4004]18
[4005]19    U_t + E_x + G_y = S
[3804]20
[4005]21    where
[3804]22
[4005]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, ...)
[3804]28
[4005]29    and _t, _x, _y denote the derivative with respect to t, x and y
30    respectively.
[3804]31
32
[4005]33    The quantities are
[3804]34
[4005]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]
[3804]45
[4005]46    eta                        mannings friction coefficient [to appear]
47    nu                         wind stress coefficient [to appear]
[3804]48
[4005]49    The conserved quantities are w, uh, vh
[3804]50
[4004]51Reference:
[4005]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
[3804]55
[7276]56    Hydrodynamic modelling of coastal inundation.
[4005]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
[3804]62
63
[4005]64SeeAlso:
65    TRAC administration of ANUGA (User Manuals etc) at
66    https://datamining.anu.edu.au/anuga and Subversion repository at
67    $HeadURL: trunk/anuga_core/source/anuga/shallow_water/shallow_water_domain.py $
[3804]68
[4004]69Constraints: See GPL license in the user guide
70Version: 1.0 ($Revision: 7810 $)
71ModifiedBy:
[4005]72    $Author: hudson $
73    $Date: 2010-06-08 04:31:04 +0000 (Tue, 08 Jun 2010) $
[3804]74"""
75
[4769]76# Subversion keywords:
[3804]77#
[4769]78# $LastChangedDate: 2010-06-08 04:31:04 +0000 (Tue, 08 Jun 2010) $
79# $LastChangedRevision: 7810 $
80# $LastChangedBy: hudson $
[3804]81
82
[7276]83import numpy as num
84
[7736]85from anuga.abstract_2d_finite_volumes.generic_domain \
86                    import Generic_Domain
[6928]87
[7733]88from anuga.shallow_water.forcing import Cross_section
[6178]89from anuga.pmesh.mesh_interface import create_mesh_from_regions
[5294]90from anuga.utilities.numerical_tools import gradient, mean, ensure_numeric
[5730]91from anuga.geospatial_data.geospatial_data import ensure_geospatial
[7776]92from anuga.file.sww import SWW_file 
[6086]93from anuga.config import netcdf_mode_r, netcdf_mode_w, netcdf_mode_a
[3804]94
[7276]95from anuga.fit_interpolate.interpolate import Modeltime_too_late, \
96                                              Modeltime_too_early
[5290]97
[7711]98from anuga.geometry.polygon import inside_polygon, polygon_area, \
[7276]99                                    is_inside_polygon
[7317]100import anuga.utilities.log as log
[5294]101
[7342]102import types
[5873]103from types import IntType, FloatType
[5874]104from warnings import warn
[5294]105
106
[7276]107################################################################################
[4769]108# Shallow water domain
[7276]109################################################################################
110
111##
112# @brief Class for a shallow water domain.
[3804]113class Domain(Generic_Domain):
114
[7342]115    #conserved_quantities = ['stage', 'xmomentum', 'ymomentum']
116    #other_quantities = ['elevation', 'friction']
[7276]117
118    ##
119    # @brief Instantiate a shallow water domain.
120    # @param coordinates
121    # @param vertices
122    # @param boundary
123    # @param tagged_elements
124    # @param geo_reference
125    # @param use_inscribed_circle
126    # @param mesh_filename
127    # @param use_cache
128    # @param verbose
[7519]129    # @param evolved_quantities
[7276]130    # @param full_send_dict
131    # @param ghost_recv_dict
132    # @param processor
133    # @param numproc
134    # @param number_of_full_nodes
135    # @param number_of_full_triangles
[3804]136    def __init__(self,
137                 coordinates=None,
138                 vertices=None,
139                 boundary=None,
140                 tagged_elements=None,
141                 geo_reference=None,
142                 use_inscribed_circle=False,
143                 mesh_filename=None,
144                 use_cache=False,
145                 verbose=False,
[7573]146                 conserved_quantities = None,
[7519]147                 evolved_quantities = None,
148                 other_quantities = None,
[3804]149                 full_send_dict=None,
150                 ghost_recv_dict=None,
151                 processor=0,
[3926]152                 numproc=1,
[3928]153                 number_of_full_nodes=None,
154                 number_of_full_triangles=None):
[3804]155
[7573]156        # Define quantities for the shallow_water domain
157        if conserved_quantities == None:
158            conserved_quantities = ['stage', 'xmomentum', 'ymomentum']
[7519]159
[7573]160        if evolved_quantities == None:
161            evolved_quantities =  ['stage', 'xmomentum', 'ymomentum']
162           
[7519]163        if other_quantities == None:
164            other_quantities = ['elevation', 'friction']
[7342]165       
[3804]166        Generic_Domain.__init__(self,
167                                coordinates,
168                                vertices,
169                                boundary,
[7342]170                                conserved_quantities,
[7573]171                                evolved_quantities,
[7342]172                                other_quantities,
[3804]173                                tagged_elements,
174                                geo_reference,
175                                use_inscribed_circle,
176                                mesh_filename,
177                                use_cache,
178                                verbose,
179                                full_send_dict,
180                                ghost_recv_dict,
181                                processor,
[3926]182                                numproc,
183                                number_of_full_nodes=number_of_full_nodes,
[7276]184                                number_of_full_triangles=number_of_full_triangles)
[3804]185
[7562]186        self.set_defaults()
187
188 
189        self.forcing_terms.append(manning_friction_implicit)
190        self.forcing_terms.append(gravity)
191
192        # Stored output
193        self.store = True
194        self.set_store_vertices_uniquely(False)
195
196        self.quantities_to_be_stored = {'elevation': 1, 
197                                        'stage': 2, 
198                                        'xmomentum': 2, 
199                                        'ymomentum': 2}
200
201
202
203
204    ##
205    # @brief Set default values, usually read in from a config file
206    # @param flag
207    def set_defaults(self):
208        """Set the default values in this routine. That way we can inherit class
209        and just over redefine the defaults for the new class
210        """
211
212        from anuga.config import minimum_storable_height
213        from anuga.config import minimum_allowed_height, maximum_allowed_speed
214        from anuga.config import g, epsilon, beta_w, beta_w_dry,\
215             beta_uh, beta_uh_dry, beta_vh, beta_vh_dry, tight_slope_limiters
216        from anuga.config import alpha_balance
217        from anuga.config import optimise_dry_cells
218        from anuga.config import optimised_gradient_limiter
219        from anuga.config import use_edge_limiter
220        from anuga.config import use_centroid_velocities
221
222
223
[7276]224        self.set_minimum_allowed_height(minimum_allowed_height)
[7562]225        self.maximum_allowed_speed = maximum_allowed_speed
[4769]226
[3804]227        self.g = g
[7276]228        self.beta_w = beta_w
229        self.beta_w_dry = beta_w_dry
230        self.beta_uh = beta_uh
[3804]231        self.beta_uh_dry = beta_uh_dry
[7276]232        self.beta_vh = beta_vh
[3804]233        self.beta_vh_dry = beta_vh_dry
[3876]234        self.alpha_balance = alpha_balance
[3804]235
[4631]236        self.tight_slope_limiters = tight_slope_limiters
[4685]237        self.optimise_dry_cells = optimise_dry_cells
[4239]238
[7562]239       
240        self.set_new_mannings_function(False)
[3804]241
242        self.minimum_storable_height = minimum_storable_height
[5162]243
[7562]244         # Limiters
[5175]245        self.use_edge_limiter = use_edge_limiter
[5162]246        self.optimised_gradient_limiter = optimised_gradient_limiter
[7562]247        self.use_centroid_velocities = use_centroid_velocities       
[3804]248
[7276]249    ##
250    # @brief
[7519]251    # @param flag
252    def set_new_mannings_function(self, flag=True):
253        """Cludge to allow unit test to pass, but to
254        also introduce new mannings friction function
255        which takes into account the slope of the bed.
256        The flag is tested in the python wrapper
257        mannings_friction_implicit
258        """
259        if flag:
260            self.use_new_mannings = True
261        else:
262            self.use_new_mannings = False
263
264
265    ##
266    # @brief
267    # @param flag
268    def set_use_edge_limiter(self, flag=True):
269        """Cludge to allow unit test to pass, but to
270        also introduce new edge limiting. The flag is
271        tested in distribute_to_vertices_and_edges
272        """
273        if flag:
274            self.use_edge_limiter = True
275        else:
276            self.use_edge_limiter = False
277
278
279         
280    ##
281    # @brief
[7276]282    # @param beta
[3848]283    def set_all_limiters(self, beta):
[7276]284        """Shorthand to assign one constant value [0,1] to all limiters.
[3847]285        0 Corresponds to first order, where as larger values make use of
[7276]286        the second order scheme.
[3847]287        """
[3804]288
[7276]289        self.beta_w = beta
290        self.beta_w_dry = beta
[5162]291        self.quantities['stage'].beta = beta
[7276]292
293        self.beta_uh = beta
[3847]294        self.beta_uh_dry = beta
[5162]295        self.quantities['xmomentum'].beta = beta
[7276]296
297        self.beta_vh = beta
[3847]298        self.beta_vh_dry = beta
[5162]299        self.quantities['ymomentum'].beta = beta
[3847]300
[7276]301    ##
302    # @brief
303    # @param flag
304    # @param reduction
[3804]305    def set_store_vertices_uniquely(self, flag, reduction=None):
306        """Decide whether vertex values should be stored uniquely as
[7312]307        computed in the model (True) or whether they should be reduced to one
308        value per vertex using self.reduction (False).
[3804]309        """
[3954]310
[7312]311        # FIXME (Ole): how about using the word "continuous vertex values" or
312        # "continuous stage surface"
[3804]313        self.smooth = not flag
314
[4733]315        # Reduction operation for get_vertex_values
[3804]316        if reduction is None:
317            self.reduction = mean
318            #self.reduction = min  #Looks better near steep slopes
319
[7276]320    ##
321    # @brief Set the minimum depth that will be written to an SWW file.
322    # @param minimum_storable_height The minimum stored height (in m).
[3804]323    def set_minimum_storable_height(self, minimum_storable_height):
[7276]324        """Set the minimum depth that will be recognised when writing
[3804]325        to an sww file. This is useful for removing thin water layers
326        that seems to be caused by friction creep.
327
328        The minimum allowed sww depth is in meters.
329        """
[7276]330
[3804]331        self.minimum_storable_height = minimum_storable_height
[4258]332
[7276]333    ##
334    # @brief
335    # @param minimum_allowed_height
[4258]336    def set_minimum_allowed_height(self, minimum_allowed_height):
[7276]337        """Set minimum depth that will be recognised in the numerical scheme.
[4258]338
339        The minimum allowed depth is in meters.
340
341        The parameter H0 (Minimal height for flux computation)
342        is also set by this function
343        """
[4438]344
345        #FIXME (Ole): rename H0 to minimum_allowed_height_in_flux_computation
[4701]346
347        #FIXME (Ole): Maybe use histogram to identify isolated extreme speeds
348        #and deal with them adaptively similarly to how we used to use 1 order
349        #steps to recover.
[7276]350
[4258]351        self.minimum_allowed_height = minimum_allowed_height
[7276]352        self.H0 = minimum_allowed_height
[3804]353
[7276]354    ##
355    # @brief
356    # @param maximum_allowed_speed
[3804]357    def set_maximum_allowed_speed(self, maximum_allowed_speed):
[7276]358        """Set the maximum particle speed that is allowed in water
[3804]359        shallower than minimum_allowed_height. This is useful for
360        controlling speeds in very thin layers of water and at the same time
361        allow some movement avoiding pooling of water.
[7276]362        """
[3804]363
364        self.maximum_allowed_speed = maximum_allowed_speed
365
[7276]366    ##
367    # @brief
368    # @param points_file_block_line_size
369    def set_points_file_block_line_size(self, points_file_block_line_size):
370        """Set the minimum depth that will be recognised when writing
[4254]371        to an sww file. This is useful for removing thin water layers
372        that seems to be caused by friction creep.
373
374        The minimum allowed sww depth is in meters.
375        """
376        self.points_file_block_line_size = points_file_block_line_size
[7276]377
[7342]378       
379    # FIXME: Probably obsolete in its curren form   
[7276]380    ##
381    # @brief Set the quantities that will be written to an SWW file.
382    # @param q The quantities to be written.
383    # @note Param 'q' may be None, single quantity or list of quantity strings.
384    # @note If 'q' is None, no quantities will be stored in the SWW file.
[3804]385    def set_quantities_to_be_stored(self, q):
[7342]386        """Specify which quantities will be stored in the sww file
387       
[3804]388        q must be either:
[7342]389          - a dictionary with quantity names
390          - a list of quantity names (for backwards compatibility)
[3804]391          - None
392
[7342]393        The format of the dictionary is as follows
394       
395        quantity_name: flag where flag must be either 1 or 2.
396        If flag is 1, the quantity is considered static and will be
397        stored once at the beginning of the simulation in a 1D array.
398       
399        If flag is 2, the quantity is considered time dependent and
400        it will be stored at each yieldstep by appending it to the
401        appropriate 2D array in the sww file.   
402       
[3804]403        If q is None, storage will be switched off altogether.
[7342]404       
405        Once the simulation has started and thw sww file opened,
406        this function will have no effect.
407       
408        The format, where q is a list of names is for backwards compatibility only.
409        It will take the specified quantities to be time dependent and assume
410        'elevation' to be static regardless.
[3804]411        """
412
413        if q is None:
[7342]414            self.quantities_to_be_stored = {}
[3804]415            self.store = False
416            return
417
[4701]418        # Check correcness
[3804]419        for quantity_name in q:
[7276]420            msg = ('Quantity %s is not a valid conserved quantity'
421                   % quantity_name)
[7342]422            assert quantity_name in self.quantities, msg
[3804]423
[7342]424        if type(q) == types.ListType:
425
426            msg = 'List arguments to set_quantities_to_be_stored '
427            msg += 'has been deprecated and will be removed in future '
428            msg += 'versions of ANUGA.'
429            msg += 'Please use dictionary argument instead'
430            warn(msg, DeprecationWarning) 
431
432       
433       
434            # FIXME: Raise deprecation warning
435            tmp = {}
436            for x in q:
437                tmp[x] = 2
438            tmp['elevation'] = 1   
439            q = tmp     
440           
441        assert type(q) == types.DictType   
[3804]442        self.quantities_to_be_stored = q
443
[7276]444    ##
445    # @brief
446    # @param indices
[3804]447    def get_wet_elements(self, indices=None):
448        """Return indices for elements where h > minimum_allowed_height
449
450        Optional argument:
451            indices is the set of element ids that the operation applies to.
452
453        Usage:
454            indices = get_wet_elements()
455
[7276]456        Note, centroid values are used for this operation
[3804]457        """
458
459        # Water depth below which it is considered to be 0 in the model
460        # FIXME (Ole): Allow this to be specified as a keyword argument as well
461        from anuga.config import minimum_allowed_height
462
463        elevation = self.get_quantity('elevation').\
[7276]464                        get_values(location='centroids', indices=indices)
465        stage = self.get_quantity('stage').\
[3804]466                    get_values(location='centroids', indices=indices)
467        depth = stage - elevation
468
469        # Select indices for which depth > 0
[6157]470        wet_indices = num.compress(depth > minimum_allowed_height,
471                                   num.arange(len(depth)))
[7276]472        return wet_indices
[3804]473
[7276]474    ##
475    # @brief
476    # @param indices
[3804]477    def get_maximum_inundation_elevation(self, indices=None):
478        """Return highest elevation where h > 0
479
480        Optional argument:
481            indices is the set of element ids that the operation applies to.
482
483        Usage:
484            q = get_maximum_inundation_elevation()
485
[7276]486        Note, centroid values are used for this operation
[3804]487        """
488
489        wet_elements = self.get_wet_elements(indices)
490        return self.get_quantity('elevation').\
[7276]491                   get_maximum_value(indices=wet_elements)
[3804]492
[7276]493    ##
494    # @brief
495    # @param indices
[3804]496    def get_maximum_inundation_location(self, indices=None):
[4554]497        """Return location of highest elevation where h > 0
[3804]498
499        Optional argument:
500            indices is the set of element ids that the operation applies to.
501
502        Usage:
[4554]503            q = get_maximum_inundation_location()
[3804]504
[7276]505        Note, centroid values are used for this operation
[3804]506        """
507
508        wet_elements = self.get_wet_elements(indices)
509        return self.get_quantity('elevation').\
[7276]510                   get_maximum_location(indices=wet_elements)
511
[7350]512
[7276]513    def get_flow_through_cross_section(self, polyline, verbose=False):
514        """Get the total flow through an arbitrary poly line.
515
516        This is a run-time equivalent of the function with same name
[7778]517        in sww_interrogate.py
[7276]518
[5729]519        Input:
[7276]520            polyline: Representation of desired cross section - it may contain
521                      multiple sections allowing for complex shapes. Assume
[5729]522                      absolute UTM coordinates.
[7276]523                      Format [[x0, y0], [x1, y1], ...]
524
525        Output:
[5729]526            Q: Total flow [m^3/s] across given segments.
[7276]527        """
528
[7350]529
530        cross_section = Cross_section(self, polyline, verbose)
531
532        return cross_section.get_flow_through_cross_section()
533
534
[7452]535    def get_energy_through_cross_section(self, polyline,
[7352]536                                         kind='total',
537                                         verbose=False):
538        """Obtain average energy head [m] across specified cross section.
539
540        Inputs:
541            polyline: Representation of desired cross section - it may contain
542                      multiple sections allowing for complex shapes. Assume
543                      absolute UTM coordinates.
544                      Format [[x0, y0], [x1, y1], ...]
545            kind:     Select which energy to compute.
546                      Options are 'specific' and 'total' (default)
547
548        Output:
549            E: Average energy [m] across given segments for all stored times.
550
551        The average velocity is computed for each triangle intersected by
552        the polyline and averaged weighted by segment lengths.
553
554        The typical usage of this function would be to get average energy of
555        flow in a channel, and the polyline would then be a cross section
556        perpendicular to the flow.
557
558        #FIXME (Ole) - need name for this energy reflecting that its dimension
559        is [m].
560        """
561
562
563
564        cross_section = Cross_section(self, polyline, verbose)
565
[7452]566        return cross_section.get_energy_through_cross_section(kind)
[7352]567
568
569    ##
[7276]570    # @brief Run integrity checks on shallow water domain.
[3804]571    def check_integrity(self):
572        Generic_Domain.check_integrity(self)
573
574        #Check that we are solving the shallow water wave equation
575        msg = 'First conserved quantity must be "stage"'
576        assert self.conserved_quantities[0] == 'stage', msg
577        msg = 'Second conserved quantity must be "xmomentum"'
578        assert self.conserved_quantities[1] == 'xmomentum', msg
579        msg = 'Third conserved quantity must be "ymomentum"'
580        assert self.conserved_quantities[2] == 'ymomentum', msg
581
[7276]582    ##
583    # @brief
[3804]584    def extrapolate_second_order_sw(self):
[7276]585        #Call correct module function (either from this module or C-extension)
[3804]586        extrapolate_second_order_sw(self)
587
[7276]588    ##
589    # @brief
[3804]590    def compute_fluxes(self):
[7276]591        #Call correct module function (either from this module or C-extension)
[3804]592        compute_fluxes(self)
593
[7276]594    ##
595    # @brief
[3804]596    def distribute_to_vertices_and_edges(self):
[4769]597        # Call correct module function
[5176]598        if self.use_edge_limiter:
[7276]599            distribute_using_edge_limiter(self)
[5175]600        else:
[5306]601            distribute_using_vertex_limiter(self)
[3804]602
[7352]603
604
[7276]605    ##
606    # @brief Evolve the model by one step.
607    # @param yieldstep
608    # @param finaltime
609    # @param duration
610    # @param skip_initial_step
[3804]611    def evolve(self,
[7276]612               yieldstep=None,
613               finaltime=None,
614               duration=None,
615               skip_initial_step=False):
616        """Specialisation of basic evolve method from parent class"""
[3804]617
[4769]618        # Call check integrity here rather than from user scripts
619        # self.check_integrity()
[3804]620
[7276]621        msg = 'Attribute self.beta_w must be in the interval [0, 2]'
[5162]622        assert 0 <= self.beta_w <= 2.0, msg
[3804]623
[4769]624        # Initial update of vertex and edge values before any STORAGE
[7276]625        # and or visualisation.
[4769]626        # This is done again in the initialisation of the Generic_Domain
[7276]627        # evolve loop but we do it here to ensure the values are ok for storage.
[3804]628        self.distribute_to_vertices_and_edges()
629
630        if self.store is True and self.time == 0.0:
631            self.initialise_storage()
632
[4769]633        # Call basic machinery from parent class
[7276]634        for t in Generic_Domain.evolve(self, yieldstep=yieldstep,
635                                       finaltime=finaltime, duration=duration,
[3804]636                                       skip_initial_step=skip_initial_step):
[4769]637            # Store model data, e.g. for subsequent visualisation
[3804]638            if self.store is True:
[7340]639                self.store_timestep()
[3804]640
[4769]641            # Pass control on to outer loop for more specific actions
[3804]642            yield(t)
643
[7276]644    ##
645    # @brief
[3804]646    def initialise_storage(self):
647        """Create and initialise self.writer object for storing data.
648        Also, save x,y and bed elevation
649        """
[7342]650       
[4769]651        # Initialise writer
[7340]652        self.writer = SWW_file(self)
[3804]653
[4769]654        # Store vertices and connectivity
[3804]655        self.writer.store_connectivity()
656
[7276]657    ##
658    # @brief
659    # @param name
[7340]660    def store_timestep(self):
661        """Store time dependent quantities and time.
[3804]662
663        Precondition:
[7340]664           self.writer has been initialised
[3804]665        """
[7276]666
[7340]667        self.writer.store_timestep()
[3804]668
[7276]669    ##
670    # @brief Get time stepping statistics string for printing.
671    # @param track_speeds
672    # @param triangle_id
[4836]673    def timestepping_statistics(self,
674                                track_speeds=False,
[7276]675                                triangle_id=None):
[4827]676        """Return string with time stepping statistics for printing or logging
[3804]677
[4827]678        Optional boolean keyword track_speeds decides whether to report
679        location of smallest timestep as well as a histogram and percentile
680        report.
681        """
682
[7276]683        from anuga.config import epsilon, g
[4827]684
685        # Call basic machinery from parent class
[7276]686        msg = Generic_Domain.timestepping_statistics(self, track_speeds,
[4836]687                                                     triangle_id)
[4827]688
689        if track_speeds is True:
690            # qwidth determines the text field used for quantities
691            qwidth = self.qwidth
[7276]692
[4836]693            # Selected triangle
[4827]694            k = self.k
695
696            # Report some derived quantities at vertices, edges and centroid
697            # specific to the shallow water wave equation
698            z = self.quantities['elevation']
[7276]699            w = self.quantities['stage']
[4827]700
701            Vw = w.get_values(location='vertices', indices=[k])[0]
702            Ew = w.get_values(location='edges', indices=[k])[0]
703            Cw = w.get_values(location='centroids', indices=[k])
704
705            Vz = z.get_values(location='vertices', indices=[k])[0]
706            Ez = z.get_values(location='edges', indices=[k])[0]
[7276]707            Cz = z.get_values(location='centroids', indices=[k])
[4827]708
709            name = 'depth'
710            Vh = Vw-Vz
711            Eh = Ew-Ez
712            Ch = Cw-Cz
[7276]713
[4827]714            s  = '    %s: vertex_values =  %.4f,\t %.4f,\t %.4f\n'\
715                 %(name.ljust(qwidth), Vh[0], Vh[1], Vh[2])
[7276]716
[4827]717            s += '    %s: edge_values =    %.4f,\t %.4f,\t %.4f\n'\
718                 %(name.ljust(qwidth), Eh[0], Eh[1], Eh[2])
[7276]719
[4827]720            s += '    %s: centroid_value = %.4f\n'\
[7276]721                 %(name.ljust(qwidth), Ch[0])
722
[4827]723            msg += s
724
725            uh = self.quantities['xmomentum']
726            vh = self.quantities['ymomentum']
727
728            Vuh = uh.get_values(location='vertices', indices=[k])[0]
729            Euh = uh.get_values(location='edges', indices=[k])[0]
730            Cuh = uh.get_values(location='centroids', indices=[k])
[7276]731
[4827]732            Vvh = vh.get_values(location='vertices', indices=[k])[0]
733            Evh = vh.get_values(location='edges', indices=[k])[0]
734            Cvh = vh.get_values(location='centroids', indices=[k])
735
736            # Speeds in each direction
737            Vu = Vuh/(Vh + epsilon)
738            Eu = Euh/(Eh + epsilon)
[7276]739            Cu = Cuh/(Ch + epsilon)
[4827]740            name = 'U'
741            s  = '    %s: vertex_values =  %.4f,\t %.4f,\t %.4f\n'\
742                 %(name.ljust(qwidth), Vu[0], Vu[1], Vu[2])
[7276]743
[4827]744            s += '    %s: edge_values =    %.4f,\t %.4f,\t %.4f\n'\
745                 %(name.ljust(qwidth), Eu[0], Eu[1], Eu[2])
[7276]746
[4827]747            s += '    %s: centroid_value = %.4f\n'\
[7276]748                 %(name.ljust(qwidth), Cu[0])
749
[4827]750            msg += s
751
752            Vv = Vvh/(Vh + epsilon)
753            Ev = Evh/(Eh + epsilon)
[7276]754            Cv = Cvh/(Ch + epsilon)
[4827]755            name = 'V'
756            s  = '    %s: vertex_values =  %.4f,\t %.4f,\t %.4f\n'\
757                 %(name.ljust(qwidth), Vv[0], Vv[1], Vv[2])
[7276]758
[4827]759            s += '    %s: edge_values =    %.4f,\t %.4f,\t %.4f\n'\
760                 %(name.ljust(qwidth), Ev[0], Ev[1], Ev[2])
[7276]761
[4827]762            s += '    %s: centroid_value = %.4f\n'\
[7276]763                 %(name.ljust(qwidth), Cv[0])
764
[4827]765            msg += s
766
767            # Froude number in each direction
768            name = 'Froude (x)'
[6157]769            Vfx = Vu/(num.sqrt(g*Vh) + epsilon)
770            Efx = Eu/(num.sqrt(g*Eh) + epsilon)
771            Cfx = Cu/(num.sqrt(g*Ch) + epsilon)
[7276]772
[4827]773            s  = '    %s: vertex_values =  %.4f,\t %.4f,\t %.4f\n'\
774                 %(name.ljust(qwidth), Vfx[0], Vfx[1], Vfx[2])
[7276]775
[4827]776            s += '    %s: edge_values =    %.4f,\t %.4f,\t %.4f\n'\
777                 %(name.ljust(qwidth), Efx[0], Efx[1], Efx[2])
[7276]778
[4827]779            s += '    %s: centroid_value = %.4f\n'\
[7276]780                 %(name.ljust(qwidth), Cfx[0])
781
[4827]782            msg += s
783
784            name = 'Froude (y)'
[6157]785            Vfy = Vv/(num.sqrt(g*Vh) + epsilon)
786            Efy = Ev/(num.sqrt(g*Eh) + epsilon)
787            Cfy = Cv/(num.sqrt(g*Ch) + epsilon)
[7276]788
[4827]789            s  = '    %s: vertex_values =  %.4f,\t %.4f,\t %.4f\n'\
790                 %(name.ljust(qwidth), Vfy[0], Vfy[1], Vfy[2])
[7276]791
[4827]792            s += '    %s: edge_values =    %.4f,\t %.4f,\t %.4f\n'\
793                 %(name.ljust(qwidth), Efy[0], Efy[1], Efy[2])
[7276]794
[4827]795            s += '    %s: centroid_value = %.4f\n'\
[7276]796                 %(name.ljust(qwidth), Cfy[0])
[4827]797
[7276]798            msg += s
[4827]799
800        return msg
[7276]801       
[4827]802       
803
[6654]804    def compute_boundary_flows(self):
805        """Compute boundary flows at current timestep.
[6647]806                       
[6648]807        Quantities computed are:
808           Total inflow across boundary
809           Total outflow across boundary
[6654]810           Flow across each tagged boundary segment
[6648]811        """
[6647]812               
[6648]813        # Run through boundary array and compute for each segment
814        # the normal momentum ((uh, vh) dot normal) times segment length.
[6653]815        # Based on sign accumulate this into boundary_inflow and boundary_outflow.
[6647]816                       
[6653]817        # Compute flows along boundary
818       
[6654]819        uh = self.get_quantity('xmomentum').get_values(location='edges')
820        vh = self.get_quantity('ymomentum').get_values(location='edges')       
[6653]821       
822        # Loop through edges that lie on the boundary and calculate
823        # flows
[6654]824        boundary_flows = {}
825        total_boundary_inflow = 0.0
826        total_boundary_outflow = 0.0
[6653]827        for vol_id, edge_id in self.boundary:
[6654]828            # Compute normal flow across edge. Since normal vector points
829            # away from triangle, a positive sign means that water
830            # flows *out* from this triangle.
831             
832            momentum = [uh[vol_id, edge_id], vh[vol_id, edge_id]]
833            normal = self.mesh.get_normal(vol_id, edge_id)
834            length = self.mesh.get_edgelength(vol_id, edge_id)           
835            normal_flow = num.dot(momentum, normal)*length
836           
837            # Reverse sign so that + is taken to mean inflow
838            # and - means outflow. This is more intuitive.
839            edge_flow = -normal_flow
840           
841            # Tally up inflows and outflows separately
842            if edge_flow > 0:
843                # Flow is inflow     
844                total_boundary_inflow += edge_flow                                 
845            else:
846                # Flow is outflow
847                total_boundary_outflow += edge_flow   
[6653]848
[6654]849            # Tally up flows by boundary tag
850            tag = self.boundary[(vol_id, edge_id)]
851           
852            if tag not in boundary_flows:
853                boundary_flows[tag] = 0.0
854            boundary_flows[tag] += edge_flow
855           
856               
857        return boundary_flows, total_boundary_inflow, total_boundary_outflow
[6653]858       
[6654]859
860    def compute_forcing_flows(self):
861        """Compute flows in and out of domain due to forcing terms.
862                       
863        Quantities computed are:
864               
[6653]865       
[6654]866           Total inflow through forcing terms
867           Total outflow through forcing terms
868           Current total volume in domain       
869
870        """
871
872        #FIXME(Ole): We need to separate what part of explicit_update was
873        # due to the normal flux calculations and what is due to forcing terms.
874       
875        pass
876                       
877       
878    def compute_total_volume(self):
879        """Compute total volume (m^3) of water in entire domain
880        """
881       
882        area = self.mesh.get_areas()
883        volume = 0.0
884       
885        stage = self.get_quantity('stage').get_values(location='centroids')
886        elevation = self.get_quantity('elevation').get_values(location='centroids')       
887        depth = stage-elevation
888       
889        return num.sum(depth*area)
890       
891       
892    def volumetric_balance_statistics(self):               
893        """Create volumetric balance report suitable for printing or logging.
894        """
895       
[7276]896        (boundary_flows, total_boundary_inflow,
897         total_boundary_outflow) = self.compute_boundary_flows() 
898       
[6654]899        s = '---------------------------\n'       
900        s += 'Volumetric balance report:\n'
901        s += '--------------------------\n'
[7276]902        s += 'Total boundary inflow [m^3/s]: %.2f\n' % total_boundary_inflow
903        s += 'Total boundary outflow [m^3/s]: %.2f\n' % total_boundary_outflow       
[6654]904        s += 'Net boundary flow by tags [m^3/s]\n'
905        for tag in boundary_flows:
[7276]906            s += '    %s [m^3/s]: %.2f\n' % (tag, boundary_flows[tag])
[6654]907       
[7276]908        s += 'Total net boundary flow [m^3/s]: %.2f\n' % (total_boundary_inflow + total_boundary_outflow) 
909        s += 'Total volume in domain [m^3]: %.2f\n' % self.compute_total_volume()
[6654]910       
911        # The go through explicit forcing update and record the rate of change for stage and
912        # record into forcing_inflow and forcing_outflow. Finally compute integral
913        # of depth to obtain total volume of domain.
914       
915        # FIXME(Ole): This part is not yet done.               
916       
917        return s       
918           
[7276]919################################################################################
920# End of class Shallow Water Domain
921################################################################################
[3804]922
[4769]923#-----------------
[3804]924# Flux computation
[4769]925#-----------------
[3804]926
[7276]927## @brief Compute fluxes and timestep suitable for all volumes in domain.
928# @param domain The domain to calculate fluxes for.
[3804]929def compute_fluxes(domain):
[7276]930    """Compute fluxes and timestep suitable for all volumes in domain.
[3804]931
932    Compute total flux for each conserved quantity using "flux_function"
933
934    Fluxes across each edge are scaled by edgelengths and summed up
935    Resulting flux is then scaled by area and stored in
936    explicit_update for each of the three conserved quantities
937    stage, xmomentum and ymomentum
938
939    The maximal allowable speed computed by the flux_function for each volume
940    is converted to a timestep that must not be exceeded. The minimum of
941    those is computed as the next overall timestep.
942
943    Post conditions:
944      domain.explicit_update is reset to computed flux values
945      domain.timestep is set to the largest step satisfying all volumes.
[4769]946
947    This wrapper calls the underlying C version of compute fluxes
[3804]948    """
949
950    import sys
[7276]951    from shallow_water_ext import compute_fluxes_ext_central \
952                                  as compute_fluxes_ext
[3804]953
[7276]954    N = len(domain)    # number_of_triangles
[3804]955
[4769]956    # Shortcuts
[3804]957    Stage = domain.quantities['stage']
958    Xmom = domain.quantities['xmomentum']
959    Ymom = domain.quantities['ymomentum']
960    Bed = domain.quantities['elevation']
961
962    timestep = float(sys.maxint)
963
[4769]964    flux_timestep = compute_fluxes_ext(timestep,
965                                       domain.epsilon,
966                                       domain.H0,
967                                       domain.g,
968                                       domain.neighbours,
969                                       domain.neighbour_edges,
970                                       domain.normals,
971                                       domain.edgelengths,
972                                       domain.radii,
973                                       domain.areas,
974                                       domain.tri_full_flag,
975                                       Stage.edge_values,
976                                       Xmom.edge_values,
977                                       Ymom.edge_values,
978                                       Bed.edge_values,
979                                       Stage.boundary_values,
980                                       Xmom.boundary_values,
981                                       Ymom.boundary_values,
982                                       Stage.explicit_update,
983                                       Xmom.explicit_update,
984                                       Ymom.explicit_update,
985                                       domain.already_computed_flux,
986                                       domain.max_speed,
987                                       int(domain.optimise_dry_cells))
[3804]988
[4769]989    domain.flux_timestep = flux_timestep
[3804]990
[7276]991################################################################################
[4769]992# Module functions for gradient limiting
[7276]993################################################################################
[3804]994
[7276]995##
996# @brief Wrapper for C version of extrapolate_second_order_sw.
997# @param domain The domain to operate on.
998# @note MH090605 The following method belongs to the shallow_water domain class
999#       see comments in the corresponding method in shallow_water_ext.c
1000def extrapolate_second_order_sw(domain):
1001    """Wrapper calling C version of extrapolate_second_order_sw"""
[3804]1002
1003    import sys
[7276]1004    from shallow_water_ext import extrapolate_second_order_sw as extrapol2
[3804]1005
[3928]1006    N = len(domain) # number_of_triangles
[3804]1007
[4710]1008    # Shortcuts
[3804]1009    Stage = domain.quantities['stage']
1010    Xmom = domain.quantities['xmomentum']
1011    Ymom = domain.quantities['ymomentum']
1012    Elevation = domain.quantities['elevation']
[4710]1013
[4769]1014    extrapol2(domain,
1015              domain.surrogate_neighbours,
1016              domain.number_of_boundaries,
1017              domain.centroid_coordinates,
1018              Stage.centroid_values,
1019              Xmom.centroid_values,
1020              Ymom.centroid_values,
1021              Elevation.centroid_values,
1022              domain.vertex_coordinates,
1023              Stage.vertex_values,
1024              Xmom.vertex_values,
1025              Ymom.vertex_values,
1026              Elevation.vertex_values,
[5315]1027              int(domain.optimise_dry_cells))
[3804]1028
[7276]1029##
1030# @brief Distribution from centroids to vertices specific to the SWW eqn.
1031# @param domain The domain to operate on.
[5306]1032def distribute_using_vertex_limiter(domain):
[7276]1033    """Distribution from centroids to vertices specific to the SWW equation.
[3804]1034
1035    It will ensure that h (w-z) is always non-negative even in the
1036    presence of steep bed-slopes by taking a weighted average between shallow
1037    and deep cases.
1038
1039    In addition, all conserved quantities get distributed as per either a
1040    constant (order==1) or a piecewise linear function (order==2).
1041
1042    FIXME: more explanation about removal of artificial variability etc
1043
1044    Precondition:
1045      All quantities defined at centroids and bed elevation defined at
1046      vertices.
1047
1048    Postcondition
1049      Conserved quantities defined at vertices
1050    """
1051
[4769]1052    # Remove very thin layers of water
[3804]1053    protect_against_infinitesimal_and_negative_heights(domain)
1054
[4769]1055    # Extrapolate all conserved quantities
[5162]1056    if domain.optimised_gradient_limiter:
[4769]1057        # MH090605 if second order,
1058        # perform the extrapolation and limiting on
1059        # all of the conserved quantities
[3804]1060
1061        if (domain._order_ == 1):
1062            for name in domain.conserved_quantities:
1063                Q = domain.quantities[name]
1064                Q.extrapolate_first_order()
1065        elif domain._order_ == 2:
1066            domain.extrapolate_second_order_sw()
1067        else:
1068            raise 'Unknown order'
1069    else:
[4769]1070        # Old code:
[3804]1071        for name in domain.conserved_quantities:
1072            Q = domain.quantities[name]
[4701]1073
[3804]1074            if domain._order_ == 1:
1075                Q.extrapolate_first_order()
1076            elif domain._order_ == 2:
[5306]1077                Q.extrapolate_second_order_and_limit_by_vertex()
[3804]1078            else:
1079                raise 'Unknown order'
1080
[5290]1081    # Take bed elevation into account when water heights are small
[3804]1082    balance_deep_and_shallow(domain)
1083
[5290]1084    # Compute edge values by interpolation
[3804]1085    for name in domain.conserved_quantities:
1086        Q = domain.quantities[name]
1087        Q.interpolate_from_vertices_to_edges()
1088
[7276]1089##
1090# @brief Distribution from centroids to edges specific to the SWW eqn.
1091# @param domain The domain to operate on.
[5306]1092def distribute_using_edge_limiter(domain):
[7276]1093    """Distribution from centroids to edges specific to the SWW eqn.
[5306]1094
1095    It will ensure that h (w-z) is always non-negative even in the
1096    presence of steep bed-slopes by taking a weighted average between shallow
1097    and deep cases.
1098
1099    In addition, all conserved quantities get distributed as per either a
1100    constant (order==1) or a piecewise linear function (order==2).
1101
1102
1103    Precondition:
1104      All quantities defined at centroids and bed elevation defined at
1105      vertices.
1106
1107    Postcondition
1108      Conserved quantities defined at vertices
1109    """
1110
1111    # Remove very thin layers of water
1112    protect_against_infinitesimal_and_negative_heights(domain)
1113
1114    for name in domain.conserved_quantities:
1115        Q = domain.quantities[name]
1116        if domain._order_ == 1:
1117            Q.extrapolate_first_order()
1118        elif domain._order_ == 2:
1119            Q.extrapolate_second_order_and_limit_by_edge()
1120        else:
1121            raise 'Unknown order'
1122
1123    balance_deep_and_shallow(domain)
1124
1125    # Compute edge values by interpolation
1126    for name in domain.conserved_quantities:
1127        Q = domain.quantities[name]
1128        Q.interpolate_from_vertices_to_edges()
1129
[7276]1130##
1131# @brief  Protect against infinitesimal heights and associated high velocities.
1132# @param domain The domain to operate on.
[3804]1133def protect_against_infinitesimal_and_negative_heights(domain):
[7276]1134    """Protect against infinitesimal heights and associated high velocities"""
[3804]1135
[7276]1136    from shallow_water_ext import protect
1137
[4769]1138    # Shortcuts
[3804]1139    wc = domain.quantities['stage'].centroid_values
1140    zc = domain.quantities['elevation'].centroid_values
1141    xmomc = domain.quantities['xmomentum'].centroid_values
1142    ymomc = domain.quantities['ymomentum'].centroid_values
1143
1144    protect(domain.minimum_allowed_height, domain.maximum_allowed_speed,
1145            domain.epsilon, wc, zc, xmomc, ymomc)
1146
[7276]1147##
1148# @brief Wrapper for C function balance_deep_and_shallow_c().
1149# @param domain The domain to operate on.
[3804]1150def balance_deep_and_shallow(domain):
1151    """Compute linear combination between stage as computed by
1152    gradient-limiters limiting using w, and stage computed by
1153    gradient-limiters limiting using h (h-limiter).
1154    The former takes precedence when heights are large compared to the
1155    bed slope while the latter takes precedence when heights are
1156    relatively small.  Anything in between is computed as a balanced
1157    linear combination in order to avoid numerical disturbances which
1158    would otherwise appear as a result of hard switching between
1159    modes.
1160
[4769]1161    Wrapper for C implementation
[3804]1162    """
1163
[7276]1164    from shallow_water_ext import balance_deep_and_shallow \
1165                                  as balance_deep_and_shallow_c
[5175]1166
[4733]1167    # Shortcuts
[3804]1168    wc = domain.quantities['stage'].centroid_values
1169    zc = domain.quantities['elevation'].centroid_values
1170    wv = domain.quantities['stage'].vertex_values
1171    zv = domain.quantities['elevation'].vertex_values
1172
[4733]1173    # Momentums at centroids
[3804]1174    xmomc = domain.quantities['xmomentum'].centroid_values
1175    ymomc = domain.quantities['ymomentum'].centroid_values
1176
[4733]1177    # Momentums at vertices
[3804]1178    xmomv = domain.quantities['xmomentum'].vertex_values
1179    ymomv = domain.quantities['ymomentum'].vertex_values
1180
[5442]1181    balance_deep_and_shallow_c(domain,
[7276]1182                               wc, zc, wv, zv, wc,
[5442]1183                               xmomc, ymomc, xmomv, ymomv)
[3804]1184
1185
1186
[7276]1187################################################################################
[4769]1188# Standard forcing terms
[7276]1189################################################################################
[4769]1190
[7276]1191##
1192# @brief Apply gravitational pull in the presence of bed slope.
1193# @param domain The domain to apply gravity to.
1194# @note Wrapper for C function gravity_c().
[3804]1195def gravity(domain):
1196    """Apply gravitational pull in the presence of bed slope
[4769]1197    Wrapper calls underlying C implementation
[3804]1198    """
1199
[7276]1200    from shallow_water_ext import gravity as gravity_c
1201
[7519]1202    xmom_update = domain.quantities['xmomentum'].explicit_update
1203    ymom_update = domain.quantities['ymomentum'].explicit_update
[3804]1204
[4687]1205    stage = domain.quantities['stage']
1206    elevation = domain.quantities['elevation']
[3804]1207
[4687]1208    h = stage.centroid_values - elevation.centroid_values
1209    z = elevation.vertex_values
1210
[3804]1211    x = domain.get_vertex_coordinates()
1212    g = domain.g
1213
[7519]1214    gravity_c(g, h, z, x, xmom_update, ymom_update)    #, 1.0e-6)
[3804]1215
[7276]1216##
1217# @brief Apply friction to a surface (implicit).
1218# @param domain The domain to apply Manning friction to.
1219# @note Wrapper for C function manning_friction_c().
[4769]1220def manning_friction_implicit(domain):
[7276]1221    """Apply (Manning) friction to water momentum
[4769]1222    Wrapper for c version
[3804]1223    """
1224
[7519]1225    from shallow_water_ext import manning_friction_old
1226    from shallow_water_ext import manning_friction_new
[3804]1227
1228    xmom = domain.quantities['xmomentum']
1229    ymom = domain.quantities['ymomentum']
1230
[7519]1231    x = domain.get_vertex_coordinates()
1232   
[3804]1233    w = domain.quantities['stage'].centroid_values
[7519]1234    z = domain.quantities['elevation'].vertex_values
[3804]1235
1236    uh = xmom.centroid_values
1237    vh = ymom.centroid_values
1238    eta = domain.quantities['friction'].centroid_values
1239
1240    xmom_update = xmom.semi_implicit_update
1241    ymom_update = ymom.semi_implicit_update
1242
[3928]1243    N = len(domain)
[3804]1244    eps = domain.minimum_allowed_height
1245    g = domain.g
1246
[7519]1247    if domain.use_new_mannings:
1248        manning_friction_new(g, eps, x, w, uh, vh, z, eta, xmom_update, ymom_update)
1249    else:
1250        manning_friction_old(g, eps, w, uh, vh, z, eta, xmom_update, ymom_update)
1251   
[3804]1252
[7276]1253##
1254# @brief Apply friction to a surface (explicit).
1255# @param domain The domain to apply Manning friction to.
1256# @note Wrapper for C function manning_friction_c().
[4769]1257def manning_friction_explicit(domain):
[7276]1258    """Apply (Manning) friction to water momentum
[4769]1259    Wrapper for c version
[3804]1260    """
1261
[7519]1262    from shallow_water_ext import manning_friction_old
1263    from shallow_water_ext import manning_friction_new
[3804]1264
1265    xmom = domain.quantities['xmomentum']
1266    ymom = domain.quantities['ymomentum']
1267
[7519]1268    x = domain.get_vertex_coordinates()
1269   
[3804]1270    w = domain.quantities['stage'].centroid_values
[7519]1271    z = domain.quantities['elevation'].vertex_values
[3804]1272
1273    uh = xmom.centroid_values
1274    vh = ymom.centroid_values
1275    eta = domain.quantities['friction'].centroid_values
1276
1277    xmom_update = xmom.explicit_update
1278    ymom_update = ymom.explicit_update
1279
[3928]1280    N = len(domain)
[3804]1281    eps = domain.minimum_allowed_height
1282    g = domain.g
1283
1284
[7519]1285    if domain.use_new_mannings:
1286        manning_friction_new(g, eps, x, w, uh, vh, z, eta, xmom_update, ymom_update)
1287    else:
1288        manning_friction_old(g, eps, w, uh, vh, z, eta, xmom_update, ymom_update)
[3804]1289
[7519]1290
1291
[4769]1292# FIXME (Ole): This was implemented for use with one of the analytical solutions (Sampson?)
[7276]1293##
1294# @brief Apply linear friction to a surface.
1295# @param domain The domain to apply Manning friction to.
1296# @note Is this still used (30 Oct 2007)?
[3804]1297def linear_friction(domain):
1298    """Apply linear friction to water momentum
1299
1300    Assumes quantity: 'linear_friction' to be present
1301    """
1302
1303    from math import sqrt
1304
1305    w = domain.quantities['stage'].centroid_values
1306    z = domain.quantities['elevation'].centroid_values
1307    h = w-z
1308
1309    uh = domain.quantities['xmomentum'].centroid_values
1310    vh = domain.quantities['ymomentum'].centroid_values
1311    tau = domain.quantities['linear_friction'].centroid_values
1312
1313    xmom_update = domain.quantities['xmomentum'].semi_implicit_update
1314    ymom_update = domain.quantities['ymomentum'].semi_implicit_update
1315
[3928]1316    N = len(domain) # number_of_triangles
[3804]1317    eps = domain.minimum_allowed_height
1318    g = domain.g #Not necessary? Why was this added?
1319
1320    for k in range(N):
1321        if tau[k] >= eps:
1322            if h[k] >= eps:
1323                S = -tau[k]/h[k]
1324
1325                #Update momentum
1326                xmom_update[k] += S*uh[k]
1327                ymom_update[k] += S*vh[k]
1328
[6644]1329def depth_dependent_friction(domain, default_friction,
1330                             surface_roughness_data,
1331                             verbose=False):
1332    """Returns an array of friction values for each wet element adjusted for depth.
1333
1334    Inputs:
1335        domain - computational domain object
1336        default_friction - depth independent bottom friction
1337        surface_roughness_data - N x 5 array of n0, d1, n1, d2, n2 values for each
1338        friction region.
1339
1340    Outputs:
1341        wet_friction - Array that can be used directly to update friction as follows:
1342                       domain.set_quantity('friction', wet_friction)
1343
1344       
1345       
1346    """
1347
[7276]1348    import numpy as num
[6644]1349   
1350    # Create a temp array to store updated depth dependent friction for wet elements
1351    # EHR this is outwardly inneficient but not obvious how to avoid recreating each call??????
1352    N=len(domain)
[7276]1353    wet_friction    = num.zeros(N, num.float)
[6644]1354    wet_friction[:] = default_n0   # Initially assign default_n0 to all array so sure have no zeros values
1355   
1356   
1357    depth = domain.create_quantity_from_expression('stage - elevation')  # create depth instance for this timestep
1358    # Recompute depth as vector 
1359    d = depth.get_values(location='centroids')
1360 
1361    # rebuild the 'friction' values adjusted for depth at this instant
1362    for i in domain.get_wet_elements():                                  # loop for each wet element in domain
1363       
1364        # Get roughness data for each element
1365        n0 = float(surface_roughness_data[i,0])
1366        d1 = float(surface_roughness_data[i,1])
1367        n1 = float(surface_roughness_data[i,2])
1368        d2 = float(surface_roughness_data[i,3])
1369        n2 = float(surface_roughness_data[i,4])
1370       
1371       
1372        # Recompute friction values from depth for this element
1373               
1374        if d[i]   <= d1: 
1375            depth_dependent_friction = n1
1376        elif d[i] >= d2:
1377            depth_dependent_friction = n2
1378        else:
1379            depth_dependent_friction = n1+((n2-n1)/(d2-d1))*(d[i]-d1)
1380           
1381        # check sanity of result
1382        if (depth_dependent_friction  < 0.010 or depth_dependent_friction > 9999.0) :
[7317]1383            log.critical('%s >>>> WARNING: computed depth_dependent friction '
1384                         'out of range, ddf%f, n1=%f, n2=%f'
1385                         % (model_data.basename,
1386                            depth_dependent_friction, n1, n2))
[6644]1387       
1388        # update depth dependent friction  for that wet element
1389        wet_friction[i] = depth_dependent_friction
1390       
1391    # EHR add code to show range of 'friction across domain at this instant as sanity check?????????
1392   
1393    if verbose :
1394        nvals=domain.get_quantity('friction').get_values(location='centroids')        # return array of domain nvals
1395        n_min=min(nvals)
1396        n_max=max(nvals)
1397       
[7317]1398        log.critical('         ++++ calculate_depth_dependent_friction - '
1399                     'Updated friction - range  %7.3f to %7.3f'
1400                     % (n_min, n_max))
[6644]1401   
1402    return wet_friction
1403
1404
1405
[7276]1406################################################################################
[4733]1407# Initialise module
[7276]1408################################################################################
[3804]1409
1410from anuga.utilities import compile
1411if compile.can_use_C_extension('shallow_water_ext.c'):
[7276]1412    # Underlying C implementations can be accessed
[7731]1413    from shallow_water_ext import assign_windfield_values
[4769]1414else:
[7276]1415    msg = 'C implementations could not be accessed by %s.\n ' % __file__
[4769]1416    msg += 'Make sure compile_all.py has been run as described in '
1417    msg += 'the ANUGA installation guide.'
1418    raise Exception, msg
[3804]1419
[4733]1420# Optimisation with psyco
[7810]1421#from anuga.config import use_psyco
1422#if use_psyco:
1423    #try:
1424        #import psyco
1425    #except:
1426        #import os
1427        #if os.name == 'posix' and os.uname()[4] in ['x86_64', 'ia64']:
1428            #pass
1429            ##Psyco isn't supported on 64 bit systems, but it doesn't matter
1430        #else:
1431            #msg = ('WARNING: psyco (speedup) could not be imported, '
1432                   #'you may want to consider installing it')
1433            #log.critical(msg)
1434    #else:
1435        #psyco.bind(Domain.distribute_to_vertices_and_edges)
1436        #psyco.bind(Domain.compute_fluxes)
[3804]1437
[7276]1438
[3804]1439if __name__ == "__main__":
1440    pass
Note: See TracBrowser for help on using the repository browser.