Changeset 6689


Ignore:
Timestamp:
Apr 1, 2009, 3:19:07 PM (15 years ago)
Author:
rwilson
Message:

Back-merge from Numeric trunk to numpy branch.

Location:
branches/numpy/anuga
Files:
4 added
27 edited

Legend:

Unmodified
Added
Removed
  • branches/numpy/anuga/abstract_2d_finite_volumes/general_mesh.py

    r6481 r6689  
    221221
    222222        return self.normals[i, 2*j:2*j+2]
     223       
     224    def get_edgelength(self, i, j):
     225        """Return length of j'th edge of the i'th triangle.
     226        Return value is the numeric array slice [x, y]
     227        """
     228        return self.edgelengths[i, j]
     229               
    223230
    224231    def get_number_of_nodes(self):
  • branches/numpy/anuga/abstract_2d_finite_volumes/pmesh2domain.py

    r6304 r6689  
    165165                     verbose=False):
    166166    """Convert a pmesh file or a pmesh mesh instance to a bunch of lists
    167     that can be used to instanciate a domain object.
     167    that can be used to instantiate a domain object.
    168168    """
    169169
     
    180180    volumes = mesh_dict['triangles']
    181181    vertex_quantity_dict = {}
     182
     183    # num.transpose(None) gives scalar array of value None
    182184    point_atts = mesh_dict['vertex_attributes']
    183     point_titles  = mesh_dict['vertex_attribute_titles']
    184     geo_reference  = mesh_dict['geo_reference']
     185
     186    point_titles = mesh_dict['vertex_attribute_titles']
     187    geo_reference = mesh_dict['geo_reference']
    185188    if point_atts is not None:
    186         point_atts = num.transpose(mesh_dict['vertex_attributes'])
     189        point_atts = num.transpose(point_atts)
    187190        for quantity, value_vector in map(None, point_titles, point_atts):
    188191            vertex_quantity_dict[quantity] = value_vector
  • branches/numpy/anuga/abstract_2d_finite_volumes/quantity.py

    r6553 r6689  
    11881188        msg = 'Indices must be a list or None'
    11891189        assert indices is None or isinstance(indices, (list, num.ndarray)), msg
    1190 #        assert type(indices) in [types.ListType, types.NoneType,
    1191 #                                 num.ndarray], \
    1192 #                   'Indices must be a list or None'     #??#
    11931190
    11941191        if location == 'centroids':
  • branches/numpy/anuga/abstract_2d_finite_volumes/test_general_mesh.py

    r6533 r6689  
    132132        domain = General_mesh(nodes, triangles)
    133133
    134 #        value = [7]
    135134        msg = ("domain.get_triangles()=\n%s\nshould be the same as "
    136135               "'triangles'=\n%s"
     
    354353                              nodes, triangles, geo_reference=geo)
    355354
    356     def test_raw_change_points_geo_ref(self):
    357         x0 = 1000.0
    358         y0 = 2000.0
    359         geo = Geo_reference(56, x0, y0)
    360        
    361        
    362 
    363 
    364 #-------------------------------------------------------------
     355################################################################################
    365356
    366357if __name__ == "__main__":
  • branches/numpy/anuga/abstract_2d_finite_volumes/test_region.py

    r6441 r6689  
    6868                    [ 0.07,  0.07,  0.07],
    6969                    [ 0.07,  0.07,  0.07],
    70                     [ 1.0,  1.0,  1.0],
    71                     [ 1.0,  1.0,  1.0]]
     70                    [ 1.0,   1.0,   1.0],
     71                    [ 1.0,   1.0,   1.0]]
    7272        msg = ("\ndomain.quantities['friction']=%s\nexpected value=%s"
    7373               % (str(domain.quantities['friction'].get_values()),
  • branches/numpy/anuga/abstract_2d_finite_volumes/util.py

    r6553 r6689  
    21222122                #add tide to stage if provided
    21232123                if quantity == 'stage':
    2124                      quantity_value[quantity] = num.array(quantity_value[quantity], num.float) \
    2125                                                           + directory_add_tide
     2124                     quantity_value[quantity] = num.array(quantity_value[quantity],
     2125                                                          num.float) + directory_add_tide
    21262126
    21272127                #condition to find max and mins for all the plots
  • branches/numpy/anuga/advection/test_advection.py

    r6304 r6689  
    142142
    143143        X = domain.quantities['stage'].explicit_update
    144         #        print 'X=%s' % str(X)
    145144        assert X[0] == -X[1]
    146145
  • branches/numpy/anuga/caching/test_caching.py

    r6553 r6689  
    390390        if os.name == 'posix' and os.uname()[4] in ['x86_64', 'ia64']:
    391391          # 64 bit hash values
    392           f1hash = 1914027059797211698
    393           f2hash = 1914027059807087171
     392          f1hash = 7079146893884768701
     393          f2hash = -6995306676314913340
     394
     395          # Prior to cluster upgrades Feb 2009
     396          #f1hash = 1914027059797211698
     397          #f2hash = 1914027059807087171
    394398        else:
     399          # 32 bit hash values
    395400          f1hash = -758136387
    396401          f2hash = -11221564     
  • branches/numpy/anuga/coordinate_transforms/geo_reference.py

    r6553 r6689  
    258258        is_list = isinstance(points, list)
    259259
    260         #points = ensure_numeric(copy.copy(points), num.float)
    261260        points = ensure_numeric(points, num.float)
    262261
  • branches/numpy/anuga/culvert_flows/culvert_routines.py

    r6553 r6689  
    2828                                   manning=0.0,
    2929                                   sum_loss=0.0,
     30                                   max_velocity=10.0,
    3031                                   log_filename=None):
    3132
     
    4748
    4849
    49     if inlet_depth > 0.01:
     50    if inlet_depth > 0.1: #this value was 0.01:
    5051        # Water has risen above inlet
    5152       
     
    116117
    117118                # Outlet control velocity using tail water
    118                 culvert_velocity = sqrt(delta_total_energy/((sum_loss/2*g)+(manning**2*culvert_length)/hyd_rad**1.33333))
     119                culvert_velocity = sqrt(delta_total_energy/((sum_loss/2/g)+(manning**2*culvert_length)/hyd_rad**1.33333))
    119120                Q_outlet_tailwater = flow_area * culvert_velocity
    120121               
     
    183184
    184185                # Outlet control velocity using tail water
    185                 culvert_velocity = sqrt(delta_total_energy/((sum_loss/2*g)+(manning**2*culvert_length)/hyd_rad**1.33333))
     186                culvert_velocity = sqrt(delta_total_energy/((sum_loss/2/g)+(manning**2*culvert_length)/hyd_rad**1.33333))
    186187                Q_outlet_tailwater = flow_area * culvert_velocity
    187188
     
    216217    else: # inlet_depth < 0.01:
    217218        Q = barrel_velocity = outlet_culvert_depth = 0.0
     219    # Temporary flow limit
     220    if barrel_velocity > max_velocity:
     221        if log_filename is not None:                           
     222            s = 'Barrel velocity was reduced from = %f m/s to %f m/s' % (barrel_velocity, max_velocity)
     223            log_to_file(log_filename, s)
     224       
     225        barrel_velocity = max_velocity
     226        Q = flow_area * barrel_velocity
     227       
     228       
     229
    218230       
    219231    return Q, barrel_velocity, outlet_culvert_depth
  • branches/numpy/anuga/culvert_flows/test_culvert_routines.py

    r6553 r6689  
    1919
    2020
    21     def NOtest_boyd_1(self):
     21    def test_boyd_1(self):
    2222        """test_boyd_1
    2323       
     
    2929        culvert_slope=0.1  # Downward
    3030
    31         inlet_depth=0.1
    32         outlet_depth=0.09
     31        inlet_depth=2.0
     32        outlet_depth=0.0
    3333
    3434        culvert_length=4.0
     
    5959                                                 sum_loss)
    6060       
    61         print Q, v, d
    62         assert num.allclose(Q, 0.1)
    63         assert num.allclose(v, 0.93)
    64         assert num.allclose(d, 0.09)
     61        #print Q, v, d
     62        assert num.allclose(Q, 3.118, rtol=1.0e-3)
    6563       
    6664
     65        #assert num.allclose(v, 0.93)
     66        #assert num.allclose(d, 0.0)
     67       
     68
     69    def test_boyd_2(self):
     70        """test_boyd_2
     71       
     72        This tests the Boyd routine with data obtained from ??? by Petar Milevski   
     73        """
     74        # FIXME(Ole): This test fails (20 Feb 2009)
     75
     76        g=9.81
     77        culvert_slope=0.1  # Downward
     78
     79        inlet_depth=0.2
     80        outlet_depth=0.0
     81
     82        culvert_length=4.0
     83        culvert_width=1.2
     84        culvert_height=0.75
     85
     86        culvert_type='box'
     87        manning=0.013
     88        sum_loss=0.0
     89
     90        inlet_specific_energy=inlet_depth #+0.5*v**2/g
     91        z_in = 0.0
     92        z_out = -culvert_length*culvert_slope/100
     93        E_in = z_in+inlet_depth # +
     94        E_out = z_out+outlet_depth # +
     95        delta_total_energy = E_in-E_out
     96
     97        Q, v, d = boyd_generalised_culvert_model(inlet_depth,
     98                                                 outlet_depth,
     99                                                 inlet_specific_energy,
     100                                                 delta_total_energy,
     101                                                 g,
     102                                                 culvert_length,
     103                                                 culvert_width,
     104                                                 culvert_height,
     105                                                 culvert_type,
     106                                                 manning,
     107                                                 sum_loss)
     108       
     109        #print Q, v, d
     110        #assert num.allclose(Q, 0.185, rtol=1.0e-3)
     111        #assert num.allclose(v, 0.93)
     112        #assert num.allclose(d, 0.0)
     113       
    67114   
    68115               
  • branches/numpy/anuga/fit_interpolate/fit.py

    r6553 r6689  
    260260            # The memory damage has been done by now.
    261261        else:
    262              AtA = self.AtA #Did this for speed, did ~nothing
     262             AtA = self.AtA # Did this for speed, did ~nothing
    263263        self.point_count += point_coordinates.shape[0]
    264264
     
    267267                                        self.mesh_boundary_polygon,
    268268                                        closed=True,
    269                                         verbose=False) # Too much output if True
    270 
     269                                        verbose=False) # Suppress output
    271270       
    272271        n = len(inside_indices)
     
    283282           
    284283            if element_found is True:
    285                 j0 = triangles[k,0] #Global vertex id for sigma0
    286                 j1 = triangles[k,1] #Global vertex id for sigma1
    287                 j2 = triangles[k,2] #Global vertex id for sigma2
     284                j0 = triangles[k,0] # Global vertex id for sigma0
     285                j1 = triangles[k,1] # Global vertex id for sigma1
     286                j2 = triangles[k,2] # Global vertex id for sigma2
    288287
    289288                sigmas = {j0:sigma0, j1:sigma1, j2:sigma2}
     
    304303                                         self.mesh_boundary_polygon,
    305304                                         closed=True,
    306                                          verbose=False) # Too much output if True               
     305                                         verbose=False) # Suppress output
    307306                msg = 'Point (%f, %f) is not inside mesh boundary' % tuple(x)
    308307                assert flag is True, msg               
     
    318317                raise RuntimeError, msg
    319318
    320             self.AtA = AtA
     319               
     320        self.AtA = AtA
    321321
    322322       
  • branches/numpy/anuga/fit_interpolate/interpolate.py

    r6553 r6689  
    835835                raise msg
    836836
     837            # Ensure 'mesh_boundary_polygon' is defined
     838            mesh_boundary_polygon = None
     839           
    837840            if triangles is not None and vertex_coordinates is not None:
    838841                # Check that all interpolation points fall within
     
    894897                                'requires the latter.')
    895898
    896             # Plot boundary and interpolation points
    897             if verbose is True:
     899            # Plot boundary and interpolation points,
     900            # but only if if 'mesh_boundary_polygon' has data.
     901            if verbose is True and mesh_boundary_polygon is not None:
    898902                import sys
    899903                if sys.platform == 'win32':
  • branches/numpy/anuga/fit_interpolate/test_interpolate.py

    r6441 r6689  
    919919        triangles = [[0,1,3], [1,0,2], [0,4,5], [0,5,2]] #abd bac aef afc
    920920
    921         point_coords_absolute = [[-2.0, 2.0],
    922                                  [-1.0, 1.0],
    923                                  [0.0, 2.0],
    924                                  [1.0, 1.0],
    925                                  [2.0, 1.0],
    926                                  [0.0, 0.0],
    927                                  [1.0, 0.0],
    928                                  [0.0, -1.0],
     921        point_coords_absolute = [[-2.0,  2.0],
     922                                 [-1.0,  1.0],
     923                                 [ 0.0, 2.0],
     924                                 [ 1.0, 1.0],
     925                                 [ 2.0, 1.0],
     926                                 [ 0.0, 0.0],
     927                                 [ 1.0, 0.0],
     928                                 [ 0.0, -1.0],
    929929                                 [-0.2, -0.5],
    930930                                 [-0.9, -1.5],
    931                                  [0.5, -1.9],
    932                                  [3.0, 1.0]]
     931                                 [ 0.5, -1.9],
     932                                 [ 3.0, 1.0]]
    933933
    934934        geo = Geo_reference(57, 100, 500)
     
    967967
    968968
    969         point_coords = [[-2.0, 2.0],
    970                         [-1.0, 1.0],
    971                         [0.0, 2.0],
    972                         [1.0, 1.0],
    973                         [2.0, 1.0],
    974                         [0.0, 0.0],
    975                         [1.0, 0.0],
    976                         [0.0, -1.0],
     969        point_coords = [[-2.0,  2.0],
     970                        [-1.0,  1.0],
     971                        [ 0.0, 2.0],
     972                        [ 1.0, 1.0],
     973                        [ 2.0, 1.0],
     974                        [ 0.0, 0.0],
     975                        [ 1.0, 0.0],
     976                        [ 0.0, -1.0],
    977977                        [-0.2, -0.5],
    978978                        [-0.9, -1.5],
    979                         [0.5, -1.9],
    980                         [3.0, 1.0]]
     979                        [ 0.5, -1.9],
     980                        [ 3.0, 1.0]]
    981981
    982982        interp = Interpolate(vertices, triangles)
     
    10301030
    10311031
    1032         point_coords = [[-2.0, 2.0],
    1033                         [-1.0, 1.0],
    1034                         [0.0, 2.0],
    1035                         [1.0, 1.0],
    1036                         [2.0, 1.0],
    1037                         [0.0, 0.0],
    1038                         [1.0, 0.0],
    1039                         [0.0, -1.0],
     1032        point_coords = [[-2.0,  2.0],
     1033                        [-1.0,  1.0],
     1034                        [ 0.0, 2.0],
     1035                        [ 1.0, 1.0],
     1036                        [ 2.0, 1.0],
     1037                        [ 0.0, 0.0],
     1038                        [ 1.0, 0.0],
     1039                        [ 0.0, -1.0],
    10401040                        [-0.2, -0.5],
    10411041                        [-0.9, -1.5],
    1042                         [0.5, -1.9],
    1043                         [3.0, 1.0]]
     1042                        [ 0.5, -1.9],
     1043                        [ 3.0, 1.0]]
    10441044
    10451045        interp = Interpolate(vertices, triangles)
  • branches/numpy/anuga/fit_interpolate/test_search_functions.py

    r6553 r6689  
    113113                if k >= 0:
    114114                    V = mesh.get_vertex_coordinates(k) # nodes for triangle k
     115                    assert is_inside_triangle(x, V, closed=True)
    115116                    assert is_inside_polygon(x, V)
    116117                    assert found is True
  • branches/numpy/anuga/geospatial_data/geospatial_data.py

    r6533 r6689  
    1111from RandomArray import randint, seed, get_seed
    1212from copy import deepcopy
     13import copy
    1314
    1415from Scientific.IO.NetCDF import NetCDFFile
     
    14181419        assert geo_reference == None, msg
    14191420    else:
    1420         points = ensure_numeric(points, num.float)
     1421        points = ensure_numeric(copy.copy(points), num.float)
    14211422
    14221423    # Sort of geo_reference and convert points
  • branches/numpy/anuga/pmesh/mesh_interface.py

    r6304 r6689  
    147147                    found = True
    148148            if found is False:
    149                 msg = 'Segment %d was not asigned a boundary_tag' %i
    150                 raise Exception, msg
     149                msg = 'Segment %d was not assigned a boundary_tag.' % i
     150                msg +=  'Default tag "exterior" will be assigned to missing segment'
     151                #raise Exception, msg
     152                # Fixme: Use proper Python warning
     153                if verbose: print 'WARNING: ', msg
    151154               
    152155
  • branches/numpy/anuga/shallow_water/data_manager.py

    r6553 r6689  
    9191from anuga.utilities.polygon import intersection
    9292
     93from anuga.utilities.system_tools import get_vars_in_expression
     94
     95
     96# Default block size for sww2dem()
     97DEFAULT_BLOCK_SIZE = 10000
    9398
    9499######
     
    21062111            origin=None,
    21072112            datum='WGS84',
    2108             format='ers'):
     2113            format='ers',
     2114            block_size=None):
    21092115    """Read SWW file and convert to Digitial Elevation model format
    21102116    (.asc or .ers)
     
    21492155
    21502156    format can be either 'asc' or 'ers'
     2157    block_size - sets the number of slices along the non-time axis to
     2158                 process in one block.
    21512159    """
    21522160
     
    21792187        number_of_decimal_places = 3
    21802188
     2189    if block_size is None:
     2190        block_size = DEFAULT_BLOCK_SIZE
     2191
     2192    # Read SWW file
    21812193    swwfile = basename_in + '.sww'
    21822194    demfile = basename_out + '.' + format
     
    22552267            if verbose: print '    %s in [%f, %f]' %(name, min(q), max(q))
    22562268
    2257     # Get quantity and reduce if applicable
    2258     if verbose: print 'Processing quantity %s' %quantity
    2259 
    2260     # Turn NetCDF objects into numeric arrays
    2261     try:
    2262         q = fid.variables[quantity][:]
    2263     except:
    2264         quantity_dict = {}
    2265         for name in fid.variables.keys():
    2266             quantity_dict[name] = fid.variables[name][:]
    2267         #Convert quantity expression to quantities found in sww file
    2268         q = apply_expression_to_dictionary(quantity, quantity_dict)
    2269 
    2270     if len(q.shape) == 2:
    2271         #q has a time component, must be reduced alongthe temporal dimension
    2272         if verbose: print 'Reducing quantity %s' %quantity
    2273         q_reduced = num.zeros(number_of_points, num.float)
    2274 
    2275         if timestep is not None:
    2276             for k in range(number_of_points):
    2277                 q_reduced[k] = q[timestep,k]
    2278         else:
    2279             for k in range(number_of_points):
    2280                 q_reduced[k] = reduction( q[:,k] )
    2281 
    2282         q = q_reduced
    2283 
     2269    # Get the variables in the supplied expression.
     2270    # This may throw a SyntaxError exception.
     2271    var_list = get_vars_in_expression(quantity)
     2272
     2273    # Check that we have the required variables in the SWW file.
     2274    missing_vars = []
     2275    for name in var_list:
     2276        try:
     2277            _ = fid.variables[name]
     2278        except:
     2279            missing_vars.append(name)
     2280    if missing_vars:
     2281        msg = ("In expression '%s', variables %s are not in the SWW file '%s'"
     2282               % (quantity, swwfile))
     2283        raise Exception, msg
     2284
     2285    # Create result array and start filling, block by block.
     2286    result = num.zeros(number_of_points, num.float)
     2287
     2288    for start_slice in xrange(0, number_of_points, block_size):
     2289        # limit slice size to array end if at last block
     2290        end_slice = min(start_slice + block_size, number_of_points)
     2291       
     2292        # get slices of all required variables
     2293        q_dict = {}
     2294        for name in var_list:
     2295            # check if variable has time axis
     2296            if len(fid.variables[name].shape) == 2:
     2297                q_dict[name] = fid.variables[name][:,start_slice:end_slice]
     2298            else:       # no time axis
     2299                q_dict[name] = fid.variables[name][start_slice:end_slice]
     2300
     2301        # Evaluate expression with quantities found in SWW file
     2302        res = apply_expression_to_dictionary(quantity, q_dict)
     2303
     2304        if len(res.shape) == 2:
     2305            new_res = num.zeros(res.shape[1], num.float)
     2306            for k in xrange(res.shape[1]):
     2307                new_res[k] = reduction(res[:,k])
     2308            res = new_res
     2309
     2310        result[start_slice:end_slice] = res
     2311                                   
    22842312    #Post condition: Now q has dimension: number_of_points
    2285     assert len(q.shape) == 1
    2286     assert q.shape[0] == number_of_points
     2313    assert len(result.shape) == 1
     2314    assert result.shape[0] == number_of_points
    22872315
    22882316    if verbose:
    22892317        print 'Processed values for %s are in [%f, %f]' % \
    2290               (quantity, min(q), max(q))
     2318              (quantity, min(result), max(result))
    22912319
    22922320    #Create grid and update xll/yll corner and x,y
     
    23162344    assert xmax >= xmin, msg
    23172345
    2318     msg = 'yax must be greater than or equal to xmin.\n'
     2346    msg = 'ymax must be greater than or equal to xmin.\n'
    23192347    msg += 'I got ymin = %f, ymax = %f' %(ymin, ymax)
    23202348    assert ymax >= ymin, msg
     
    23612389    #Interpolate using quantity values
    23622390    if verbose: print 'Interpolating'
    2363     grid_values = interp.interpolate(q, grid_points).flatten()
     2391    grid_values = interp.interpolate(result, grid_points).flatten()
    23642392
    23652393    if verbose:
     
    26722700
    26732701    if verbose:
    2674         print 'Interpolated values are in [%f, %f]' \
    2675               % (num.min(interpolated_values), num.max(interpolated_values))
     2702        print ('Interpolated values are in [%f, %f]'
     2703               % (num.min(interpolated_values), num.max(interpolated_values)))
    26762704
    26772705    # Assign NODATA_value to all points outside bounding polygon
     
    54755503            origin=None,
    54765504            zone=None,
     5505            central_meridian=None,           
    54775506            mean_stage=0.0,
    54785507            zscale=1.0,
     
    55145543    output:
    55155544      basename_out: name of sts file in which mux2 data is stored.
     5545     
     5546     
     5547     
     5548    NOTE: South is positive in mux files so sign of y-component of velocity is reverted
    55165549    """
    55175550
     
    56875720    # Check zone boundaries
    56885721    if zone is None:
    5689         refzone, _, _ = redfearn(latitudes[0], longitudes[0])
     5722        refzone, _, _ = redfearn(latitudes[0], longitudes[0],
     5723                                 central_meridian=central_meridian)
    56905724    else:
    56915725        refzone = zone
     
    56945728
    56955729    for i in range(number_of_points):
    5696         zone, easting, northing = redfearn(latitudes[i], longitudes[i],
    5697                                            zone=zone)
     5730        computed_zone, easting, northing = redfearn(latitudes[i], longitudes[i],
     5731                                                    zone=zone,
     5732                                                    central_meridian=central_meridian)
    56985733        x[i] = easting
    56995734        y[i] = northing
    5700         if zone != refzone:
     5735        if computed_zone != refzone:
    57015736            msg = 'All sts gauges need to be in the same zone. \n'
    57025737            msg += 'offending gauge:Zone %d,%.4f, %4f\n' \
    5703                    % (zone, easting, northing)
     5738                   % (computed_zone, easting, northing)
    57045739            msg += 'previous gauge:Zone %d,%.4f, %4f' \
    57055740                   % (old_zone, old_easting, old_northing)
    57065741            raise Exception, msg
    5707         old_zone = zone
     5742        old_zone = computed_zone
    57085743        old_easting = easting
    57095744        old_northing = northing
     
    57445779
    57455780            xmomentum[j,i] = ua * h
    5746             ymomentum[j,i] = va * h
     5781            ymomentum[j,i] = -va * h # South is positive in mux files
     5782
    57475783
    57485784    outfile.close()
     
    58635899        # This is being used to seperate one number from a list.
    58645900        # what it is actually doing is sorting lists from numeric arrays.
    5865         if type(times) is list or isinstance(times, num.ndarray):
     5901        if isinstance(times, (list, num.ndarray)):
    58665902            number_of_times = len(times)
    58675903            times = ensure_numeric(times)
     
    59305966            #outfile.variables[q+Write_sww.RANGE][1] = -max_float # Max
    59315967
    5932         if type(times) is list or isinstance(times, num.ndarray):
     5968        if isinstance(times, (list, num.ndarray)):
    59335969            outfile.variables['time'][:] = times    #Store time relative
    59345970
     
    62656301        # This is being used to seperate one number from a list.
    62666302        # what it is actually doing is sorting lists from numeric arrays.
    6267         if type(times) is list or isinstance(times, num.ndarray):
     6303        if isinstance(times, (list, num.ndarray)):
    62686304            number_of_times = len(times)
    62696305            times = ensure_numeric(times)
     
    63126348            outfile.variables[q + Write_sts.RANGE][1] = -max_float # Max
    63136349
    6314         if type(times) is list or isinstance(times, num.ndarray):
     6350        if isinstance(times, (list, num.ndarray)):
    63156351            outfile.variables['time'][:] = times    #Store time relative
    63166352
  • branches/numpy/anuga/shallow_water/shallow_water_domain.py

    r6553 r6689  
    812812
    813813        return msg
    814 
     814       
     815       
     816
     817    def compute_boundary_flows(self):
     818        """Compute boundary flows at current timestep.
     819                       
     820        Quantities computed are:
     821           Total inflow across boundary
     822           Total outflow across boundary
     823           Flow across each tagged boundary segment
     824        """
     825               
     826        # Run through boundary array and compute for each segment
     827        # the normal momentum ((uh, vh) dot normal) times segment length.
     828        # Based on sign accumulate this into boundary_inflow and boundary_outflow.
     829                       
     830        # Compute flows along boundary
     831       
     832        uh = self.get_quantity('xmomentum').get_values(location='edges')
     833        vh = self.get_quantity('ymomentum').get_values(location='edges')       
     834       
     835        # Loop through edges that lie on the boundary and calculate
     836        # flows
     837        boundary_flows = {}
     838        total_boundary_inflow = 0.0
     839        total_boundary_outflow = 0.0
     840        for vol_id, edge_id in self.boundary:
     841            # Compute normal flow across edge. Since normal vector points
     842            # away from triangle, a positive sign means that water
     843            # flows *out* from this triangle.
     844             
     845            momentum = [uh[vol_id, edge_id], vh[vol_id, edge_id]]
     846            normal = self.mesh.get_normal(vol_id, edge_id)
     847            length = self.mesh.get_edgelength(vol_id, edge_id)           
     848            normal_flow = num.dot(momentum, normal)*length
     849           
     850            # Reverse sign so that + is taken to mean inflow
     851            # and - means outflow. This is more intuitive.
     852            edge_flow = -normal_flow
     853           
     854            # Tally up inflows and outflows separately
     855            if edge_flow > 0:
     856                # Flow is inflow     
     857                total_boundary_inflow += edge_flow                                 
     858            else:
     859                # Flow is outflow
     860                total_boundary_outflow += edge_flow   
     861
     862            # Tally up flows by boundary tag
     863            tag = self.boundary[(vol_id, edge_id)]
     864           
     865            if tag not in boundary_flows:
     866                boundary_flows[tag] = 0.0
     867            boundary_flows[tag] += edge_flow
     868           
     869               
     870        return boundary_flows, total_boundary_inflow, total_boundary_outflow
     871       
     872
     873    def compute_forcing_flows(self):
     874        """Compute flows in and out of domain due to forcing terms.
     875                       
     876        Quantities computed are:
     877               
     878       
     879           Total inflow through forcing terms
     880           Total outflow through forcing terms
     881           Current total volume in domain       
     882
     883        """
     884
     885        #FIXME(Ole): We need to separate what part of explicit_update was
     886        # due to the normal flux calculations and what is due to forcing terms.
     887       
     888        pass
     889                       
     890       
     891    def compute_total_volume(self):
     892        """Compute total volume (m^3) of water in entire domain
     893        """
     894       
     895        area = self.mesh.get_areas()
     896        volume = 0.0
     897       
     898        stage = self.get_quantity('stage').get_values(location='centroids')
     899        elevation = self.get_quantity('elevation').get_values(location='centroids')       
     900        depth = stage-elevation
     901       
     902        #print 'z', elevation
     903        #print 'w', stage
     904        #print 'h', depth
     905        return num.sum(depth*area)
     906       
     907       
     908    def volumetric_balance_statistics(self):               
     909        """Create volumetric balance report suitable for printing or logging.
     910        """
     911       
     912        boundary_flows, total_boundary_inflow, total_boundary_outflow = self.compute_boundary_flows()
     913       
     914        s = '---------------------------\n'       
     915        s += 'Volumetric balance report:\n'
     916        s += '--------------------------\n'
     917        s += 'Total boundary inflow [m^3/s]: %.2f\n' % total_boundary_inflow
     918        s += 'Total boundary outflow [m^3/s]: %.2f\n' % total_boundary_outflow       
     919        s += 'Net boundary flow by tags [m^3/s]\n'
     920        for tag in boundary_flows:
     921            s += '    %s [m^3/s]: %.2f\n' % (tag, boundary_flows[tag])
     922       
     923        s += 'Total net boundary flow [m^3/s]: %.2f\n' % (total_boundary_inflow + total_boundary_outflow)
     924        s += 'Total volume in domain [m^3]: %.2f\n' % self.compute_total_volume()
     925       
     926        # The go through explicit forcing update and record the rate of change for stage and
     927        # record into forcing_inflow and forcing_outflow. Finally compute integral
     928        # of depth to obtain total volume of domain.
     929       
     930        # FIXME(Ole): This part is not yet done.               
     931       
     932        return s       
     933           
    815934################################################################################
    816935# End of class Shallow Water Domain
     
    13861505        domain: pointer to shallow water domain for which the boundary applies
    13871506        mean_stage: The mean water level which will be added to stage derived
    1388                     from the sww file
     1507                    from the boundary condition
    13891508        time_thinning: Will set how many time steps from the sww file read in
    13901509                       will be interpolated to the boundary. For example if
     
    14021521                          that available in the underlying data.
    14031522
     1523                          Note that mean_stage will also be added to this.
     1524                                               
    14041525        use_cache:
    14051526        verbose:
     
    15741695                xmom_update[k] += S*uh[k]
    15751696                ymom_update[k] += S*vh[k]
     1697
     1698def depth_dependent_friction(domain, default_friction,
     1699                             surface_roughness_data,
     1700                             verbose=False):
     1701    """Returns an array of friction values for each wet element adjusted for depth.
     1702
     1703    Inputs:
     1704        domain - computational domain object
     1705        default_friction - depth independent bottom friction
     1706        surface_roughness_data - N x 5 array of n0, d1, n1, d2, n2 values for each
     1707        friction region.
     1708
     1709    Outputs:
     1710        wet_friction - Array that can be used directly to update friction as follows:
     1711                       domain.set_quantity('friction', wet_friction)
     1712
     1713       
     1714       
     1715    """
     1716
     1717    import Numeric as num
     1718   
     1719    # Create a temp array to store updated depth dependent friction for wet elements
     1720    # EHR this is outwardly inneficient but not obvious how to avoid recreating each call??????
     1721    N=len(domain)
     1722    wet_friction    = num.zeros(N, num.Float)
     1723    wet_friction[:] = default_n0   # Initially assign default_n0 to all array so sure have no zeros values
     1724   
     1725   
     1726    depth = domain.create_quantity_from_expression('stage - elevation')  # create depth instance for this timestep
     1727    # Recompute depth as vector 
     1728    d = depth.get_values(location='centroids')
     1729 
     1730    # rebuild the 'friction' values adjusted for depth at this instant
     1731    for i in domain.get_wet_elements():                                  # loop for each wet element in domain
     1732       
     1733        # Get roughness data for each element
     1734        n0 = float(surface_roughness_data[i,0])
     1735        d1 = float(surface_roughness_data[i,1])
     1736        n1 = float(surface_roughness_data[i,2])
     1737        d2 = float(surface_roughness_data[i,3])
     1738        n2 = float(surface_roughness_data[i,4])
     1739       
     1740       
     1741        # Recompute friction values from depth for this element
     1742               
     1743        if d[i]   <= d1:
     1744            depth_dependent_friction = n1
     1745        elif d[i] >= d2:
     1746            depth_dependent_friction = n2
     1747        else:
     1748            depth_dependent_friction = n1+((n2-n1)/(d2-d1))*(d[i]-d1)
     1749           
     1750        # check sanity of result
     1751        if (depth_dependent_friction  < 0.010 or depth_dependent_friction > 9999.0) :
     1752            print model_data.basename+' >>>> WARNING: computed depth_dependent friction out of range ddf,n1,n2 ', depth_dependent_friction, n1,n2
     1753       
     1754        # update depth dependent friction  for that wet element
     1755        wet_friction[i] = depth_dependent_friction
     1756       
     1757    # EHR add code to show range of 'friction across domain at this instant as sanity check?????????
     1758   
     1759    if verbose :
     1760        nvals=domain.get_quantity('friction').get_values(location='centroids')        # return array of domain nvals
     1761        n_min=min(nvals)
     1762        n_max=max(nvals)
     1763       
     1764        print "         ++++ calculate_depth_dependent_friction  - Updated friction - range  %7.3f to %7.3f" %(n_min,n_max)
     1765   
     1766    return wet_friction
     1767
    15761768
    15771769################################################################################
     
    19012093                assert is_inside_polygon(point, bounding_polygon), msg
    19022094
    1903             # Compute area and check that it is greater than 0
    1904             self.exchange_area = polygon_area(self.polygon)
    1905 
    1906             msg = 'Polygon %s in forcing term' % str(self.polygon)
    1907             msg += ' has area = %f' % self.exchange_area
    1908             assert self.exchange_area > 0.0
    1909 
    19102095        # Pointer to update vector
    19112096        self.update = domain.quantities[self.quantity_name].explicit_update
     
    19312116        if self.polygon is not None:
    19322117            # Inlet is polygon
    1933             inlet_region = ('polygon=\n%s, area=%f m^2' %
    1934                             (self.polygon, self.exchange_area))
    1935 
     2118            inlet_region = 'polygon=%s' % (self.polygon)
    19362119            self.exchange_indices = inside_polygon(points, self.polygon)
    19372120
    1938         if self.exchange_indices is not None:
     2121        if self.exchange_indices is None:
     2122            self.exchange_area = polygon_area(bounding_polygon)
     2123        else:   
    19392124            if len(self.exchange_indices) == 0:
    19402125                msg = 'No triangles have been identified in '
     
    19422127                raise Exception, msg
    19432128
     2129            # Compute exchange area as the sum of areas of triangles identified
     2130            # by circle or polygon
     2131            self.exchange_area = 0.0
     2132            for i in self.exchange_indices:
     2133                self.exchange_area += domain.areas[i]
     2134           
     2135
     2136        msg = 'Exchange area in forcing term'
     2137        msg += ' has area = %f' %self.exchange_area
     2138        assert self.exchange_area > 0.0           
     2139           
     2140               
     2141
     2142           
    19442143        # Check and store default_rate
    19452144        msg = ('Keyword argument default_rate must be either None '
  • branches/numpy/anuga/shallow_water/test_data_manager.py

    r6553 r6689  
    12211221
    12221222    def test_sww2dem_asc_elevation_depth(self):
    1223         """
    1224         test_sww2dem_asc_elevation_depth(self):
     1223        """test_sww2dem_asc_elevation_depth
     1224       
    12251225        Test that sww information can be converted correctly to asc/prj
    12261226        format readable by e.g. ArcView
     1227       
     1228        Also check geo_reference is correct
    12271229        """
    12281230
     
    12491251        sww.store_timestep()
    12501252
    1251         #self.domain.tight_slope_limiters = 1
    12521253
    12531254        self.domain.evolve_to_end(finaltime = 0.01)
     
    20002001                number_of_decimal_places = 9,
    20012002                verbose = self.verbose,
    2002                 format = 'asc')
     2003                format = 'asc',
     2004                block_size=2)
    20032005
    20042006
     
    50645066        #print "longitudes_news",longitudes_news
    50655067
    5066         self.failUnless(latitudes_new == [-35, -40, -45] and \
     5068        self.failUnless(latitudes_new == [-35, -40, -45] and
    50675069                        longitudes_news == [148, 149,150],
    50685070                         'failed')
     
    55195521        na and va quantities will be the Easting values.
    55205522        Depth and ua will be the Northing value.
     5523       
     5524        The mux file format has south as positive so
     5525        this function will swap the sign for va. 
    55215526        """
    55225527
     
    60366041                quantities_init[2].append(num.ones(time_step_count,num.float)*this_va) #
    60376042            else:
    6038                 quantities_init[2].append(va[i])           
     6043                quantities_init[2].append(-va[i]) # South is negative in MUX
    60396044
    60406045        file_handle, base_name = tempfile.mkstemp("")
     
    61586163        assert num.allclose(xvelocity,ua)
    61596164        msg='incorrect gauge va time series returned'
    6160         assert num.allclose(yvelocity,va)
     6165        assert num.allclose(yvelocity, -va)
    61616166
    61626167    def test_urs2sts_read_mux2_pyII(self):
     
    62166221        assert num.allclose(xvelocity,ua)
    62176222        msg='incorrect gauge va time series returned'
    6218         assert num.allclose(yvelocity,va)
     6223        assert num.allclose(yvelocity, -va) # South is positive in MUX
    62196224
    62206225    def test_urs2sts_read_mux2_pyIII(self):
     
    62886293        assert num.allclose(xvelocity,ua)
    62896294        msg='incorrect gauge va time series returned'
    6290         assert num.allclose(yvelocity,va)
     6295        assert num.allclose(yvelocity, -va) # South is positive in mux
    62916296       
    62926297
     
    63906395                if j == 0: assert num.allclose(data[i][:parameters_index], ha0[permutation[i], :])
    63916396                if j == 1: assert num.allclose(data[i][:parameters_index], ua0[permutation[i], :])
    6392                 if j == 2: assert num.allclose(data[i][:parameters_index], va0[permutation[i], :])
     6397                if j == 2: assert num.allclose(data[i][:parameters_index], -va0[permutation[i], :])
    63936398       
    63946399
     
    65966601                        x = unpack('f',f.read(4))[0]
    65976602                        #print time, x, q_time[time, point_i]
     6603                        if q == 'VA': x = -x # South is positive in MUX
    65986604                        assert abs(q_time[time, point_i]-x)<epsilon
    65996605
     
    66696675                    #print 'v ', data[i][:parameters_index][8]                   
    66706676               
    6671                     assert num.allclose(data[i][:parameters_index], va1[permutation[i], :])
     6677                    # South is positive in MUX
     6678                    assert num.allclose(data[i][:parameters_index], -va1[permutation[i], :])
    66726679                   
    66736680       
     
    67876794        os.remove(sts_file)
    67886795
    6789     def test_urs2sts_nonstandard_projection(self):
     6796    def test_urs2sts_nonstandard_meridian(self):
    67906797        """
    6791         Test single source
     6798        Test single source using the meridian from zone 50 as a nonstandard meridian
    67926799        """
    67936800        tide=0
     
    68216828        urs2sts(base_name,
    68226829                basename_out=base_name,
    6823                 zone=50,
    6824                 mean_stage=tide,verbose=False)
     6830                central_meridian=123,
     6831                mean_stage=tide,
     6832                verbose=False)
    68256833
    68266834        # now I want to check the sts file ...
     
    68456853        #Using the non standard projection (50)
    68466854        for i in range(4):
    6847             zone, e, n = redfearn(lat_long_points[i][0], lat_long_points[i][1], zone=50)
     6855            zone, e, n = redfearn(lat_long_points[i][0],
     6856                                  lat_long_points[i][1],
     6857                                  central_meridian=123)
    68486858            assert num.allclose([x[i],y[i]], [e,n])
    6849             assert zone==geo_reference.zone
    6850 
     6859            assert zone==-1
     6860
     6861           
    68516862    def test_urs2sts_nonstandard_projection_reverse(self):
    68526863        """
     
    69076918        #Using the non standard projection (50)
    69086919        for i in range(4):
    6909             zone, e, n = redfearn(lat_long_points[i][0], lat_long_points[i][1], zone=50)
     6920            zone, e, n = redfearn(lat_long_points[i][0], lat_long_points[i][1],
     6921                                  zone=50)
    69106922            assert num.allclose([x[i],y[i]], [e,n])
    69116923            assert zone==geo_reference.zone
     
    71917203                msg = 'urs stage is not equal to sts stage for for source %i and station %i' %(source_number,j)
    71927204                assert num.allclose(urs_stage[index_start_urs_z:index_end_urs_z],
    7193                                 sts_stage[index_start_stage:index_end_stage],
    7194                                 rtol=1.0e-6, atol=1.0e-5 ), msg                               
     7205                                    sts_stage[index_start_stage:index_end_stage],
     7206                                    rtol=1.0e-6, atol=1.0e-5 ), msg                               
    71957207                               
    71967208                # check that urs e velocity and sts xmomentum are the same
     
    72057217                msg = 'urs n velocity is not equivalent to sts y momentum for source %i and station %i' %(source_number,j)
    72067218                assert num.allclose(urs_n[index_start_urs_n:index_end_urs_n]*(urs_stage[index_start_urs_n:index_end_urs_n]-elevation[j]),
    7207                                 sts_ymom[index_start_y:index_end_y],
     7219                                -sts_ymom[index_start_y:index_end_y],
    72087220                                rtol=1.0e-5, atol=1.0e-4 ), msg
    72097221                                               
     
    1059610608        ymomentum = fid.variables['ymomentum'][:]       
    1059710609
     10610       
     10611       
    1059810612        for i in range(stage.shape[0]):
    1059910613            h = stage[i]-z # depth vector at time step i
     
    1064610660
    1064710661        # Check runup restricted to a polygon
    10648         p = num.array([[50,1], [99,1], [99,49], [50,49]], num.int) + \
    10649                 num.array([E, N], num.int)      #array default#
     10662        p = num.array([[50,1], [99,1], [99,49], [50,49]], num.int) + num.array([E, N], num.int)      #array default#
    1065010663
    1065110664        runup = get_maximum_inundation_elevation(swwfile, polygon=p)
  • branches/numpy/anuga/shallow_water/test_shallow_water_domain.py

    r6553 r6689  
    22
    33import unittest, os
     4import os.path
    45from math import pi, sqrt
    56import tempfile
     
    1415from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular_cross
    1516
     17from anuga.utilities.system_tools import get_pathname_from_package
    1618from shallow_water_domain import *
    1719
     
    720722        # Check that flux seen from other triangles is inverse
    721723        (ql, qr) = (qr, ql)
    722         #tmp = qr
    723         #qr = ql
    724         #ql = tmp
    725724        normal = domain.get_normal(2, 2)
    726725        max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux,
     
    738737        # Check that flux seen from other triangles is inverse
    739738        (ql, qr) = (qr, ql)
    740         #tmp = qr
    741         #qr = ql
    742         #ql = tmp
    743739        normal = domain.get_normal(3, 0)
    744740        max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux,
     
    756752        # Check that flux seen from other triangles is inverse
    757753        (ql, qr) = (qr, ql)
    758         #tmp = qr
    759         #qr=ql
    760         #ql=tmp
    761754        normal = domain.get_normal(0, 1)
    762755        max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux,
     
    24722465        R = Rainfall(domain, rate=2.0, polygon=[[1,1], [2,1], [2,2], [1,2]])
    24732466
    2474         assert num.allclose(R.exchange_area, 1)
     2467        assert num.allclose(R.exchange_area, 2)
    24752468
    24762469        domain.forcing_terms.append(R)
     
    25082501        # restricted to a polygon enclosing triangle #1 (bce)
    25092502        domain.forcing_terms = []
    2510         R = Rainfall(domain, rate=lambda t: 3*t + 7,
    2511                      polygon=[[1,1], [2,1], [2,2], [1,2]])
    2512 
    2513         assert num.allclose(R.exchange_area, 1)
    2514 
     2503        R = Rainfall(domain,
     2504                     rate=lambda t: 3*t + 7,
     2505                     polygon = [[1,1], [2,1], [2,2], [1,2]])
     2506
     2507        assert num.allclose(R.exchange_area, 2)
     2508       
    25152509        domain.forcing_terms.append(R)
    25162510
     
    25512545        # restricted to a polygon enclosing triangle #1 (bce)
    25522546        domain.forcing_terms = []
    2553         R = Rainfall(domain, rate=lambda t: 3*t + 7,
    2554                      polygon=rainfall_poly)
    2555 
    2556         assert num.allclose(R.exchange_area, 1)
    2557 
     2547        R = Rainfall(domain,
     2548                     rate=lambda t: 3*t + 7,
     2549                     polygon=rainfall_poly)                     
     2550
     2551        assert num.allclose(R.exchange_area, 2)
     2552       
    25582553        domain.forcing_terms.append(R)
    25592554
     
    26162611        # restricted to a polygon enclosing triangle #1 (bce)
    26172612        domain.forcing_terms = []
    2618         R = Rainfall(domain, rate=lambda t: 3*t + 7, polygon=rainfall_poly)
    2619 
    2620         assert num.allclose(R.exchange_area, 1)
    2621 
     2613        R = Rainfall(domain,
     2614                     rate=lambda t: 3*t + 7,
     2615                     polygon=rainfall_poly)
     2616
     2617        assert num.allclose(R.exchange_area, 2)
     2618       
    26222619        domain.forcing_terms.append(R)
    26232620
     
    26802677
    26812678        domain.forcing_terms = []
    2682         R = Rainfall(domain, rate=main_rate,
    2683                      polygon=[[1,1], [2,1], [2,2], [1,2]], default_rate=5.0)
    2684 
    2685         assert num.allclose(R.exchange_area, 1)
    2686 
     2679        R = Rainfall(domain,
     2680                     rate=main_rate,
     2681                     polygon = [[1,1], [2,1], [2,2], [1,2]],
     2682                     default_rate=5.0)
     2683
     2684        assert num.allclose(R.exchange_area, 2)
     2685       
    26872686        domain.forcing_terms.append(R)
    26882687
     
    27462745
    27472746        domain.forcing_terms = []
    2748         R = Rainfall(domain, rate=main_rate,
    2749                      polygon=[[1,1], [2,1], [2,2], [1,2]], default_rate=5.0)
    2750 
    2751         assert num.allclose(R.exchange_area, 1)
    2752 
     2747        R = Rainfall(domain,
     2748                     rate=main_rate,
     2749                     polygon=[[1,1], [2,1], [2,2], [1,2]],
     2750                     default_rate=5.0)
     2751
     2752        assert num.allclose(R.exchange_area, 2)
     2753       
    27532754        domain.forcing_terms.append(R)
    27542755
     
    27852786        # on a circle affecting triangles #0 and #1 (bac and bce)
    27862787        domain.forcing_terms = []
    2787         domain.forcing_terms.append(Inflow(domain, rate=2.0,
    2788                                            center=(1,1), radius=1))
    2789 
     2788       
     2789        I = Inflow(domain, rate=2.0, center=(1,1), radius=1)
     2790        domain.forcing_terms.append(I)       
    27902791        domain.compute_forcing_terms()
    27912792
    2792         assert num.allclose(domain.quantities['stage'].explicit_update[1],
    2793                             2.0/pi)
    2794         assert num.allclose(domain.quantities['stage'].explicit_update[0],
    2795                             2.0/pi)
    2796         assert num.allclose(domain.quantities['stage'].explicit_update[2:], 0)
     2793       
     2794        A = I.exchange_area
     2795        assert num.allclose(A, 4) # Two triangles       
     2796       
     2797        assert num.allclose(domain.quantities['stage'].explicit_update[1], 2.0/A)
     2798        assert num.allclose(domain.quantities['stage'].explicit_update[0], 2.0/A)
     2799        assert num.allclose(domain.quantities['stage'].explicit_update[2:], 0)       
     2800
    27972801
    27982802    def test_inflow_using_circle_function(self):
     
    28232827        # on a circle affecting triangles #0 and #1 (bac and bce)
    28242828        domain.forcing_terms = []
    2825         domain.forcing_terms.append(Inflow(domain, rate=lambda t: 2.,
    2826                                            center=(1,1), radius=1))
    2827 
    2828         domain.compute_forcing_terms()
    2829 
    2830         assert num.allclose(domain.quantities['stage'].explicit_update[1],
    2831                             2.0/pi)
    2832         assert num.allclose(domain.quantities['stage'].explicit_update[0],
    2833                             2.0/pi)
    2834         assert num.allclose(domain.quantities['stage'].explicit_update[2:], 0)
     2829        I = Inflow(domain, rate=lambda t: 2., center=(1,1), radius=1)
     2830        domain.forcing_terms.append(I)
     2831       
     2832        domain.compute_forcing_terms()       
     2833
     2834        A = I.exchange_area
     2835        assert num.allclose(A, 4) # Two triangles
     2836       
     2837        assert num.allclose(domain.quantities['stage'].explicit_update[1], 2.0/A)
     2838        assert num.allclose(domain.quantities['stage'].explicit_update[0], 2.0/A)
     2839        assert num.allclose(domain.quantities['stage'].explicit_update[2:], 0)       
     2840       
     2841
     2842
    28352843
    28362844    def test_inflow_catch_too_few_triangles(self):
     
    62336241        from anuga.geospatial_data.geospatial_data import Geospatial_data
    62346242
    6235         #--------------------------------------------------------------------
     6243       
     6244        # Get path where this test is run
     6245        path = get_pathname_from_package('anuga.shallow_water')       
     6246       
     6247       
     6248        #----------------------------------------------------------------------
    62366249        # Create domain
    62376250        #--------------------------------------------------------------------
     
    62766289        interior_regions = []
    62776290        for polygon in offending_regions:
    6278             interior_regions.append( [polygon, 100] )
    6279 
    6280         meshname = 'offending_mesh.msh'
     6291            interior_regions.append( [polygon, 100] ) 
     6292
     6293        meshname = os.path.join(path, 'offending_mesh.msh')
    62816294        create_mesh_from_regions(bounding_polygon,
    62826295                                 boundary_tags={'south': [0], 'east': [1],
     
    62906303        domain = Domain(meshname, use_cache=False, verbose=verbose)
    62916304
    6292         #--------------------------------------------------------------------
     6305        #----------------------------------------------------------------------
    62936306        # Fit data point to mesh
    6294         #--------------------------------------------------------------------
    6295         points_file = 'offending_point.pts'
    6296 
    6297         # This is the offending point
     6307        #----------------------------------------------------------------------
     6308
     6309        points_file = os.path.join(path, 'offending_point.pts')
     6310
     6311        # Offending point
    62986312        G = Geospatial_data(data_points=[[306953.344, 6194461.5]],
    62996313                            attributes=[1])
    63006314        G.export_points_file(points_file)
    6301 
    6302         domain.set_quantity('elevation', filename=points_file, use_cache=False,
    6303                             verbose=verbose, alpha=0.01)
     6315       
     6316        try:
     6317            domain.set_quantity('elevation', filename=points_file,
     6318                                use_cache=False, verbose=verbose, alpha=0.01)
     6319        except RuntimeError, e:
     6320            msg = 'Test failed: %s' % str(e)
     6321            raise Exception, msg
     6322            # clean up in case raise fails
     6323            os.remove(meshname)
     6324            os.remove(points_file)
     6325        else:
     6326            os.remove(meshname)
     6327            os.remove(points_file)           
     6328       
     6329        #finally:
     6330            # Cleanup regardless
     6331            #FIXME(Ole): Finally does not work like this in python2.4
     6332            #FIXME(Ole): Reinstate this when Python2.4 is out of the way
     6333            #FIXME(Ole): Python 2.6 apparently introduces something called 'with'           
     6334            #os.remove(meshname)
     6335            #os.remove(points_file)
     6336
     6337
     6338    def test_fitting_example_that_crashed_2(self):
     6339        """test_fitting_example_that_crashed_2
     6340       
     6341        This unit test has been derived from a real world example
     6342        (the JJKelly study, by Petar Milevski).
     6343       
     6344        It shows a condition where set_quantity crashes due to AtA
     6345        not being built properly
     6346       
     6347        See ticket:314
     6348        """
     6349
     6350        verbose = False       
     6351
     6352        from anuga.shallow_water import Domain
     6353        from anuga.pmesh.mesh_interface import create_mesh_from_regions
     6354        from anuga.geospatial_data import Geospatial_data
     6355       
     6356        # Get path where this test is run
     6357        path = get_pathname_from_package('anuga.shallow_water')       
     6358
     6359        meshname = os.path.join(path, 'test_mesh.msh')
     6360        W = 304180
     6361        S = 6185270
     6362        E = 307650
     6363        N = 6189040
     6364        maximum_triangle_area = 1000000
     6365
     6366        bounding_polygon = [[W, S], [E, S], [E, N], [W, N]]
     6367
     6368        create_mesh_from_regions(bounding_polygon,
     6369                                 boundary_tags={'south': [0],
     6370                                                'east': [1],
     6371                                                'north': [2],
     6372                                                'west': [3]},
     6373                                 maximum_triangle_area=maximum_triangle_area,
     6374                                 filename=meshname,
     6375                                 use_cache=False,
     6376                                 verbose=verbose)
     6377
     6378        domain = Domain(meshname, use_cache=True, verbose=verbose)
     6379       
     6380        # Large test set revealed one problem
     6381        points_file = os.path.join(path, 'test_points_large.csv')
     6382        try:
     6383            domain.set_quantity('elevation', filename=points_file,
     6384                                use_cache=False, verbose=verbose)
     6385        except AssertionError, e:
     6386            msg = 'Test failed: %s' % str(e)
     6387            raise Exception, msg
     6388            # Cleanup in case this failed
     6389            os.remove(meshname)
     6390
     6391        # Small test set revealed another problem
     6392        points_file = os.path.join(path, 'test_points_small.csv')
     6393        try:   
     6394            domain.set_quantity('elevation', filename=points_file,
     6395                                use_cache=False, verbose=verbose)                           
     6396        except AssertionError, e:
     6397            msg = 'Test failed: %s' % str(e)
     6398            raise Exception, msg
     6399            # Cleanup in case this failed
     6400            os.remove(meshname)
     6401        else:
     6402            os.remove(meshname)
     6403
     6404           
     6405    def test_total_volume(self):       
     6406        """test_total_volume
     6407       
     6408        Test that total volume can be computed correctly
     6409        """           
     6410
     6411        #----------------------------------------------------------------------
     6412        # Import necessary modules
     6413        #----------------------------------------------------------------------
     6414        from anuga.abstract_2d_finite_volumes.mesh_factory \
     6415                import rectangular_cross
     6416        from anuga.shallow_water import Domain
     6417
     6418        #----------------------------------------------------------------------
     6419        # Setup computational domain
     6420        #----------------------------------------------------------------------
     6421
     6422        length = 100.
     6423        width  = 20.
     6424        dx = dy = 5       # Resolution: of grid on both axes
     6425       
     6426        points, vertices, boundary = rectangular_cross(int(length/dx),
     6427                                                       int(width/dy),
     6428                                                       len1=length,
     6429                                                       len2=width)
     6430        domain = Domain(points, vertices, boundary)   
     6431
     6432        #----------------------------------------------------------------------
     6433        # Simple flat bottom bathtub
     6434        #----------------------------------------------------------------------
     6435
     6436        d = 1.0
     6437        domain.set_quantity('elevation', 0.0)
     6438        domain.set_quantity('stage', d)
     6439       
     6440        assert num.allclose(domain.compute_total_volume(), length*width*d)
     6441
     6442        #----------------------------------------------------------------------
     6443        # Slope
     6444        #----------------------------------------------------------------------
     6445               
     6446        slope = 1.0/10          # RHS drops to -10m
     6447        def topography(x, y):
     6448            return -x * slope
     6449
     6450        domain.set_quantity('elevation', topography)
     6451        domain.set_quantity('stage', 0.0)       # Domain full
     6452       
     6453        V = domain.compute_total_volume()
     6454        assert num.allclose(V, length*width*10/2)
     6455
     6456        domain.set_quantity('stage', -5.0)      # Domain 'half' full
     6457       
     6458        # IMPORTANT: Adjust stage to match elevation
     6459        domain.distribute_to_vertices_and_edges()
     6460       
     6461        V = domain.compute_total_volume()
     6462        assert num.allclose(V, width*(length/2)*5.0/2)
     6463
     6464
     6465    def test_volumetric_balance_computation(self):
     6466        """test_volumetric_balance_computation
     6467       
     6468        Test that total in and out flows are computed correctly
     6469        FIXME(Ole): This test is more about looking at the printed report
     6470        """
     6471
     6472        verbose = False
     6473
     6474        #----------------------------------------------------------------------
     6475        # Import necessary modules
     6476        #----------------------------------------------------------------------
     6477
     6478        from anuga.abstract_2d_finite_volumes.mesh_factory \
     6479                import rectangular_cross
     6480        from anuga.shallow_water import Domain
     6481        from anuga.shallow_water.shallow_water_domain import Reflective_boundary
     6482        from anuga.shallow_water.shallow_water_domain import Dirichlet_boundary
     6483        from anuga.shallow_water.shallow_water_domain import Inflow
     6484        from anuga.shallow_water.data_manager \
     6485                import get_flow_through_cross_section
     6486
     6487        #----------------------------------------------------------------------
     6488        # Setup computational domain
     6489        #----------------------------------------------------------------------
     6490
     6491        finaltime = 300.0
     6492        length = 300.
     6493        width  = 20.
     6494        dx = dy = 5       # Resolution: of grid on both axes
     6495       
     6496        # Input parameters
     6497        uh = 1.0
     6498        vh = 0.0
     6499        d = 1.0
     6500       
     6501        # 20 m^3/s in the x direction across entire domain
     6502        ref_flow = uh*d*width
     6503
     6504        points, vertices, boundary = rectangular_cross(int(length/dx),
     6505                                                       int(width/dy),
     6506                                                       len1=length,
     6507                                                       len2=width)
     6508
     6509        domain = Domain(points, vertices, boundary)   
     6510        domain.set_name('Inflow_flowline_test')              # Output name
     6511
     6512        #----------------------------------------------------------------------
     6513        # Setup initial conditions
     6514        #----------------------------------------------------------------------
     6515
     6516        slope = 0.0
     6517        def topography(x, y):
     6518            z=-x * slope
     6519            return z
     6520
     6521        # Use function for elevation
     6522        domain.set_quantity('elevation', topography)
     6523        domain.set_quantity('friction', 0.0)        # Constant friction
     6524               
     6525        domain.set_quantity('stage', expression='elevation')
     6526
     6527        #----------------------------------------------------------------------
     6528        # Setup boundary conditions
     6529        #----------------------------------------------------------------------
     6530
     6531        Br = Reflective_boundary(domain)      # Solid reflective wall
     6532               
     6533        # Constant flow into domain
     6534        # Depth = 1m, uh=1 m/s, i.e. a flow of 20 m^3/s
     6535       
     6536        Bi = Dirichlet_boundary([d, uh, vh])
     6537        Bo = Dirichlet_boundary([0, 0, 0])
     6538
     6539        domain.set_boundary({'left': Bi, 'right': Bo, 'top': Br, 'bottom': Br})
     6540
     6541        #----------------------------------------------------------------------
     6542        # Evolve system through time
     6543        #----------------------------------------------------------------------
     6544
     6545        for t in domain.evolve(yieldstep=100.0, finaltime=finaltime):
     6546            S = domain.volumetric_balance_statistics()
     6547            if verbose :
     6548                print domain.timestepping_statistics()
     6549                print S
     6550
     6551
     6552    def Xtest_inflow_using_flowline(self):
     6553        """test_inflow_using_flowline
     6554
     6555        Test the ability of a flowline to match inflow above the flowline by
     6556        creating constant inflow onto a circle at the head of a 20m wide by 300m
     6557        long plane dipping at various slopes with a perpendicular flowline and
     6558        gauge downstream of the inflow and a 45 degree flowlines at 200m
     6559        downstream.
     6560        """
     6561
     6562        # FIXME(Ole): Move this and the following test to validate_all.py as
     6563        # they are rather time consuming
     6564       
     6565        verbose = True
     6566
     6567        #----------------------------------------------------------------------
     6568        # Import necessary modules
     6569        #----------------------------------------------------------------------
     6570
     6571        from anuga.abstract_2d_finite_volumes.mesh_factory \
     6572                import rectangular_cross
     6573        from anuga.shallow_water import Domain
     6574        from anuga.shallow_water.shallow_water_domain import Reflective_boundary
     6575        from anuga.shallow_water.shallow_water_domain import Dirichlet_boundary
     6576        from anuga.shallow_water.shallow_water_domain import Inflow
     6577        from anuga.shallow_water.data_manager \
     6578                import get_flow_through_cross_section
     6579        from anuga.abstract_2d_finite_volumes.util \
     6580                import sww2csv_gauges, csv2timeseries_graphs
     6581
     6582        #----------------------------------------------------------------------
     6583        # Setup computational domain
     6584        #----------------------------------------------------------------------
     6585
     6586        number_of_inflows = 2       # Number of inflows on top of each other
     6587        finaltime = 1000.0          # 6000.0
     6588
     6589        length = 300.
     6590        width  = 20.
     6591        dx = dy = 5                 # Resolution: of grid on both axes
     6592
     6593        points, vertices, boundary = rectangular_cross(int(length/dx),
     6594                                                       int(width/dy),
     6595                                                       len1=length,
     6596                                                       len2=width)
     6597
     6598        for mannings_n in [0.150, 0.07, 0.035]:    #[0.012, 0.035, 0.070, 0.150]:
     6599            for slope in [1.0/300, 1.0/150, 1.0/75]:
     6600                # Loop over a range of bedslopes representing sub to
     6601                # super critical flows
     6602                if verbose:
     6603                    print
     6604                    print 'Slope:', slope, 'Mannings n:', mannings_n
     6605                domain = Domain(points, vertices, boundary)   
     6606                domain.set_name('Inflow_flowline_test')     # Output name
     6607
     6608                #--------------------------------------------------------------
     6609                # Setup initial conditions
     6610                #--------------------------------------------------------------
     6611
     6612                def topography(x, y):
     6613                    z = -x * slope
     6614                    return z
     6615
     6616                # Use function for elevation
     6617                domain.set_quantity('elevation', topography)
     6618                # Constant friction of conc surface   
     6619                domain.set_quantity('friction', mannings_n)
     6620                # Dry initial condition
     6621                domain.set_quantity('stage', expression='elevation')
     6622
     6623                #--------------------------------------------------------------
     6624                # Setup boundary conditions
     6625                #--------------------------------------------------------------
     6626
     6627                Br = Reflective_boundary(domain)      # Solid reflective wall
     6628                # Outflow stage at -0.9m d=0.1m
     6629                Bo = Dirichlet_boundary([-100, 0, 0])
     6630
     6631                domain.set_boundary({'left': Br, 'right': Bo,
     6632                                     'top': Br,  'bottom': Br})
     6633
     6634                #--------------------------------------------------------------
     6635                # Setup Inflow
     6636                #--------------------------------------------------------------
     6637
     6638                # Fixed Flowrate onto Area
     6639                fixed_inflow = Inflow(domain, center=(10.0, 10.0),
     6640                                      radius=5.00, rate=10.00)   
     6641
     6642                # Stack this flow
     6643                for i in range(number_of_inflows):
     6644                    domain.forcing_terms.append(fixed_inflow)
     6645               
     6646                ref_flow = fixed_inflow.rate*number_of_inflows
     6647
     6648                #--------------------------------------------------------------
     6649                # Evolve system through time
     6650                #--------------------------------------------------------------
     6651
     6652                for t in domain.evolve(yieldstep=100.0, finaltime=finaltime):
     6653                    pass
     6654                    #if verbose :
     6655                    #    print domain.timestepping_statistics()
     6656
     6657                if verbose:
     6658                    print domain.volumetric_balance_statistics()                                                           
     6659                #--------------------------------------------------------------
     6660                # Compute flow thru flowlines ds of inflow
     6661                #--------------------------------------------------------------
     6662                   
     6663                # Square on flowline at 200m
     6664                q = domain.get_flow_through_cross_section([[200.0,  0.0],
     6665                                                           [200.0, 20.0]])
     6666                if verbose:
     6667                    print ('90 degree flowline: ANUGA = %f, Ref = %f'
     6668                           % (q, ref_flow))
     6669
     6670                msg = ('Predicted flow was %f, should have been %f'
     6671                       % (q, ref_flow))
     6672                assert num.allclose(q, ref_flow, rtol=1.0e-2), msg         
     6673
     6674                           
     6675                # 45 degree flowline at 200m
     6676                q = domain.get_flow_through_cross_section([[200.0,  0.0],
     6677                                                           [220.0, 20.0]])
     6678                if verbose:
     6679                    print ('45 degree flowline: ANUGA = %f, Ref = %f'
     6680                           % (q, ref_flow))
     6681                   
     6682                msg = ('Predicted flow was %f, should have been %f'
     6683                       % (q, ref_flow))
     6684                assert num.allclose(q, ref_flow, rtol=1.0e-2), msg         
     6685
     6686                # Stage recorder (gauge) in middle of plane at 200m
     6687                x = 200.0
     6688                y = 10.00
     6689                w = domain.get_quantity('stage').\
     6690                        get_values(interpolation_points=[[x, y]])[0]
     6691                z = domain.get_quantity('elevation').\
     6692                        get_values(interpolation_points=[[x, y]])[0]
     6693                domain_depth = w-z
     6694
     6695                # Compute normal depth at gauge location using Manning equation
     6696                # v=1/n*(r^2/3)*(s^0.5) or r=(Q*n/(s^0.5*W))^0.6
     6697                normal_depth = (ref_flow*mannings_n/(slope**0.5*width))**0.6
     6698                msg = ('Predicted depth of flow was %f, should have been %f'
     6699                       % (domain_depth, normal_depth))
     6700                if verbose:
     6701                    diff = abs(domain_depth-normal_depth)
     6702                    print ('Depth: ANUGA = %f, Mannings = %f, E=%f'
     6703                           % (domain_depth, normal_depth, diff/domain_depth))
     6704                    assert diff < 0.1
     6705                   
     6706                if slope >= 1.0/100:
     6707                    # Really super critical flow is not as stable.
     6708                    #assert num.allclose(domain_depth, normal_depth,
     6709                    #                    rtol=1.0e-1), msg
     6710                    pass
     6711                else:
     6712                    pass
     6713                    #assert num.allclose(domain_depth, normal_depth,
     6714                    #                    rtol=1.0e-2), msg
     6715
     6716
     6717    def Xtest_friction_dependent_flow_using_flowline(self):
     6718        """test_friction_dependent_flow_using_flowline
     6719       
     6720        Test the internal flow (using flowline) as a function of
     6721        different values of Mannings n and different slopes.
     6722       
     6723        Flow is applied in the form of boundary conditions with fixed momentum.
     6724        """
     6725
     6726        verbose = True
     6727
     6728        #----------------------------------------------------------------------
     6729        # Import necessary modules
     6730        #----------------------------------------------------------------------
     6731
     6732        from anuga.abstract_2d_finite_volumes.mesh_factory \
     6733                import rectangular_cross
     6734        from anuga.shallow_water import Domain
     6735        from anuga.shallow_water.shallow_water_domain import Reflective_boundary
     6736        from anuga.shallow_water.shallow_water_domain import Dirichlet_boundary
     6737        from anuga.shallow_water.shallow_water_domain import Inflow
     6738        from anuga.shallow_water.data_manager \
     6739                import get_flow_through_cross_section
     6740        from anuga.abstract_2d_finite_volumes.util \
     6741                import sww2csv_gauges, csv2timeseries_graphs
     6742
     6743
     6744        #----------------------------------------------------------------------
     6745        # Setup computational domain
     6746        #----------------------------------------------------------------------
     6747
     6748        finaltime = 1000.0
     6749
     6750        length = 300.
     6751        width  = 20.
     6752        dx = dy = 5       # Resolution: of grid on both axes
     6753       
     6754        # Input parameters
     6755        uh = 1.0
     6756        vh = 0.0
     6757        d = 1.0
     6758       
     6759        ref_flow = uh*d*width # 20 m^3/s in the x direction across entire domain
     6760
     6761        points, vertices, boundary = rectangular_cross(int(length/dx),
     6762                                                       int(width/dy),
     6763                                                       len1=length,
     6764                                                       len2=width)
     6765
     6766        for mannings_n in [0.035]:          #[0.0, 0.012, 0.035]:
     6767            for slope in [1.0/300]:         #[0.0, 1.0/300, 1.0/150]:
     6768                # Loop over a range of bedslopes representing
     6769                # sub to super critical flows
     6770                if verbose:
     6771                    print
     6772                    print 'Slope:', slope, 'Mannings n:', mannings_n
     6773                domain = Domain(points, vertices, boundary)   
     6774                domain.set_name('Inflow_flowline_test')     # Output name
     6775
     6776                #--------------------------------------------------------------
     6777                # Setup initial conditions
     6778                #--------------------------------------------------------------
     6779
     6780                def topography(x, y):
     6781                    z = -x * slope
     6782                    return z
     6783
     6784                # Use function for elevation
     6785                domain.set_quantity('elevation', topography)
     6786                # Constant friction
     6787                domain.set_quantity('friction', mannings_n)
     6788               
     6789                #domain.set_quantity('stage', expression='elevation')
     6790                     
     6791                # Set initial flow as depth=1m, uh=1.0 m/s, vh = 0.0
     6792                # making it 20 m^3/s across entire domain
     6793                domain.set_quantity('stage', expression='elevation + %f' % d)
     6794                domain.set_quantity('xmomentum', uh)
     6795                domain.set_quantity('ymomentum', vh)               
     6796
     6797                #--------------------------------------------------------------
     6798                # Setup boundary conditions
     6799                #--------------------------------------------------------------
     6800
     6801                Br = Reflective_boundary(domain)      # Solid reflective wall
     6802               
     6803                # Constant flow in and out of domain
     6804                # Depth = 1m, uh=1 m/s, i.e. a flow of 20 m^3/s
     6805                # across boundaries
     6806                Bi = Dirichlet_boundary([d, uh, vh])
     6807                Bo = Dirichlet_boundary([-length*slope+d, uh, vh])
     6808                #Bo = Dirichlet_boundary([-100, 0, 0])
     6809
     6810                domain.set_boundary({'left': Bi, 'right': Bo,
     6811                                     'top': Br,  'bottom': Br})
     6812
     6813                #--------------------------------------------------------------
     6814                # Evolve system through time
     6815                #--------------------------------------------------------------
     6816
     6817                for t in domain.evolve(yieldstep=100.0, finaltime=finaltime):
     6818                    if verbose :
     6819                        print domain.timestepping_statistics()
     6820                        print domain.volumetric_balance_statistics()
     6821
     6822                # 90 degree flowline at 200m
     6823                q = domain.get_flow_through_cross_section([[200.0,  0.0],
     6824                                                           [200.0, 20.0]])
     6825                msg = ('Predicted flow was %f, should have been %f'
     6826                       % (q, ref_flow))
     6827                if verbose:
     6828                    print ('90 degree flowline: ANUGA = %f, Ref = %f'
     6829                           % (q, ref_flow))
     6830
     6831                # 45 degree flowline at 200m
     6832                q = domain.get_flow_through_cross_section([[200.0,  0.0],
     6833                                                           [220.0, 20.0]])
     6834                msg = ('Predicted flow was %f, should have been %f'
     6835                       % (q, ref_flow))
     6836                if verbose:
     6837                    print ('45 degree flowline: ANUGA = %f, Ref = %f'
     6838                           % (q, ref_flow))
     6839
     6840                # Stage recorder (gauge) in middle of plane at 200m
     6841                x = 200.0
     6842                y = 10.00               
     6843                w = domain.get_quantity('stage').\
     6844                        get_values(interpolation_points=[[x, y]])[0]
     6845                z = domain.get_quantity('elevation').\
     6846                        get_values(interpolation_points=[[x, y]])[0]
     6847                domain_depth = w-z
     6848
     6849                xmom = domain.get_quantity('xmomentum').\
     6850                        get_values(interpolation_points=[[x, y]])[0]
     6851                ymom = domain.get_quantity('ymomentum').\
     6852                        get_values(interpolation_points=[[x, y]])[0]           
     6853                if verbose:                   
     6854                    print ('At interpolation point (h, uh, vh): ',
     6855                           domain_depth, xmom, ymom)
     6856                    print 'uh * d * width = ', xmom*domain_depth*width
     6857                               
     6858                if slope > 0.0:
     6859                    # Compute normal depth at gauge location using Manning eqn
     6860                    # v=1/n*(r^2/3)*(s^0.5) or r=(Q*n/(s^0.5*W))^0.6
     6861                    normal_depth = (ref_flow*mannings_n/(slope**0.5*width))**0.6
     6862                    if verbose:
     6863                        print ('Depth: ANUGA = %f, Mannings = %f'
     6864                               % (domain_depth, normal_depth))
     6865
     6866#################################################################################
    63046867
    63056868if __name__ == "__main__":
  • branches/numpy/anuga/utilities/log.py

    r6553 r6689  
    77Use it this way:
    88    import anuga.utilities.log as log
     9    log.console_logging_level = log.DEBUG
    910    log.debug('A message at DEBUG level')
    1011    log.info('Another message, INFO level')
     
    1718Until the first call to log() the user is free to play with the module data
    1819to configure the logging.
     20Note that this module uses features of the logging package that were introduced
     21in python2.5.  If running on earlier versions, these features are disables:
     22    . Module name + line number
    1923'''
    2024
    2125import sys
    22 import os.path
     26import os
     27import sys
    2328import traceback
    2429import logging
     
    2732################################################################################
    2833# Module variables - only one copy of these, ever.
     34#
     35# The console logging level is set to a high level, like CRITICAL.  The logfile
     36# logging is set lower, between DEBUG and CRITICAL.  The idea is to log least to
     37# the console, but ensure that everything that goes to the console *will* also
     38# appear in the log file.  There is code to ensure log <= console levels.
    2939################################################################################
    3040
     
    3343
    3444# logging level for the console
    35 console_logging_level = logging.INFO
     45console_logging_level = logging.CRITICAL
    3646
    3747# logging level for the logfile
    38 log_logging_level = logging.DEBUG
    39 
    40 # The name of the file to log to.
    41 log_filename = './anuga.log'
     48log_logging_level = logging.INFO
     49
     50# The default name of the file to log to.
     51log_filename = os.path.join('.', 'anuga.log')
     52
     53# set module variables so users don't have to do 'import logging'.
     54CRITICAL = logging.CRITICAL
     55ERROR = logging.ERROR
     56WARNING = logging.WARNING
     57INFO = logging.INFO
     58DEBUG = logging.DEBUG
     59NOTSET = logging.NOTSET
    4260
    4361
     
    5876    '''
    5977
    60     global _setup
     78    global _setup, log_logging_level
     79
     80    # get running python version for later
     81    (version_major, version_minor, _, _, _) = sys.version_info
    6182
    6283    # have we been setup?
    6384    if not _setup:
     85        # sanity check the logging levels, require console >= file
     86        if console_logging_level < log_logging_level:
     87            log_logging_level = console_logging_level
     88
    6489        # setup the file logging system
    65         fmt = '%(asctime)s %(levelname)-8s %(mname)25s:%(lnum)-4d|%(message)s'
     90        if version_major >= 2 and version_minor >= 5:
     91            fmt = '%(asctime)s %(levelname)-8s %(mname)25s:%(lnum)-4d|%(message)s'
     92        else:
     93            fmt = '%(asctime)s %(levelname)-8s|%(message)s'
    6694        logging.basicConfig(level=log_logging_level, format=fmt,
    6795                            filename=log_filename, filemode='w')
     
    80108                        logging.getLevelName(log_logging_level),
    81109                        logging.getLevelName(console_logging_level)))
    82         logging.log(logging.CRITICAL, start_msg,
    83                     extra={'mname': __name__, 'lnum': 0})
     110        if version_major >= 2 and version_minor >= 5:
     111            logging.log(logging.CRITICAL, start_msg,
     112                        extra={'mname': __name__, 'lnum': 0})
     113        else:
     114            logging.log(logging.CRITICAL, start_msg)
    84115
    85116        # mark module as *setup*
     
    89120    frames = traceback.extract_stack()
    90121    frames.reverse()
    91     for (mname, lnum, _, _) in frames:
    92         mname = os.path.basename(mname).rsplit('.', 1)[0]
     122    for (fpath, lnum, mname, _) in frames:
     123        fname = os.path.basename(mname).rsplit('.', 1)[0]
    93124        if mname != __name__:
    94125            break
    95126
    96     logging.log(level, msg, extra={'mname': mname, 'lnum': lnum})
     127    if version_major >= 2 and version_minor >= 5:
     128        logging.log(level, msg, extra={'mname': fname, 'lnum': lnum})
     129    else:
     130        logging.log(level, msg)
    97131
    98132################################################################################
     
    130164    log(logging.CRITICAL, msg)
    131165
     166def resource_usage(level=logging.CRITICAL):
     167    '''Log resource usage at given log level.'''
     168
     169    if sys.platform != 'win32':
     170        _proc_status = '/proc/%d/status' % os.getpid()
     171        _scale = {'KB': 1024.0, 'MB': 1024.0*1024.0, 'GB': 1024.0*1024.0*1024.0,
     172                  'kB': 1024.0, 'mB': 1024.0*1024.0, 'gB': 1024.0*1024.0*1024.0}
     173
     174        def _VmB(VmKey):
     175            '''Get number of virtual bytes used.'''
     176
     177            # get pseudo file /proc/<pid>/status
     178            try:
     179                t = open(_proc_status)
     180                v = t.read()
     181                t.close()
     182            except IOError:
     183                return 0.0
     184
     185            # get VmKey line, eg: 'VmRSS: 999 kB\n ...
     186            i = v.index(VmKey)
     187            v = v[i:].split(None, 3)
     188            if len(v) < 3:
     189                return 0.0
     190
     191            # convert Vm value to bytes
     192            return float(v[1]) * _scale[v[2]]
     193
     194        def memory(since=0.0):
     195            '''Get virtual memory usage in bytes.'''
     196
     197            return _VmB('VmSize:') - since
     198
     199        def resident(since=0.0):
     200            '''Get resident memory usage in bytes.'''
     201
     202            return _VmB('VmRSS:') - since
     203
     204        def stacksize(since=0.0):
     205            '''Get stack size in bytes.'''
     206
     207            return _VmB('VmStk:') - since
     208
     209        msg = ('Resource usage: memory=%.1f resident=%.1f stacksize=%.1f'
     210               % (memory()/_scale['GB'], resident()/_scale['GB'],
     211                  stacksize()/_scale['GB']))
     212        log(level, msg)
     213    else:
     214        pass
  • branches/numpy/anuga/utilities/polygon.py

    r6553 r6689  
    10701070    class Found(exceptions.Exception): pass
    10711071
     1072    polygon = ensure_numeric(polygon)
     1073   
    10721074    point_in = False
    10731075    while not point_in:
     
    11191121    # TO DO check if any of the regions fall inside one another
    11201122
    1121     print '--------------------------------------------------------------------'
    1122     print ('Polygon   Max triangle area (m^2)   Total area (km^2)   '
    1123            'Estimated #triangles')
    1124     print '--------------------------------------------------------------------'
    1125 
     1123    print '----------------------------------------------------------------------------'
     1124    print 'Polygon   Max triangle area (m^2)   Total area (km^2)   Estimated #triangles'
     1125    print '----------------------------------------------------------------------------'   
     1126       
    11261127    no_triangles = 0.0
    11271128    area = polygon_area(bounding_poly)
  • branches/numpy/anuga/utilities/system_tools.py

    r6428 r6689  
    77import os
    88import string
     9import urllib
     10import urllib2
     11import getpass
     12import tarfile
     13import md5
     14
    915
    1016def log_to_file(filename, s, verbose=False):
     
    4955    """
    5056
    51     # Create dummy info
    52     #info = 'Revision: Version info could not be obtained.'
    53     #info += 'A command line version of svn must be availbable '
    54     #info += 'on the system PATH, access to the subversion '
    55     #info += 'repository is necessary and the output must '
    56     #info += 'contain a line starting with "Revision:"'
    57    
    58 
    59     #FIXME (Ole): Change this so that svn info is attempted first.
    60     # If that fails, try to read a stored file with that same info (this would be created by e.g. the release script). Failing that, throw an exception.
    61 
    62     #FIXME (Ole): Move this and store_version_info to utilities
    63 
    64 
     57    # FIXME (Ole): Change this so that svn info is attempted first.
     58    # If that fails, try to read a stored file with that same info (this would
     59    # be created by e.g. the release script). Failing that, throw an exception.
     60    # FIXME (Ole): Move this and store_version_info to utilities
    6561    try:
    6662        from anuga.stored_version_info import version_info
    6763    except:
    68         msg = 'No version info stored and command "svn" is not '
    69         msg += 'recognised on the system PATH.\n\n'
    70         msg += 'If ANUGA has been installed from a distribution e.g. as '
    71         msg += 'obtained from SourceForge,\n'
    72         msg += 'the version info should be '
    73         msg += 'available in the automatically generated file '
    74         msg += 'stored_version_info.py\n'
    75         msg += 'in the anuga root directory.\n'
    76         msg += 'If run from a Subversion sandpit, '
    77         msg += 'ANUGA will try to obtain the version info '
    78         msg += 'by using the command: "svn info".\n'
    79         msg += 'In this case, make sure svn is accessible on the system path. '
    80         msg += 'Simply aliasing svn to the binary will not work. '
    81         msg += 'Good luck!'
     64        msg = '''
     65No version info stored and command 'svn' is not recognised on the system PATH.
     66
     67If ANUGA has been installed from a distribution e.g. as obtained from
     68SourceForge, the version info should be available in the automatically
     69generated file 'stored_version_info.py' in the anuga root directory.
     70
     71If run from a Subversion sandpit, ANUGA will try to obtain the version info
     72by using the command 'svn info'.  In this case, make sure 'svn' is accessible
     73on the system PATH.  Simply aliasing 'svn' to the binary will not work.
     74
     75If you are using Windows, you have to install the file svn.exe which can be
     76obtained from http://www.collab.net/downloads/subversion.
     77
     78Good luck!
     79'''
    8280
    8381        # No file available - try using Subversion
     
    9189            else:
    9290                fid = os.popen('svn info 2>/dev/null')
    93        
    9491        except:
    9592            raise Exception(msg)
    9693        else:
    97             #print 'Got version from svn'           
    9894            version_info = fid.read()
    99            
    100             if version_info == '':
    101                 raise Exception(msg)   
     95
     96        if version_info == '':
     97            raise Exception, msg
    10298    else:
    10399        pass
    104100        #print 'Got version from file'
    105 
    106101           
    107102    for line in version_info.split('\n'):
     
    110105
    111106    fields = line.split(':')
    112     msg = 'Keyword "Revision" was not found anywhere in text: %s' %version_info
     107    msg = 'Keyword "Revision" was not found anywhere in text: %s' % version_info
    113108    assert fields[0].startswith('Revision'), msg           
    114109
     
    118113        msg = 'Revision number must be an integer. I got %s' %fields[1]
    119114        msg += 'Check that the command svn is on the system path'
    120         raise Exception(msg)               
    121        
     115        raise Exception, msg
     116
    122117    return revision_number
    123118
     
    181176        if verbose is True:
    182177            print 'Version info stored to %s' %filename
     178
    183179
    184180def safe_crc(string):
     
    319315
    320316
     317##
     318# @brief Get list of variable names in a python expression string.
     319# @param source A string containing a python expression.
     320# @return A list of variable name strings.
     321# @note Throws SyntaxError exception if not a valid expression.
     322def get_vars_in_expression(source):
     323    '''Get list of variable names in a python expression.'''
     324
     325    import compiler
     326    from compiler.ast import Node
     327
     328    ##
     329    # @brief Internal recursive function.
     330    # @param node An AST parse Node.
     331    # @param var_list Input list of variables.
     332    # @return An updated list of variables.
     333    def get_vars_body(node, var_list=[]):
     334        if isinstance(node, Node):
     335            if node.__class__.__name__ == 'Name':
     336                for child in node.getChildren():
     337                    if child not in var_list:
     338                        var_list.append(child)
     339            for child in node.getChildren():
     340                if isinstance(child, Node):
     341                    for child in node.getChildren():
     342                        var_list = get_vars_body(child, var_list)
     343                    break
     344
     345        return var_list
     346
     347    return get_vars_body(compiler.parse(source))
     348
     349
     350##
     351# @brief Get a file from the web.
     352# @param file_url URL of the file to fetch.
     353# @param file_name Path to file to create in the filesystem.
     354# @param auth Auth tuple (httpproxy, proxyuser, proxypass).
     355# @param blocksize Read file in this block size.
     356# @return 'auth' tuple for subsequent calls, if successful.
     357# @note If 'auth' not supplied, will prompt user.
     358# @note Will try using environment variable HTTP_PROXY for proxy server.
     359# @note Will try using environment variable PROXY_USERNAME for proxy username.
     360# @note Will try using environment variable PROXY_PASSWORD for proxy password.
     361def get_web_file(file_url, file_name, auth=None, blocksize=1024*1024):
     362    '''Get a file from the web.
     363
     364    file_url:  The URL of the file to get
     365    file_name: Local path to save loaded file in
     366    auth:      A tuple (httpproxy, proxyuser, proxypass)
     367    blocksize: Block size of file reads
     368   
     369    Will try simple load through urllib first.  Drop down to urllib2
     370    if there is a proxy and it requires authentication.
     371
     372    Environment variable HTTP_PROXY can be used to supply proxy information.
     373    PROXY_USERNAME is used to supply the authentication username.
     374    PROXY_PASSWORD supplies the password, if you dare!
     375    '''
     376
     377    # Simple fetch, if fails, check for proxy error
     378    try:
     379        urllib.urlretrieve(file_url, file_name)
     380        return None     # no proxy, no auth required
     381    except IOError, e:
     382        if e[1] != 407:
     383            raise       # raise error if *not* proxy auth error
     384
     385    # We get here if there was a proxy error, get file through the proxy
     386    # unpack auth info
     387    try:
     388        (httpproxy, proxyuser, proxypass) = auth
     389    except:
     390        (httpproxy, proxyuser, proxypass) = (None, None, None)
     391
     392    # fill in any gaps from the environment
     393    if httpproxy is None:
     394        httpproxy = os.getenv('HTTP_PROXY')
     395    if proxyuser is None:
     396        proxyuser = os.getenv('PROXY_USERNAME')
     397    if proxypass is None:
     398        proxypass = os.getenv('PROXY_PASSWORD')
     399
     400    # Get auth info from user if still not supplied
     401    if httpproxy is None or proxyuser is None or proxypass is None:
     402        print '-'*80
     403        print ('You need to supply proxy authentication information.  '
     404               'Use environment variables\n'
     405               'HTTP_PROXY, PROXY_USERNAME and PROXY_PASSWORD to bypass '
     406               'entry here:')
     407        if httpproxy is None:
     408            httpproxy = raw_input('  proxy server: ')
     409        if proxyuser is None:
     410            proxyuser = raw_input('proxy username: ')
     411        if proxypass is None:
     412            proxypass = getpass.getpass('proxy password: ')
     413        print '-'*80
     414
     415    # the proxy URL cannot start with 'http://', we add that later
     416    httpproxy = httpproxy.lower()
     417    if httpproxy.startswith('http://'):
     418        httpproxy = httpproxy.replace('http://', '', 1)
     419
     420    # open remote file
     421    proxy = urllib2.ProxyHandler({'http': 'http://' + proxyuser
     422                                              + ':' + proxypass
     423                                              + '@' + httpproxy})
     424    authinfo = urllib2.HTTPBasicAuthHandler()
     425    opener = urllib2.build_opener(proxy, authinfo, urllib2.HTTPHandler)
     426    urllib2.install_opener(opener)
     427    webget = urllib2.urlopen(file_url)
     428
     429    # transfer file to local filesystem
     430    fd = open(file_name, 'wb')
     431    while True:
     432        data = webget.read(blocksize)
     433        if len(data) == 0:
     434            break
     435        fd.write(data)
     436    fd.close
     437    webget.close()
     438
     439    # return successful auth info
     440    return (httpproxy, proxyuser, proxypass)
     441
     442
     443##
     444# @brief Tar a file (or directory) into a tarfile.
     445# @param files A list of files (or directories) to tar.
     446# @param tarfile The created tarfile name.
     447# @note We use gzip compression.
     448def tar_file(files, tarname):
     449    '''Compress a file or directory into a tar file.'''
     450
     451    o = tarfile.open(tarname, 'w:gz')
     452    for file in files:
     453        o.add(file)
     454    o.close()
     455
     456
     457##
     458# @brief Untar a file into an optional target directory.
     459# @param tarname Name of the file to untar.
     460# @param target_dir Directory to untar into.
     461def untar_file(tarname, target_dir='.'):
     462    '''Uncompress a tar file.'''
     463
     464    o = tarfile.open(tarname, 'r:gz')
     465    members = o.getmembers()
     466    for member in members:
     467        o.extract(member, target_dir)
     468    o.close()
     469
     470
     471##
     472# @brief Return a hex digest string for a given file.
     473# @param filename Path to the file of interest.
     474# @param blocksize Size of data blocks to read.
     475# @return A hex digest string (16 bytes).
     476# @note Uses MD5 digest.
     477def get_file_hexdigest(filename, blocksize=1024*1024*10):
     478    '''Get a hex digest of a file.'''
     479   
     480    m = md5.new()
     481    fd = open(filename, 'r')
     482           
     483    while True:
     484        data = fd.read(blocksize)
     485        if len(data) == 0:
     486            break
     487        m.update(data)
     488                                                               
     489    fd.close()
     490    return m.hexdigest()
     491
     492    fd = open(filename, 'r')
     493
     494
     495##
     496# @brief Create a file containing a hexdigest string of a data file.
     497# @param data_file Path to the file to get the hexdigest from.
     498# @param digest_file Path to hexdigest file to create.
     499# @note Uses MD5 digest.
     500def make_digest_file(data_file, digest_file):
     501    '''Create a file containing the hex digest string of a data file.'''
     502   
     503    hexdigest = get_file_hexdigest(data_file)
     504    fd = open(digest_file, 'w')
     505    fd.write(hexdigest)
     506    fd.close()
     507
     508
     509##
     510# @brief Function to return the length of a file.
     511# @param in_file Path to file to get length of.
     512# @return Number of lines in file.
     513# @note Doesn't count '\n' characters.
     514# @note Zero byte file, returns 0.
     515# @note No \n in file at all, but >0 chars, returns 1.
     516def file_length(in_file):
     517    '''Function to return the length of a file.'''
     518
     519    fid = open(in_file)
     520    data = fid.readlines()
     521    fid.close()
     522    return len(data)
     523
     524
  • branches/numpy/anuga/utilities/test_polygon.py

    r6553 r6689  
    13811381            assert status == 1
    13821382
    1383             msg = ('Orientation of line shouldn not matter.\n'
     1383            msg = ('Orientation of line should not matter.\n'
    13841384                   'However, segment [%f,%f], [%f, %f]' %
    13851385                   (x, 0, common_end_point[0], common_end_point[1]))
  • branches/numpy/anuga/utilities/test_system_tools.py

    r6428 r6689  
    44import unittest
    55import numpy as num
     6import random
     7import tempfile
    68import zlib
    79import os
     
    8587            pass
    8688        else:
    87             test_array = num.array([[7.0, 3.14], [-31.333, 0.0]],
    88                                    dtype=num.float)
     89            test_array = num.array([[7.0, 3.14], [-31.333, 0.0]])
    8990
    9091            # First file
     
    305306        os.remove(FILENAME)
    306307
    307 #        print ('test_string_to_char: str_list=%s, new_str_list=%s'
    308                 #               % (str(str_list), str(new_str_list)))
     308
     309    def test_get_vars_in_expression(self):
     310        '''Test the 'get vars from expression' code.'''
     311
     312        def test_it(source, expected):
     313            result = get_vars_in_expression(source)
     314            result.sort()
     315            expected.sort()
     316            msg = ("Source: '%s'\nResult: %s\nExpected: %s"
     317                   % (source, str(result), str(expected)))
     318            self.failUnlessEqual(result, expected, msg)
     319               
     320        source = 'fred'
     321        expected = ['fred']
     322        test_it(source, expected)
     323
     324        source = 'tom + dick'
     325        expected = ['tom', 'dick']
     326        test_it(source, expected)
     327
     328        source = 'tom * (dick + harry)'
     329        expected = ['tom', 'dick', 'harry']
     330        test_it(source, expected)
     331
     332        source = 'tom + dick**0.5 / (harry - tom)'
     333        expected = ['tom', 'dick', 'harry']
     334        test_it(source, expected)
     335
     336
     337    def test_tar_untar_files(self):
     338        '''Test that tarring & untarring files is OK.'''
     339
     340        num_lines = 100
     341        line_size = 100
     342
     343        # these test files must exist in the current directory
     344        # create them with random data
     345        files = ('alpha', 'beta', 'gamma')
     346        for file in files:
     347            fd = open(file, 'w')
     348            line = ''
     349            for i in range(num_lines):
     350                for j in range(line_size):
     351                    line += chr(random.randint(ord('A'), ord('Z')))
     352                line += '\n'
     353                fd.write(line)
     354            fd.close()
     355
     356        # name of tar file and test (temp) directory
     357        tar_filename = 'test.tgz'
     358        tmp_dir = tempfile.mkdtemp()
     359
     360        # tar and untar the test files into a temporary directory
     361        tar_file(files, tar_filename)
     362        untar_file(tar_filename, tmp_dir)
     363
     364        # see if original files and untarred ones are the same
     365        for file in files:
     366            fd = open(file, 'r')
     367            orig = fd.readlines()
     368            fd.close()
     369
     370            fd = open(os.path.join(tmp_dir, file), 'r')
     371            copy = fd.readlines()
     372            fd.close()
     373
     374            msg = "Original file %s isn't the same as untarred copy?" % file
     375            self.failUnless(orig == copy, msg)
     376
     377        # clean up
     378        for file in files:
     379            os.remove(file)
     380        os.remove(tar_filename)
     381
     382
     383    def test_file_digest(self):
     384        '''Test that file digest functions give 'correct' answer.
     385       
     386        Not a good test as we get 'expected_digest' from a digest file,
     387        but *does* alert us if the digest algorithm ever changes.
     388        '''
     389
     390        # we expect this digest string from the data file
     391        expected_digest = '831a1dde6edd365ec4163a47871fa21b'
     392
     393        # prepare test directory and filenames
     394        tmp_dir = tempfile.mkdtemp()
     395        data_file = os.path.join(tmp_dir, 'test.data')
     396        digest_file = os.path.join(tmp_dir, 'test.digest')
     397
     398        # create the data file
     399        data_line = 'The quick brown fox jumps over the lazy dog. 0123456789\n'
     400        fd = open(data_file, 'w')
     401        for line in range(100):
     402            fd.write(data_line)
     403        fd.close()
     404
     405        # create the digest file
     406        make_digest_file(data_file, digest_file)
     407
     408        # get digest string for the data file
     409        digest = get_file_hexdigest(data_file)
     410
     411        # check that digest is as expected, string
     412        msg = ("Digest string wrong, got '%s', expected '%s'"
     413               % (digest, expected_digest))
     414        self.failUnless(expected_digest == digest, msg)
     415
     416        # check that digest is as expected, file
     417        msg = ("Digest file wrong, got '%s', expected '%s'"
     418               % (digest, expected_digest))
     419        fd = open(digest_file, 'r')
     420        digest = fd.readline()
     421        fd.close()
     422        self.failUnless(expected_digest == digest, msg)
     423
     424
     425    def test_file_length_function(self):
     426        '''Test that file_length() give 'correct' answer.'''
     427
     428        # prepare test directory and filenames
     429        tmp_dir = tempfile.mkdtemp()
     430        test_file1 = os.path.join(tmp_dir, 'test.file1')
     431        test_file2 = os.path.join(tmp_dir, 'test.file2')
     432        test_file3 = os.path.join(tmp_dir, 'test.file3')
     433        test_file4 = os.path.join(tmp_dir, 'test.file4')
     434
     435        # create files of known length
     436        fd = open(test_file1, 'w')      # 0 lines
     437        fd.close
     438        fd = open(test_file2, 'w')      # 5 lines, all '\n'
     439        for i in range(5):
     440            fd.write('\n')
     441        fd.close()
     442        fd = open(test_file3, 'w')      # 25 chars, no \n, 1 lines
     443        fd.write('no newline at end of line')
     444        fd.close()
     445        fd = open(test_file4, 'w')      # 1000 lines
     446        for i in range(1000):
     447            fd.write('The quick brown fox jumps over the lazy dog.\n')
     448        fd.close()
     449
     450        # use file_length() to get and check lengths
     451        size1 = file_length(test_file1)
     452        msg = 'Expected file_length() to return 0, but got %d' % size1
     453        self.failUnless(size1 == 0, msg)
     454        size2 = file_length(test_file2)
     455        msg = 'Expected file_length() to return 5, but got %d' % size2
     456        self.failUnless(size2 == 5, msg)
     457        size3 = file_length(test_file3)
     458        msg = 'Expected file_length() to return 1, but got %d' % size3
     459        self.failUnless(size3 == 1, msg)
     460        size4 = file_length(test_file4)
     461        msg = 'Expected file_length() to return 1000, but got %d' % size4
     462        self.failUnless(size4 == 1000, msg)
     463
    309464################################################################################
    310465
  • branches/numpy/anuga/utilities/util_ext.h

    r6410 r6689  
    346346
    347347
    348 // check that numpy array objects arrays are C contiguous memory
     348// check that numpy array objects are C contiguous memory
    349349#define CHECK_C_CONTIG(varname) if (!PyArray_ISCONTIGUOUS(varname)) { \
    350350                                    char msg[1024]; \
Note: See TracChangeset for help on using the changeset viewer.