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

Last change on this file since 7870 was 7870, checked in by hudson, 9 years ago

Cleaned up pylint warnings.

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