Changeset 7035 for branches/numpy/anuga


Ignore:
Timestamp:
May 14, 2009, 2:24:54 PM (16 years ago)
Author:
rwilson
Message:

Back-merge from Numeric trunk.

Location:
branches/numpy/anuga
Files:
1 deleted
26 edited

Legend:

Unmodified
Added
Removed

  • branches/numpy/anuga/abstract_2d_finite_volumes/domain.py

    r6904 r7035  
    1818from anuga.abstract_2d_finite_volumes.generic_boundary_conditions\
    1919     import File_boundary
     20from anuga.abstract_2d_finite_volumes.generic_boundary_conditions\
     21     import AWI_boundary
    2022from anuga.abstract_2d_finite_volumes.generic_boundary_conditions\
    2123     import Dirichlet_boundary
    … …  
    311313    def get_normal(self, *args, **kwargs):
    312314        return self.mesh.get_normal(*args, **kwargs)
     315
    313316    def get_triangle_containing_point(self, *args, **kwargs):
    314317        return self.mesh.get_triangle_containing_point(*args, **kwargs)

  • branches/numpy/anuga/abstract_2d_finite_volumes/generic_boundary_conditions.py

    r6904 r7035  
    399399            msg += 'vol_id=%s, edge_id=%s' %(str(vol_id), str(edge_id))
    400400            raise Exception, msg
     401
     402class AWI_boundary(Boundary):
     403    """The AWI_boundary reads values for the conserved
     404    quantities (only STAGE) from an sww NetCDF file, and returns interpolated values
     405    at the midpoints of each associated boundary segment.
     406    Time dependency is interpolated linearly.
     407
     408    Assumes that file contains a time series and possibly
     409    also spatial info. See docstring for File_function in util.py
     410    for details about admissible file formats
     411
     412    AWI boundary must read and interpolate from *smoothed* version
     413    as stored in sww and cannot work with the discontinuos triangles.
     414
     415    Example:
     416    Ba = AWI_boundary('source_file.sww', domain)
     417
     418
     419    Note that the resulting solution history is not exactly the same as if
     420    the models were coupled as there is no feedback into the source model.
     421
     422    This was added by Nils Goseberg et al in April 2009
     423    """
     424
     425    def __init__(self, filename, domain, time_thinning=1,
     426                 use_cache=False, verbose=False):
     427        import time
     428        from anuga.config import time_format
     429        from anuga.abstract_2d_finite_volumes.util import file_function
     430
     431        Boundary.__init__(self)
     432
     433        # Get x,y vertex coordinates for all triangles
     434        V = domain.vertex_coordinates
     435
     436        # Compute midpoint coordinates for all boundary elements
     437        # Only a subset may be invoked when boundary is evaluated but
     438        # we don't know which ones at this stage since this object can
     439        # be attached to
     440        # any tagged boundary later on.
     441
     442        if verbose: print 'Find midpoint coordinates of entire boundary'
     443        self.midpoint_coordinates = num.zeros((len(domain.boundary), 2), num.float)
     444        boundary_keys = domain.boundary.keys()
     445
     446        xllcorner = domain.geo_reference.get_xllcorner()
     447        yllcorner = domain.geo_reference.get_yllcorner()       
     448
     449        # Make ordering unique #FIXME: should this happen in domain.py?
     450        boundary_keys.sort()
     451
     452        # Record ordering #FIXME: should this also happen in domain.py?
     453        self.boundary_indices = {}
     454        for i, (vol_id, edge_id) in enumerate(boundary_keys):
     455            base_index = 3*vol_id
     456            x0, y0 = V[base_index, :]
     457            x1, y1 = V[base_index+1, :]
     458            x2, y2 = V[base_index+2, :]
     459           
     460            # Compute midpoints
     461            if edge_id == 0: m = num.array([(x1 + x2)/2, (y1 + y2)/2])
     462            if edge_id == 1: m = num.array([(x0 + x2)/2, (y0 + y2)/2])
     463            if edge_id == 2: m = num.array([(x1 + x0)/2, (y1 + y0)/2])
     464
     465            # Convert to absolute UTM coordinates
     466            m[0] += xllcorner
     467            m[1] += yllcorner
     468           
     469            # Register point and index
     470            self.midpoint_coordinates[i,:] = m
     471
     472            # Register index of this boundary edge for use with evaluate
     473            self.boundary_indices[(vol_id, edge_id)] = i
     474
     475
     476        if verbose: print 'Initialise file_function'
     477        self.F = file_function(filename, domain,
     478                                   interpolation_points=self.midpoint_coordinates,
     479                               time_thinning=time_thinning,
     480                               use_cache=use_cache,
     481                               verbose=verbose)
     482        self.domain = domain
     483
     484        # Test
     485
     486        # Here we'll flag indices outside the mesh as a warning
     487        # as suggested by Joaquim Luis in sourceforge posting
     488        # November 2007
     489        # We won't make it an error as it is conceivable that
     490        # only part of mesh boundary is actually used with a given
     491        # file boundary sww file.
     492        if hasattr(self.F, 'indices_outside_mesh') and
     493               len(self.F.indices_outside_mesh) > 0:
     494            msg = 'WARNING: File_boundary has points outside the mesh '
     495            msg += 'given in %s. ' % filename
     496            msg += 'See warning message issued by Interpolation_function '
     497            msg += 'for details (should appear above somewhere if '
     498            msg += 'verbose is True).\n'
     499            msg += 'This is perfectly OK as long as the points that are '
     500            msg += 'outside aren\'t used on the actual boundary segment.'
     501            if verbose is True:           
     502                print msg
     503            #raise Exception(msg)
     504
     505        q = self.F(0, point_id=0)
     506
     507        d = len(domain.conserved_quantities)
     508        msg = 'Values specified in file %s must be ' % filename
     509        msg += 'a list or an array of length %d' % d
     510        assert len(q) == d, msg
     511
     512
     513    def __repr__(self):
     514        return 'File boundary'
     515
     516
     517    def evaluate(self, vol_id=None, edge_id=None):
     518        """Return linearly interpolated values based on domain.time
     519               at midpoint of segment defined by vol_id and edge_id.
     520        """
     521
     522        q = self.domain.get_conserved_quantities(vol_id, edge=edge_id)
     523        t = self.domain.time
     524
     525        if vol_id is not None and edge_id is not None:
     526            i = self.boundary_indices[vol_id, edge_id]
     527            res = self.F(t, point_id=i)
     528
     529            if res == NAN:
     530                x,y = self.midpoint_coordinates[i,:]
     531                msg = 'NAN value found in file_boundary at '
     532                msg += 'point id #%d: (%.2f, %.2f).\n' % (i, x, y)
     533
     534                if (hasattr(self.F, 'indices_outside_mesh') and
     535                       len(self.F.indices_outside_mesh) > 0):
     536                    # Check if NAN point is due it being outside
     537                    # boundary defined in sww file.
     538
     539                    if i in self.F.indices_outside_mesh:
     540                        msg += 'This point refers to one outside the '
     541                        msg += 'mesh defined by the file %s.\n' % self.F.filename
     542                        msg += 'Make sure that the file covers '
     543                        msg += 'the boundary segment it is assigned to '
     544                        msg += 'in set_boundary.'
     545                    else:
     546                        msg += 'This point is inside the mesh defined '
     547                        msg += 'the file %s.\n' % self.F.filename
     548                        msg += 'Check this file for NANs.'
     549                raise Exception, msg
     550           
     551            q[0] = res[0] # Take stage, leave momentum alone
     552            return q
     553            #return res
     554        else:
     555            # raise 'Boundary call without point_id not implemented'
     556            # FIXME: What should the semantics be?
     557            return self.F(t)
     558
     559

  • branches/numpy/anuga/abstract_2d_finite_volumes/quantity.py

    r6902 r7035  
    459459            raise Exception, msg
    460460
    461         msg = 'Indices must be a list or None'
    462         assert (indices is None
    463                 or isinstance(indices, (list, num.ndarray))), msg
     461        msg = 'Indices must be a list, array or None'
     462        #assert (indices is None
     463        assert isinstance(indices, (None, list, num.ndarray))), msg
    464464
    465465        # Determine which 'set_values_from_...' to use
    … …  
    11871187        import types
    11881188
    1189         msg = 'Indices must be a list or None'
    1190         assert indices is None or isinstance(indices, (list, num.ndarray)), msg
     1189        msg = 'Indices must be a list, array or None'
     1190        #assert indices is None or isinstance(indices, (list, num.ndarray)), msg
     1191        assert isinstance(indices, (list, None, num.ndarray)), msg
    11911192
    11921193        if location == 'centroids':

  • branches/numpy/anuga/abstract_2d_finite_volumes/test_domain.py

    r6902 r7035  
    614614
    615615        sem = num.array([1.,1.,1.,1.])/num.array([1, 2, 3, 4])
    616         denom = num.ones(4, num.float)-domain.timestep*sem
     616        denom = num.ones(4, num.float) - domain.timestep*sem
    617617
    618618#        x = array([1, 2, 3, 4]) + array( [.4,.3,.2,.1] )

  • branches/numpy/anuga/coordinate_transforms/geo_reference.py

    r6689 r7035  
    319319        """
    320320
    321         is_list = False
    322         if type(points) == types.ListType:
    323             is_list = True
     321        # remember if we got a list
     322        is_list = isinstance(points, list)
    324323
    325324        points = ensure_numeric(points, num.float)

  • branches/numpy/anuga/culvert_flows/culvert_class.py

    r6902 r7035  
    1515
    1616import numpy as num
    17 
     17from math import sqrt
    1818
    1919class Below_interval(Exception): pass
    … …  
    137137                 manning=None,          # Mannings Roughness for Culvert
    138138                 sum_loss=None,
    139                  use_velocity_head=False,
     139                 use_velocity_head=False, # FIXME(Ole): Get rid of - always True
    140140                 use_momentum_jet=False, # FIXME(Ole): Not yet implemented
    141141                 label=None,
    … …  
    464464            log_filename = self.log_filename
    465465           
    466         # Compute stage and energy at the
     466        # Compute stage, energy and velocity at the
    467467        # enquiry points at each end of the culvert
    468468        openings = self.openings
    … …  
    474474            depth = h = stage-opening.elevation
    475475                                                           
    476                                                
     476           
     477            # Get velocity                                 
     478            xmomentum = dq['xmomentum'].get_values(location='centroids',
     479                                                   indices=[idx])[0]
     480            ymomentum = dq['xmomentum'].get_values(location='centroids',
     481                                                   indices=[idx])[0]
     482           
     483            if h > minimum_allowed_height:
     484                u = xmomentum/(h + velocity_protection/h)
     485                v = ymomentum/(h + velocity_protection/h)
     486            else:
     487                u = v = 0.0
     488               
     489            v_squared = u*u + v*v
     490           
    477491            if self.use_velocity_head is True:
    478                 xmomentum = dq['xmomentum'].get_values(location='centroids',
    479                                                        indices=[idx])[0]
    480                 ymomentum = dq['xmomentum'].get_values(location='centroids',
    481                                                        indices=[idx])[0]
    482            
    483                 if h > minimum_allowed_height:
    484                     u = xmomentum/(h + velocity_protection/h)
    485                     v = ymomentum/(h + velocity_protection/h)
    486                 else:
    487                     u = v = 0.0
    488                
    489                 velocity_head = 0.5*(u*u + v*v)/g   
     492                velocity_head = 0.5*v_squared/g   
    490493            else:
    491494                velocity_head = 0.0
    … …  
    495498            opening.stage = stage
    496499            opening.depth = depth
     500            opening.velocity = sqrt(v_squared)
    497501           
    498502
    … …  
    583587                    self.culvert_routine(inlet.depth,
    584588                                         outlet.depth,
     589                                         inlet.velocity,
     590                                         outlet.velocity,
    585591                                         inlet.specific_energy,
    586592                                         delta_total_energy,

  • branches/numpy/anuga/culvert_flows/culvert_routines.py

    r6902 r7035  
    1919def boyd_generalised_culvert_model(inlet_depth,
    2020                                   outlet_depth,
     21                                   inlet_velocity,
     22                                   outlet_velocity,
    2123                                   inlet_specific_energy,
    2224                                   delta_total_energy,

  • branches/numpy/anuga/culvert_flows/test_culvert_routines.py

    r6902 r7035  
    3131        inlet_depth=2.0
    3232        outlet_depth=0.0
     33       
     34        inlet_velocity=0.0,
     35        outlet_velocity=0.0,       
    3336
    3437        culvert_length=4.0
    … …  
    4952        Q, v, d = boyd_generalised_culvert_model(inlet_depth,
    5053                                                 outlet_depth,
     54                                                 inlet_velocity,
     55                                                 outlet_velocity,
    5156                                                 inlet_specific_energy,
    5257                                                 delta_total_energy,
    … …  
    8085        outlet_depth=0.0
    8186
     87        inlet_velocity=0.0,
     88        outlet_velocity=0.0,               
     89       
    8290        culvert_length=4.0
    8391        culvert_width=1.2
    … …  
    97105        Q, v, d = boyd_generalised_culvert_model(inlet_depth,
    98106                                                 outlet_depth,
     107                                                 inlet_velocity,
     108                                                 outlet_velocity,
    99109                                                 inlet_specific_energy,
    100110                                                 delta_total_energy,
    … …  
    187197        Q, v, d = boyd_generalised_culvert_model(inlet_depth,
    188198                                                 outlet_depth,
     199                                                 inlet_velocity,
     200                                                 outlet_velocity,
    189201                                                 inlet_specific_energy,
    190202                                                 delta_total_energy,

  • branches/numpy/anuga/fit_interpolate/interpolate.py

    r6689 r7035  
    705705                              imposing a particular order on the output vector.
    706706        vertex_coordinates:   mx2 array of coordinates (float)
    707         triangles:            nx3 array of indices into vertex_coordinates (Int)
     707        triangles:            nx3 array of indices into vertex_coordinates (int)
    708708        interpolation_points: Nx2 array of coordinates to be interpolated to
    709709        verbose:              Level of reporting

  • branches/numpy/anuga/geospatial_data/test_geospatial_data.py

    r7009 r7035  
    15301530                        'Bad zone error!')
    15311531
     1532        # use self.failUnlessRaises(ValueError, Geospatial_data(latitudes=lats))
    15321533        try:
    15331534            results = Geospatial_data(latitudes=lats)

  • branches/numpy/anuga/load_mesh/loadASCII.py

    r6904 r7035  
    11011101def take_points(dict, indices_to_keep):
    11021102    dict = point_atts2array(dict)
     1103    # FIXME maybe the points data structure should become a class?
    11031104    dict['pointlist'] = num.take(dict['pointlist'], indices_to_keep, axis=0)
    11041105

  • branches/numpy/anuga/load_mesh/test_loadASCII.py

    r6904 r7035  
    260260    def test_load_bad_no_file_tsh(self):
    261261        fileName = tempfile.mktemp('.tsh')
     262        # use self.faileUnlessRaises(IOError, import_mesh_file(fileName))
    262263        try:
    263264            dict = import_mesh_file(fileName)

  • branches/numpy/anuga/mesh_engine/mesh_engine.py

    r6304 r7035  
    6060       
    6161    try:
    62         # If int is used, instead of int32, it fails in Linux
     62        # If num.int is used, instead of num.int32, it fails in Linux
    6363        segments = ensure_numeric(segments, num.int32)
    6464       

  • branches/numpy/anuga/pmesh/mesh.py

    r6304 r7035  
    134134        self.attributes = []
    135135        if not attributes is None:
    136             self.attributes=attributes
     136            self.attributes = attributes
    137137   
    138138    def setAttributes(self,attributes):

  • branches/numpy/anuga/pmesh/test_mesh_interface.py

    r6360 r7035  
    589589        inner_polygon = [[800,400], [900,500], [800,600]]
    590590        interior_regions = [(inner_polygon, 5)]
     591
     592        m = create_mesh_from_regions(polygon,
     593                                     boundary_tags,
     594                                     10000000,
     595                                     interior_regions=interior_regions)
    591596        try:
    592597            m = create_mesh_from_regions(polygon,
    … …  
    600605            msg += 'does not cause an Exception to be raised'
    601606            raise Exception, msg
     607
     608           
     609    def test_create_mesh_with_segments_out_of_bounds(self):
     610        """Test that create_mesh_from_regions fails when a segment is out of bounds.
     611        """
     612       
     613        # These are the absolute values
     614        min_x = 10
     615        min_y = 88
     616        polygon = [[min_x,min_y],[1000,100],[1000,1000],[100,1000]]
     617
     618       
     619        boundary_tags = {'walls': [0,1], 'bom':[2,3], 'out': [5]}
     620
     621        # This one is inside bounding polygon - should pass
     622        inner_polygon = [[800,400],[900,500],[800,600]]
     623       
     624        interior_regions = [(inner_polygon, 5)]
     625
     626
     627        try:
     628            m = create_mesh_from_regions(polygon,
     629                                         boundary_tags,
     630                                         10000000,
     631                                         interior_regions=interior_regions)
     632        except:
     633            pass
     634        else:
     635            msg = 'Tags are listed repeatedly, but create mesh from regions '
     636            msg += 'does not cause an Exception to be raised'
     637            raise Exception, msg
     638           
     639                   
     640
    602641
    603642    def test_create_mesh_from_regions_with_duplicate_verts(self):

  • branches/numpy/anuga/shallow_water/__init__.py

    r6553 r7035  
    1313     Dirichlet_discharge_boundary,\
    1414     Field_boundary,\
    15      Transmissive_stage_zero_momentum_boundary
     15     Transmissive_stage_zero_momentum_boundary,\
     16     AWI_boundary
    1617
    1718

  • branches/numpy/anuga/shallow_water/data_manager.py

    r7009 r7035  
    27822782# @param verbose True if this function is to be verbose.
    27832783def _convert_dem_from_ascii2netcdf(basename_in, basename_out = None,
    2784                                   verbose = False):
     2784                                   verbose = False):
    27852785    """Read Digital Elevation model from the following ASCII format (.asc)
    27862786
    … …  
    29202920        if verbose and i % ((n+10)/10) == 0:
    29212921            print 'Processing row %d of %d' % (i, nrows)
     2922           
     2923        if len(fields) != ncols:
     2924            msg = 'Wrong number of columns in file "%s" line %d\n' % (basename_in + '.asc', i)
     2925            msg += 'I got %d elements, but there should have been %d\n' % (len(fields), ncols)
     2926            raise Exception, msg
     2927
    29222928        elevation[i, :] = num.array([float(x) for x in fields])
    29232929
    … …  
    53895395
    53905396    if weights is None:
    5391         weights = num.ones(numSrc, num.int)     #array default#
     5397        weights = num.ones(numSrc)
    53925398
    53935399    if permutation is None:

  • branches/numpy/anuga/shallow_water/shallow_water_domain.py

    r6689 r7035  
    9595from anuga.abstract_2d_finite_volumes.generic_boundary_conditions\
    9696     import Transmissive_boundary
     97from anuga.abstract_2d_finite_volumes.generic_boundary_conditions\
     98     import AWI_boundary
     99
    97100from anuga.pmesh.mesh_interface import create_mesh_from_regions
    98101from anuga.utilities.numerical_tools import gradient, mean, ensure_numeric
    … …  
    271274    # @param minimum_allowed_height
    272275    def set_minimum_allowed_height(self, minimum_allowed_height):
    273         """Set the minimum depth recognised in the numerical scheme.
     276        """Set minimum depth that will be recognised in the numerical scheme.
    274277
    275278        The minimum allowed depth is in meters.
    … …  
    19521955                            s_vec, phi_vec, const):
    19531956    """Python version of assigning wind field to update vectors.
    1954     A c version also exists (for speed)
     1957    A C version also exists (for speed)
    19551958    """
    19561959
    … …  
    20572060
    20582061        # Update area if applicable
    2059         self.exchange_area = None
    20602062        if center is not None and radius is not None:
    20612063            assert len(center) == 2
    20622064            msg = 'Polygon cannot be specified when center and radius are'
    20632065            assert polygon is None, msg
    2064 
    2065             self.exchange_area = radius**2*pi
    20662066
    20672067            # Check that circle center lies within the mesh.

  • branches/numpy/anuga/shallow_water/test_all.py

    r6410 r7035  
    7979
    8080if __name__ == '__main__':   
    81 
    8281    suite = regressionTest()
    8382    runner = unittest.TextTestRunner() #verbosity=2)

  • branches/numpy/anuga/shallow_water/test_data_manager.py

    r7009 r7035  
    49604960
    49614961        ## 7th test
    4962         m2d = num.array([[0,1,2,3],[4,5,6,7],[8,9,10,11],[12,13,14,15]])
     4962        m2d = num.array([[0,1,2,3],[4,5,6,7],[8,9,10,11],[12,13,14,15]], num.int)    #array default#
    49634963        kmin, kmax, lmin, lmax = data_manager._get_min_max_indexes(
    49644964            latitudes,longitudes,
    … …  
    50195019        longitudes = [148,149,150,151]
    50205020
    5021         m2d = num.array([[0,1,2,3],[4,5,6,7],[8,9,10,11],[12,13,14,15]])
     5021        m2d = num.array([[0,1,2,3],[4,5,6,7],[8,9,10,11],[12,13,14,15]], num.int)    #array default#
    50225022
    50235023        # k - lat
    … …  
    60436043                quantities_init[2].append(-va[i]) # South is negative in MUX
    60446044
    6045         file_handle, base_name = tempfile.mkstemp("")
     6045        file_handle, base_name = tempfile.mkstemp("write_mux2")
    60466046        os.close(file_handle)
    60476047        os.remove(base_name)

  • branches/numpy/anuga/shallow_water/test_shallow_water_domain.py

    r6902 r7035  
    33023302        domain.tight_slope_limiters = 0
    33033303        domain.distribute_to_vertices_and_edges()
    3304         assert num.allclose(L[1], [4.1, 16.1, 20.1])
     3304        assert num.allclose(L[1], [4.1, 16.1, 20.1])       
    33053305        for i in range(len(L)):
    33063306            assert num.allclose(volumes[i], num.sum(L[i])/3)
    3307 
    3308         # Allow triangle to be flatter (closer to bed)
    3309         domain.tight_slope_limiters = 1
    3310 
     3307       
     3308               
     3309        domain.tight_slope_limiters = 1 # Allow triangle to be flatter (closer to bed)
    33113310        domain.distribute_to_vertices_and_edges()
    33123311        assert num.allclose(L[1], [4.2386, 16.0604, 20.001])
    33133312        for i in range(len(L)):
    3314             assert num.allclose(volumes[i], num.sum(L[i])/3)
     3313            assert num.allclose(volumes[i], num.sum(L[i])/3)   
     3314       
    33153315
    33163316        domain._order_ = 2
    3317 
    3318         domain.tight_slope_limiters = 0
     3317       
     3318        domain.tight_slope_limiters = 0   
    33193319        domain.distribute_to_vertices_and_edges()
    33203320        assert num.allclose(L[1], [4.1, 16.1, 20.1])
    … …  
    33263326
    33273327        domain.distribute_to_vertices_and_edges()
    3328         assert (num.allclose(L[1], [4.23370103, 16.06529897, 20.001]) or
    3329                 num.allclose(L[1], [4.18944138, 16.10955862, 20.001]) or
    3330                 num.allclose(L[1], [4.19351461, 16.10548539, 20.001]))
    3331         # old limiters
    3332 
     3328        assert num.allclose(L[1], [4.23370103, 16.06529897, 20.001]) or\
     3329               num.allclose(L[1], [4.18944138, 16.10955862, 20.001]) or\
     3330               num.allclose(L[1], [4.19351461, 16.10548539, 20.001]) # old limiters
     3331       
    33333332        for i in range(len(L)):
    33343333            assert num.allclose(volumes[i], num.sum(L[i])/3)
     3334
    33353335
    33363336    def test_second_order_distribute_real_data(self):
    … …  
    61176117
    61186118         msg = "test_tags_to_boundaries failed. Single boundary wasn't added."
    6119          self.failUnless( domain.boundary[(1, 0)]  == '1', msg)
    6120          self.failUnless( domain.boundary[(1, 2)]  == '2', msg)
    6121          self.failUnless( domain.boundary[(0, 1)]  == '3', msg)
    6122          self.failUnless( domain.boundary[(0, 0)]  == 'exterior', msg)
     6119         self.failUnless(domain.boundary[(1, 0)]  == '1', msg)
     6120         self.failUnless(domain.boundary[(1, 2)]  == '2', msg)
     6121         self.failUnless(domain.boundary[(0, 1)]  == '3', msg)
     6122         self.failUnless(domain.boundary[(0, 0)]  == 'exterior', msg)
    61236123         msg = "test_pmesh2Domain Too many boundaries"
    6124          self.failUnless( len(domain.boundary)  == 4, msg)
     6124         self.failUnless(len(domain.boundary)  == 4, msg)
    61256125
    61266126         # FIXME change to use get_xllcorner
    … …  
    61796179                           [ 2.66666667, 0.66666667],
    61806180                           [ 0.66666667, 2.66666667],
    6181                            [ 0.0, 1.0],
    6182                            [ 0.0, 3.0],
    6183                            [ 1.0, 0.0],
    6184                            [ 1.0, 1.0],
    6185                            [ 1.0, 2.0],
    6186                            [ 1.0, 3.0],
    6187                            [ 2.0, 1.0],
    6188                            [ 3.0, 0.0],
    6189                            [ 3.0, 1.0]]
     6181                           [ 0.0,        1.0],
     6182                           [ 0.0,        3.0],
     6183                           [ 1.0,        0.0],
     6184                           [ 1.0,        1.0],
     6185                           [ 1.0,        2.0],
     6186                           [ 1.0,        3.0],
     6187                           [ 2.0,        1.0],
     6188                           [ 3.0,        0.0],
     6189                           [ 3.0,        1.0]]
    61906190
    61916191        data_geo_spatial = Geospatial_data(data_points_rel,
    … …  
    62016201        file = open(ptsfile, "w")
    62026202        file.write(" x,y," + att + " \n")
    6203         for data_point,attribute in map(None, data_points_absolute, attributes):
     6203        for data_point, attribute in map(None, data_points_absolute, attributes):
    62046204            row = (str(data_point[0]) + ',' +
    62056205                   str(data_point[1]) + ',' +
    … …  
    65486548
    65496549
    6550     def Xtest_inflow_using_flowline(self):
     6550    def test_inflow_using_flowline(self):
    65516551        """test_inflow_using_flowline
    65526552
    65536553        Test the ability of a flowline to match inflow above the flowline by
    6554         creating constant inflow onto a circle at the head of a 20m wide by 300m
    6555         long plane dipping at various slopes with a perpendicular flowline and
    6556         gauge downstream of the inflow and a 45 degree flowlines at 200m
    6557         downstream.
     6554        creating constant inflow onto a circle at the head of a 20m
     6555        wide by 300m long plane dipping at various slopes with a
     6556        perpendicular flowline and gauge downstream of the inflow and
     6557        a 45 degree flowlines at 200m downstream.
     6558
     6559        A more substantial version of this test with finer resolution and
     6560        including the depth calculation using Manning's equation is
     6561        available under the validate_all suite in the directory
     6562        anuga_validation/automated_validation_tests/flow_tests.
    65586563        """
    65596564
    6560         # FIXME(Ole): Move this and the following test to validate_all.py as
    6561         # they are rather time consuming
    65626565       
    6563         verbose = True
     6566        verbose = False
     6567       
    65646568
    65656569        #----------------------------------------------------------------------
    … …  
    65816585        # Setup computational domain
    65826586        #----------------------------------------------------------------------
    6583 
    6584         number_of_inflows = 2       # Number of inflows on top of each other
    6585         finaltime = 1000.0          # 6000.0
    6586 
    6587         length = 300.
     6587        number_of_inflows = 2 # Number of inflows on top of each other
     6588        finaltime = 500 #700.0 # If this is too short, steady state will not be achieved
     6589
     6590        length = 250.
    65886591        width  = 20.
    65896592        dx = dy = 5                 # Resolution: of grid on both axes
    … …  
    65946597                                                       len2=width)
    65956598
    6596         for mannings_n in [0.150, 0.07, 0.035]:    #[0.012, 0.035, 0.070, 0.150]:
    6597             for slope in [1.0/300, 1.0/150, 1.0/75]:
    6598                 # Loop over a range of bedslopes representing sub to
    6599                 # super critical flows
    6600                 if verbose:
    6601                     print
    6602                     print 'Slope:', slope, 'Mannings n:', mannings_n
     6599        for mannings_n in [0.1, 0.01]:
     6600            # Loop over a range of roughnesses             
     6601           
     6602            for slope in [1.0/300, 1.0/100]:
     6603                # Loop over a range of bedslopes representing
     6604                # sub to super critical flows
     6605
     6606
    66036607                domain = Domain(points, vertices, boundary)   
    6604                 domain.set_name('Inflow_flowline_test')     # Output name
     6608                domain.set_name('inflow_flowline_test')     # Output name
    66056609
    66066610                #--------------------------------------------------------------
    … …  
    66206624
    66216625                #--------------------------------------------------------------
    6622                 # Setup boundary conditions
    6623                 #--------------------------------------------------------------
    6624 
    6625                 Br = Reflective_boundary(domain)      # Solid reflective wall
    6626                 # Outflow stage at -0.9m d=0.1m
    6627                 Bo = Dirichlet_boundary([-100, 0, 0])
    6628 
    6629                 domain.set_boundary({'left': Br, 'right': Bo,
    6630                                      'top': Br,  'bottom': Br})
    6631 
    6632                 #--------------------------------------------------------------
    66336626                # Setup Inflow
    66346627                #--------------------------------------------------------------
    66356628
    66366629                # Fixed Flowrate onto Area
    6637                 fixed_inflow = Inflow(domain, center=(10.0, 10.0),
    6638                                       radius=5.00, rate=10.00)   
     6630                fixed_inflow = Inflow(domain,
     6631                                      center=(10.0, 10.0),
     6632                                      radius=5.00,
     6633                                      rate=10.00)   
    66396634
    66406635                # Stack this flow
    … …  
    66446639                ref_flow = fixed_inflow.rate*number_of_inflows
    66456640
     6641                # Compute normal depth on plane using Mannings equation
     6642                # v=1/n*(r^2/3)*(s^0.5) or r=(Q*n/(s^0.5*W))^0.6
     6643                normal_depth=(ref_flow*mannings_n/(slope**0.5*width))**0.6
     6644                if verbose:
     6645                    print
     6646                    print 'Slope:', slope, 'Mannings n:', mannings_n
     6647                   
     6648                   
     6649                #--------------------------------------------------------------
     6650                # Setup boundary conditions
     6651                #--------------------------------------------------------------
     6652
     6653                Br = Reflective_boundary(domain)
     6654               
     6655                # Define downstream boundary based on predicted depth
     6656                def normal_depth_stage_downstream(t):
     6657                    return (-slope*length) + normal_depth
     6658               
     6659                Bt = Transmissive_momentum_set_stage_boundary(domain=domain,
     6660                                                              function=normal_depth_stage_downstream)
     6661               
     6662
     6663               
     6664
     6665                domain.set_boundary({'left': Br,
     6666                                     'right': Bt,
     6667                                     'top': Br,
     6668                                     'bottom': Br})
     6669
     6670
     6671
    66466672                #--------------------------------------------------------------
    66476673                # Evolve system through time
    … …  
    66536679                    #    print domain.timestepping_statistics()
    66546680
    6655                 if verbose:
    6656                     print domain.volumetric_balance_statistics()                                                           
     6681                    #    print domain.volumetric_balance_statistics()                                   
     6682
     6683
    66576684                #--------------------------------------------------------------
    66586685                # Compute flow thru flowlines ds of inflow
    … …  
    66826709                assert num.allclose(q, ref_flow, rtol=1.0e-2), msg         
    66836710
    6684                 # Stage recorder (gauge) in middle of plane at 200m
    6685                 x = 200.0
    6686                 y = 10.00
    6687                 w = domain.get_quantity('stage').\
    6688                         get_values(interpolation_points=[[x, y]])[0]
    6689                 z = domain.get_quantity('elevation').\
    6690                         get_values(interpolation_points=[[x, y]])[0]
    6691                 domain_depth = w-z
    6692 
    6693                 # Compute normal depth at gauge location using Manning equation
    6694                 # v=1/n*(r^2/3)*(s^0.5) or r=(Q*n/(s^0.5*W))^0.6
    6695                 normal_depth = (ref_flow*mannings_n/(slope**0.5*width))**0.6
    6696                 msg = ('Predicted depth of flow was %f, should have been %f'
    6697                        % (domain_depth, normal_depth))
    6698                 if verbose:
    6699                     diff = abs(domain_depth-normal_depth)
    6700                     print ('Depth: ANUGA = %f, Mannings = %f, E=%f'
    6701                            % (domain_depth, normal_depth, diff/domain_depth))
    6702                     assert diff < 0.1
    6703                    
    6704                 if slope >= 1.0/100:
    6705                     # Really super critical flow is not as stable.
    6706                     #assert num.allclose(domain_depth, normal_depth,
    6707                     #                    rtol=1.0e-1), msg
    6708                     pass
    6709                 else:
    6710                     pass
    6711                     #assert num.allclose(domain_depth, normal_depth,
    6712                     #                    rtol=1.0e-2), msg
    6713 
    6714 
     6711       
     6712       
    67156713    def Xtest_friction_dependent_flow_using_flowline(self):
    67166714        """test_friction_dependent_flow_using_flowline

  • branches/numpy/anuga/utilities/compile.py

    r6304 r7035  
    264264            % (compiler, FN, I_dirs, python_include, utilities_include_dir, root)
    265265      elif FN == "polygon_ext.c":
    266         # gcc 4.3.x screws up in this file if '-O3' is used
     266        # gcc 4.3.x is problematic with this file if '-O3' is used
    267267        s = '%s -c %s %s -I"%s" -I"%s" -o "%s.o" -Wall'\
    268268            %(compiler, FN, I_dirs, python_include, utilities_include_dir, root)

  • branches/numpy/anuga/utilities/log.py

    r6902 r7035  
    2424'''
    2525
    26 import sys
    2726import os
    2827import sys
    … …  
    6059NOTSET = logging.NOTSET
    6160
    62 # get True if python version 2.5 or later
     61# set new_python to True if python version 2.5 or later
    6362(version_major, version_minor, _, _, _) = sys.version_info
    6463new_python = ((version_major == 2 and version_minor >= 5) or version_major > 2)
    … …  
    103102        console.setFormatter(formatter)
    104103        logging.getLogger('').addHandler(console)
     104
     105        # catch exceptions
     106        sys.excepthook = log_exception_hook
    105107
    106108        # tell the world how we are set up
    … …  
    122124    frames = traceback.extract_stack()
    123125    frames.reverse()
    124     (_, mod_name) = __name__.rsplit('.', 1)
     126    try:
     127        (_, mod_name) = __name__.rsplit('.', 1)
     128    except ValueError:
     129        mod_name = __name__
    125130    for (fpath, lnum, mname, _) in frames:
    126131        (fname, _) = os.path.basename(fpath).rsplit('.', 1)
    … …  
    133138        logging.log(level, msg)
    134139
     140##
     141# @brief Hook function to process uncaught exceptions.
     142# @param type
     143# @param value
     144# @param traceback
     145# @note Same interface as sys.excepthook()
     146def log_exception_hook(type, value, tb):
     147    msg = ''.join(traceback.format_exception(type, value, tb))
     148    critical(msg)
     149
     150   
    135151################################################################################
    136152# Shortcut routines to make for simpler user code.
    … …  
    176192    if sys.platform != 'win32':
    177193        _proc_status = '/proc/%d/status' % os.getpid()
    178         _scale = {'KB': 1024.0, 'MB': 1024.0*1024.0, 'GB': 1024.0*1024.0*1024.0,
    179                   'kB': 1024.0, 'mB': 1024.0*1024.0, 'gB': 1024.0*1024.0*1024.0}
     194        _scale = {'KB': 1024, 'MB': 1024*1024, 'GB': 1024*1024*1024,
     195                  'kB': 1024, 'mB': 1024*1024, 'gB': 1024*1024*1024}
    180196
    181197        def _VmB(VmKey):
    … …  
    214230            return _VmB('VmStk:') - since
    215231
    216         msg = ('Resource usage: memory=%.1f resident=%.1f stacksize=%.1f'
    217                % (memory()/_scale['GB'], resident()/_scale['GB'],
    218                   stacksize()/_scale['GB']))
     232        msg = ('Resource usage: memory=%.1fMB resident=%.1fMB stacksize=%.1fMB'
     233               % (memory()/_scale['MB'], resident()/_scale['MB'],
     234                  stacksize()/_scale['MB']))
    219235        log(level, msg)
    220236    else:
    221         pass
     237        msg = ('Sorry, no memory statistics for Windows (yet).')
     238        log(level, msg)
     239
     240
     241if __name__ == '__main__':
     242##    critical('Testing exception capturing')
     243    def test_it(num=100):
     244        if num > 0:
     245            test_it(num-1)
     246        else:
     247            resource_usage()
     248
     249    import numpy as num
     250   
     251    a = num.zeros((1000,1000), num.float)
     252
     253    test_it()

  • branches/numpy/anuga/utilities/numerical_tools.py

    r6462 r7035  
    8383    # Compute angle
    8484    p = num.inner(v1, v2)
    85     c = num.inner(v1, normal_vector(v2)) # Projection onto normal
     85    c = num.inner(v1, normal_vector(v2))    # Projection onto normal
    8686                                            # (negative cross product)
    8787       
    … …  
    238238def ensure_numeric(A, typecode=None):
    239239    """Ensure that sequence is a numeric array.
     240
    240241    Inputs:
    241242        A: Sequence. If A is already a numeric array it will be returned

  • branches/numpy/anuga/utilities/system_tools.py

    r7024 r7035  
    1111import getpass
    1212import tarfile
    13 import md5
     13try:
     14    import hashlib
     15except ImportError:
     16    import md5 as hashlib
    1417
    1518
    … …  
    556559def get_file_hexdigest(filename, blocksize=1024*1024*10):
    557560    '''Get a hex digest of a file.'''
    558    
    559     m = md5.new()
     561
     562    if hashlib.__name__ == 'hashlib':
     563        m = hashlib.md5()       # new - 'hashlib' module
     564    else:
     565        m = hashlib.new()       # old - 'md5' module - remove once py2.4 gone
    560566    fd = open(filename, 'r')
    561567           

  • branches/numpy/anuga/utilities/test_polygon.py

    r6689 r7035  
    15051505
    15061506
    1507     def zzztest_inside_polygon_main(self):
     1507    def test_inside_polygon_main(self):
    15081508        # FIXME (Ole): Why is this disabled?
    1509         print "inside",inside
    1510         print "outside",outside
     1509        #print "inside",inside
     1510        #print "outside",outside
    15111511
    15121512        assert not inside_polygon((0.5, 1.5), polygon)
Note: See TracChangeset for help on using the changeset viewer.