Changeset 7870


Ignore:
Timestamp:
Jun 23, 2010, 2:33:54 PM (14 years ago)
Author:
hudson
Message:

Cleaned up pylint warnings.

Location:
trunk/anuga_core/source/anuga
Files:
3 edited

Legend:

Unmodified
Added
Removed
  • trunk/anuga_core/source/anuga/file/test_csv.py

    r7866 r7870  
    33import tempfile
    44import numpy as num
     5
     6from anuga.utilities.system_tools import get_pathname_from_package
    57
    68from csv_file import load_csv_as_array, load_csv_as_dict, store_parameters, \
  • trunk/anuga_core/source/anuga/file/test_sww.py

    r7866 r7870  
    99from anuga.shallow_water.shallow_water_domain import Domain
    1010from sww import load_sww_as_domain, weed, get_mesh_and_quantities_from_file
     11from Scientific.IO.NetCDF import NetCDFFile
    1112
    1213# boundary functions
  • trunk/anuga_core/source/anuga/shallow_water/shallow_water_domain.py

    r7810 r7870  
    6565    TRAC administration of ANUGA (User Manuals etc) at
    6666    https://datamining.anu.edu.au/anuga and Subversion repository at
    67     $HeadURL$
     67    $HeadURL: https://datamining.anu.edu.au/svn/anuga/trunk/anuga_core/source/
     68                anuga/shallow_water/shallow_water_domain.py $
    6869
    6970Constraints: See GPL license in the user guide
     
    7475"""
    7576
    76 # Subversion keywords:
    77 #
    78 # $LastChangedDate$
    79 # $LastChangedRevision$
    80 # $LastChangedBy$
    81 
    8277
    8378import numpy as num
     
    8782
    8883from anuga.shallow_water.forcing import Cross_section
    89 from anuga.pmesh.mesh_interface import create_mesh_from_regions
    90 from anuga.utilities.numerical_tools import gradient, mean, ensure_numeric
    91 from anuga.geospatial_data.geospatial_data import ensure_geospatial
     84from anuga.utilities.numerical_tools import mean
    9285from anuga.file.sww import SWW_file 
    93 from anuga.config import netcdf_mode_r, netcdf_mode_w, netcdf_mode_a
    94 
    95 from anuga.fit_interpolate.interpolate import Modeltime_too_late, \
    96                                               Modeltime_too_early
    97 
    98 from anuga.geometry.polygon import inside_polygon, polygon_area, \
    99                                     is_inside_polygon
     86           
    10087import anuga.utilities.log as log
    10188
    10289import types
    103 from types import IntType, FloatType
    104 from warnings import warn
    105 
    106 
    107 ################################################################################
    108 # Shallow water domain
    109 ################################################################################
    110 
    111 ##
    112 # @brief Class for a shallow water domain.
     90
    11391class Domain(Generic_Domain):
    114 
    115     #conserved_quantities = ['stage', 'xmomentum', 'ymomentum']
    116     #other_quantities = ['elevation', 'friction']
    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
    129     # @param evolved_quantities
    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
     92    """ Class for a shallow water domain."""
    13693    def __init__(self,
    13794                 coordinates=None,
     
    153110                 number_of_full_nodes=None,
    154111                 number_of_full_triangles=None):
     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        """
    155131
    156132        # Define quantities for the shallow_water domain
     
    165141       
    166142        Generic_Domain.__init__(self,
    167                                 coordinates,
    168                                 vertices,
    169                                 boundary,
    170                                 conserved_quantities,
    171                                 evolved_quantities,
    172                                 other_quantities,
    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,
    182                                 numproc,
    183                                 number_of_full_nodes=number_of_full_nodes,
    184                                 number_of_full_triangles=number_of_full_triangles)
     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)
    185161
    186162        self.set_defaults()
     
    200176
    201177
    202 
    203 
    204     ##
    205     # @brief Set default values, usually read in from a config file
    206     # @param flag
    207178    def set_defaults(self):
    208179        """Set the default values in this routine. That way we can inherit class
     
    212183        from anuga.config import minimum_storable_height
    213184        from anuga.config import minimum_allowed_height, maximum_allowed_speed
    214         from anuga.config import g, epsilon, beta_w, beta_w_dry,\
     185        from anuga.config import g, beta_w, beta_w_dry, \
    215186             beta_uh, beta_uh_dry, beta_vh, beta_vh_dry, tight_slope_limiters
    216187        from anuga.config import alpha_balance
     
    219190        from anuga.config import use_edge_limiter
    220191        from anuga.config import use_centroid_velocities
    221 
    222 
    223192
    224193        self.set_minimum_allowed_height(minimum_allowed_height)
     
    247216        self.use_centroid_velocities = use_centroid_velocities       
    248217
    249     ##
    250     # @brief
    251     # @param flag
     218
    252219    def set_new_mannings_function(self, flag=True):
    253220        """Cludge to allow unit test to pass, but to
     
    263230
    264231
    265     ##
    266     # @brief
    267     # @param flag
    268232    def set_use_edge_limiter(self, flag=True):
    269233        """Cludge to allow unit test to pass, but to
     
    277241
    278242
    279          
    280     ##
    281     # @brief
    282     # @param beta
    283243    def set_all_limiters(self, beta):
    284244        """Shorthand to assign one constant value [0,1] to all limiters.
     
    299259        self.quantities['ymomentum'].beta = beta
    300260
    301     ##
    302     # @brief
    303     # @param flag
    304     # @param reduction
     261
    305262    def set_store_vertices_uniquely(self, flag, reduction=None):
    306263        """Decide whether vertex values should be stored uniquely as
     
    406363        this function will have no effect.
    407364       
    408         The format, where q is a list of names is for backwards compatibility only.
     365        The format, where q is a list of names is for backwards compatibility
     366        only.
    409367        It will take the specified quantities to be time dependent and assume
    410368        'elevation' to be static regardless.
     
    416374            return
    417375
    418         # Check correcness
     376        # Check correctness
    419377        for quantity_name in q:
    420378            msg = ('Quantity %s is not a valid conserved quantity'
     
    422380            assert quantity_name in self.quantities, msg
    423381
    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            
    441382        assert type(q) == types.DictType   
    442383        self.quantities_to_be_stored = q
     
    567508
    568509
    569     ##
    570     # @brief Run integrity checks on shallow water domain.
    571510    def check_integrity(self):
     511        """ Run integrity checks on shallow water domain. """
    572512        Generic_Domain.check_integrity(self)
    573513
     
    580520        assert self.conserved_quantities[2] == 'ymomentum', msg
    581521
    582     ##
    583     # @brief
    584522    def extrapolate_second_order_sw(self):
    585         #Call correct module function (either from this module or C-extension)
     523        """Call correct module function
     524            (either from this module or C-extension)"""
    586525        extrapolate_second_order_sw(self)
    587526
    588     ##
    589     # @brief
    590527    def compute_fluxes(self):
    591         #Call correct module function (either from this module or C-extension)
     528        """Call correct module function
     529            (either from this module or C-extension)"""
    592530        compute_fluxes(self)
    593531
    594     ##
    595     # @brief
    596532    def distribute_to_vertices_and_edges(self):
    597         # Call correct module function
     533        """ Call correct module function """
    598534        if self.use_edge_limiter:
    599535            distribute_using_edge_limiter(self)
     
    603539
    604540
    605     ##
    606     # @brief Evolve the model by one step.
    607     # @param yieldstep
    608     # @param finaltime
    609     # @param duration
    610     # @param skip_initial_step
    611541    def evolve(self,
    612542               yieldstep=None,
     
    614544               duration=None,
    615545               skip_initial_step=False):
    616         """Specialisation of basic evolve method from parent class"""
     546        """Specialisation of basic evolve method from parent class.
     547       
     548            Evolve the model by 1 step.
     549        """
    617550
    618551        # Call check integrity here rather than from user scripts
     
    642575            yield(t)
    643576
    644     ##
    645     # @brief
     577
    646578    def initialise_storage(self):
    647579        """Create and initialise self.writer object for storing data.
     
    655587        self.writer.store_connectivity()
    656588
    657     ##
    658     # @brief
    659     # @param name
     589
    660590    def store_timestep(self):
    661591        """Store time dependent quantities and time.
     
    667597        self.writer.store_timestep()
    668598
    669     ##
    670     # @brief Get time stepping statistics string for printing.
    671     # @param track_speeds
    672     # @param triangle_id
     599
    673600    def timestepping_statistics(self,
    674601                                track_speeds=False,
     
    712639            Ch = Cw-Cz
    713640
    714             s  = '    %s: vertex_values =  %.4f,\t %.4f,\t %.4f\n'\
    715                  %(name.ljust(qwidth), Vh[0], Vh[1], Vh[2])
    716 
    717             s += '    %s: edge_values =    %.4f,\t %.4f,\t %.4f\n'\
    718                  %(name.ljust(qwidth), Eh[0], Eh[1], Eh[2])
    719 
    720             s += '    %s: centroid_value = %.4f\n'\
    721                  %(name.ljust(qwidth), Ch[0])
    722 
    723             msg += s
     641            message  = '    %s: vertex_values =  %.4f,\t %.4f,\t %.4f\n'\
     642                 % (name.ljust(qwidth), Vh[0], Vh[1], Vh[2])
     643
     644            message += '    %s: edge_values =    %.4f,\t %.4f,\t %.4f\n'\
     645                 % (name.ljust(qwidth), Eh[0], Eh[1], Eh[2])
     646
     647            message += '    %s: centroid_value = %.4f\n'\
     648                 % (name.ljust(qwidth), Ch[0])
     649
     650            msg += message
    724651
    725652            uh = self.quantities['xmomentum']
     
    739666            Cu = Cuh/(Ch + epsilon)
    740667            name = 'U'
    741             s  = '    %s: vertex_values =  %.4f,\t %.4f,\t %.4f\n'\
    742                  %(name.ljust(qwidth), Vu[0], Vu[1], Vu[2])
    743 
    744             s += '    %s: edge_values =    %.4f,\t %.4f,\t %.4f\n'\
    745                  %(name.ljust(qwidth), Eu[0], Eu[1], Eu[2])
    746 
    747             s += '    %s: centroid_value = %.4f\n'\
    748                  %(name.ljust(qwidth), Cu[0])
    749 
    750             msg += s
     668            message  = '    %s: vertex_values =  %.4f,\t %.4f,\t %.4f\n' \
     669                 % (name.ljust(qwidth), Vu[0], Vu[1], Vu[2])
     670
     671            message += '    %s: edge_values =    %.4f,\t %.4f,\t %.4f\n' \
     672                 % (name.ljust(qwidth), Eu[0], Eu[1], Eu[2])
     673
     674            message += '    %s: centroid_value = %.4f\n' \
     675                 % (name.ljust(qwidth), Cu[0])
     676
     677            msg += message
    751678
    752679            Vv = Vvh/(Vh + epsilon)
     
    754681            Cv = Cvh/(Ch + epsilon)
    755682            name = 'V'
    756             s  = '    %s: vertex_values =  %.4f,\t %.4f,\t %.4f\n'\
    757                  %(name.ljust(qwidth), Vv[0], Vv[1], Vv[2])
    758 
    759             s += '    %s: edge_values =    %.4f,\t %.4f,\t %.4f\n'\
    760                  %(name.ljust(qwidth), Ev[0], Ev[1], Ev[2])
    761 
    762             s += '    %s: centroid_value = %.4f\n'\
     683            message  = '    %s: vertex_values =  %.4f,\t %.4f,\t %.4f\n' \
     684                 % (name.ljust(qwidth), Vv[0], Vv[1], Vv[2])
     685
     686            message += '    %s: edge_values =    %.4f,\t %.4f,\t %.4f\n' \
     687                 % (name.ljust(qwidth), Ev[0], Ev[1], Ev[2])
     688
     689            message += '    %s: centroid_value = %.4f\n'\
    763690                 %(name.ljust(qwidth), Cv[0])
    764691
    765             msg += s
     692            msg += message
    766693
    767694            # Froude number in each direction
     
    771698            Cfx = Cu/(num.sqrt(g*Ch) + epsilon)
    772699
    773             s  = '    %s: vertex_values =  %.4f,\t %.4f,\t %.4f\n'\
    774                  %(name.ljust(qwidth), Vfx[0], Vfx[1], Vfx[2])
    775 
    776             s += '    %s: edge_values =    %.4f,\t %.4f,\t %.4f\n'\
    777                  %(name.ljust(qwidth), Efx[0], Efx[1], Efx[2])
    778 
    779             s += '    %s: centroid_value = %.4f\n'\
    780                  %(name.ljust(qwidth), Cfx[0])
    781 
    782             msg += s
     700            message  = '    %s: vertex_values =  %.4f,\t %.4f,\t %.4f\n'\
     701                 % (name.ljust(qwidth), Vfx[0], Vfx[1], Vfx[2])
     702
     703            message += '    %s: edge_values =    %.4f,\t %.4f,\t %.4f\n'\
     704                 % (name.ljust(qwidth), Efx[0], Efx[1], Efx[2])
     705
     706            message += '    %s: centroid_value = %.4f\n'\
     707                 % (name.ljust(qwidth), Cfx[0])
     708
     709            msg += message
    783710
    784711            name = 'Froude (y)'
     
    787714            Cfy = Cv/(num.sqrt(g*Ch) + epsilon)
    788715
    789             s  = '    %s: vertex_values =  %.4f,\t %.4f,\t %.4f\n'\
    790                  %(name.ljust(qwidth), Vfy[0], Vfy[1], Vfy[2])
    791 
    792             s += '    %s: edge_values =    %.4f,\t %.4f,\t %.4f\n'\
    793                  %(name.ljust(qwidth), Efy[0], Efy[1], Efy[2])
    794 
    795             s += '    %s: centroid_value = %.4f\n'\
    796                  %(name.ljust(qwidth), Cfy[0])
    797 
    798             msg += s
     716            message  = '    %s: vertex_values =  %.4f,\t %.4f,\t %.4f\n'\
     717                 % (name.ljust(qwidth), Vfy[0], Vfy[1], Vfy[2])
     718
     719            message += '    %s: edge_values =    %.4f,\t %.4f,\t %.4f\n'\
     720                 % (name.ljust(qwidth), Efy[0], Efy[1], Efy[2])
     721
     722            message += '    %s: centroid_value = %.4f\n'\
     723                 % (name.ljust(qwidth), Cfy[0])
     724
     725            msg += message
    799726
    800727        return msg
     
    813740        # Run through boundary array and compute for each segment
    814741        # the normal momentum ((uh, vh) dot normal) times segment length.
    815         # Based on sign accumulate this into boundary_inflow and boundary_outflow.
     742        # Based on sign accumulate this into boundary_inflow and
     743        # boundary_outflow.
    816744                       
    817745        # Compute flows along boundary
     
    842770            if edge_flow > 0:
    843771                # Flow is inflow     
    844                 total_boundary_inflow += edge_flow                                  
     772                total_boundary_inflow += edge_flow       
    845773            else:
    846774                # Flow is outflow
     
    881809       
    882810        area = self.mesh.get_areas()
    883         volume = 0.0
    884811       
    885812        stage = self.get_quantity('stage').get_values(location='centroids')
    886         elevation = self.get_quantity('elevation').get_values(location='centroids')       
     813        elevation = \
     814            self.get_quantity('elevation').get_values(location='centroids')
    887815        depth = stage-elevation
    888816       
     
    897825         total_boundary_outflow) = self.compute_boundary_flows()
    898826       
    899         s = '---------------------------\n'       
    900         s += 'Volumetric balance report:\n'
    901         s += '--------------------------\n'
    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       
    904         s += 'Net boundary flow by tags [m^3/s]\n'
     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'
    905833        for tag in boundary_flows:
    906             s += '    %s [m^3/s]: %.2f\n' % (tag, boundary_flows[tag])
    907        
    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()
    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.
     834            message += '    %s [m^3/s]: %.2f\n' % (tag, boundary_flows[tag])
     835       
     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()
     840       
     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.
    914845       
    915846        # FIXME(Ole): This part is not yet done.               
    916847       
    917         return s       
     848        return message       
    918849           
    919850################################################################################
     
    951882    from shallow_water_ext import compute_fluxes_ext_central \
    952883                                  as compute_fluxes_ext
    953 
    954     N = len(domain)    # number_of_triangles
    955884
    956885    # Shortcuts
     
    1001930    """Wrapper calling C version of extrapolate_second_order_sw"""
    1002931
    1003     import sys
    1004932    from shallow_water_ext import extrapolate_second_order_sw as extrapol2
    1005 
    1006     N = len(domain) # number_of_triangles
    1007933
    1008934    # Shortcuts
     
    1066992            domain.extrapolate_second_order_sw()
    1067993        else:
    1068             raise 'Unknown order'
     994            raise Exception('Unknown order')
    1069995    else:
    1070996        # Old code:
     
    10771003                Q.extrapolate_second_order_and_limit_by_vertex()
    10781004            else:
    1079                 raise 'Unknown order'
     1005                raise Exception('Unknown order')
    10801006
    10811007    # Take bed elevation into account when water heights are small
     
    11191045            Q.extrapolate_second_order_and_limit_by_edge()
    11201046        else:
    1121             raise 'Unknown order'
     1047            raise Exception('Unknown order')
    11221048
    11231049    balance_deep_and_shallow(domain)
     
    12061132    elevation = domain.quantities['elevation']
    12071133
    1208     h = stage.centroid_values - elevation.centroid_values
    1209     z = elevation.vertex_values
    1210 
    1211     x = domain.get_vertex_coordinates()
    1212     g = domain.g
    1213 
    1214     gravity_c(g, h, z, x, xmom_update, ymom_update)    #, 1.0e-6)
     1134    height = stage.centroid_values - elevation.centroid_values
     1135    elevation = elevation.vertex_values
     1136
     1137    point = domain.get_vertex_coordinates()
     1138
     1139    gravity_c(domain.g, height, elevation, point, xmom_update, ymom_update)
    12151140
    12161141##
     
    12411166    ymom_update = ymom.semi_implicit_update
    12421167
    1243     N = len(domain)
    12441168    eps = domain.minimum_allowed_height
    12451169    g = domain.g
    12461170
    12471171    if domain.use_new_mannings:
    1248         manning_friction_new(g, eps, x, w, uh, vh, z, eta, xmom_update, ymom_update)
     1172        manning_friction_new(g, eps, x, w, uh, vh, z, eta, xmom_update, \
     1173                                ymom_update)
    12491174    else:
    1250         manning_friction_old(g, eps, w, uh, vh, z, eta, xmom_update, ymom_update)
     1175        manning_friction_old(g, eps, w, uh, vh, z, eta, xmom_update, \
     1176                                ymom_update)
    12511177   
    12521178
     
    12781204    ymom_update = ymom.explicit_update
    12791205
    1280     N = len(domain)
    12811206    eps = domain.minimum_allowed_height
    1282     g = domain.g
    1283 
    12841207
    12851208    if domain.use_new_mannings:
    1286         manning_friction_new(g, eps, x, w, uh, vh, z, eta, xmom_update, ymom_update)
     1209        manning_friction_new(domain.g, eps, x, w, uh, vh, z, eta, xmom_update, \
     1210                            ymom_update)
    12871211    else:
    1288         manning_friction_old(g, eps, w, uh, vh, z, eta, xmom_update, ymom_update)
    1289 
    1290 
    1291 
    1292 # FIXME (Ole): This was implemented for use with one of the analytical solutions (Sampson?)
     1212        manning_friction_old(domain.g, eps, w, uh, vh, z, eta, xmom_update, \
     1213                            ymom_update)
     1214
     1215
     1216
     1217# FIXME (Ole): This was implemented for use with one of the analytical solutions
    12931218##
    12941219# @brief Apply linear friction to a surface.
     
    13011226    """
    13021227
    1303     from math import sqrt
    1304 
    13051228    w = domain.quantities['stage'].centroid_values
    13061229    z = domain.quantities['elevation'].centroid_values
     
    13141237    ymom_update = domain.quantities['ymomentum'].semi_implicit_update
    13151238
    1316     N = len(domain) # number_of_triangles
     1239    num_tris = len(domain)
    13171240    eps = domain.minimum_allowed_height
    1318     g = domain.g #Not necessary? Why was this added?
    1319 
    1320     for k in range(N):
     1241
     1242    for k in range(num_tris):
    13211243        if tau[k] >= eps:
    13221244            if h[k] >= eps:
     
    13301252                             surface_roughness_data,
    13311253                             verbose=False):
    1332     """Returns an array of friction values for each wet element adjusted for depth.
     1254    """Returns an array of friction values for each wet element adjusted for
     1255            depth.
    13331256
    13341257    Inputs:
    13351258        domain - computational domain object
    13361259        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.
     1260        surface_roughness_data - N x 5 array of n0, d1, n1, d2, n2 values
     1261        for each friction region.
    13391262
    13401263    Outputs:
    1341         wet_friction - Array that can be used directly to update friction as follows:
     1264        wet_friction - Array that can be used directly to update friction as
     1265                        follows:
    13421266                       domain.set_quantity('friction', wet_friction)
    13431267
     
    13451269       
    13461270    """
    1347 
    1348     import numpy as num
    13491271   
    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)
    1353     wet_friction    = num.zeros(N, num.float)
    1354     wet_friction[:] = default_n0   # Initially assign default_n0 to all array so sure have no zeros values
     1272    default_n0 = 0  # James - this was missing, don't know what it should be
    13551273   
     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
    13561282   
    1357     depth = domain.create_quantity_from_expression('stage - elevation')  # create depth instance for this timestep
     1283    # create depth instance for this timestep   
     1284    depth = domain.create_quantity_from_expression('stage - elevation') 
    13581285    # Recompute depth as vector 
    1359     d = depth.get_values(location='centroids')
     1286    d_vals = depth.get_values(location='centroids')
    13601287 
    13611288    # 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        
     1289    # loop for each wet element in domain
     1290   
     1291    for i in domain.get_wet_elements():       
    13641292        # 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])
     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])
    13701297       
    13711298       
    13721299        # Recompute friction values from depth for this element
    13731300               
    1374         if d[i]   <= d1:
    1375             depth_dependent_friction = n1
    1376         elif d[i] >= d2:
    1377             depth_dependent_friction = n2
     1301        if d_vals[i]   <= d1:
     1302            ddf = n1
     1303        elif d_vals[i] >= d2:
     1304            ddf = n2
    13781305        else:
    1379             depth_dependent_friction = n1+((n2-n1)/(d2-d1))*(d[i]-d1)
     1306            ddf = n1 + ((n2-n1)/(d2-d1))*(d_vals[i]-d1)
    13801307           
    13811308        # check sanity of result
    1382         if (depth_dependent_friction  < 0.010 or depth_dependent_friction > 9999.0) :
    1383             log.critical('%s >>>> WARNING: computed depth_dependent friction '
     1309        if (ddf  < 0.010 or \
     1310                            ddf > 9999.0) :
     1311            log.critical('>>>> WARNING: computed depth_dependent friction '
    13841312                         'out of range, ddf%f, n1=%f, n2=%f'
    1385                          % (model_data.basename,
    1386                             depth_dependent_friction, n1, n2))
     1313                         % (ddf, n1, n2))
    13871314       
    13881315        # 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?????????
     1316        wet_friction[i] = ddf
     1317       
     1318    # EHR add code to show range of 'friction across domain at this instant as
     1319    # sanity check?????????
    13921320   
    13931321    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)
     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)
    13971326       
    13981327        log.critical('         ++++ calculate_depth_dependent_friction - '
     
    14081337################################################################################
    14091338
    1410 from anuga.utilities import compile
    1411 if compile.can_use_C_extension('shallow_water_ext.c'):
    1412     # Underlying C implementations can be accessed
    1413     from shallow_water_ext import assign_windfield_values
    1414 else:
     1339def _raise_compile_exception():
     1340    """ Raise exception if compiler not available. """
    14151341    msg = 'C implementations could not be accessed by %s.\n ' % __file__
    14161342    msg += 'Make sure compile_all.py has been run as described in '
    14171343    msg += 'the ANUGA installation guide.'
    1418     raise Exception, msg
    1419 
    1420 # Optimisation with psyco
    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)
    1437 
     1344    raise Exception(msg)
     1345
     1346from anuga.utilities import compile
     1347if not compile.can_use_C_extension('shallow_water_ext.c'):
     1348    _raise_compile_exception()
    14381349
    14391350if __name__ == "__main__":
Note: See TracChangeset for help on using the changeset viewer.