Changeset 6517


Ignore:
Timestamp:
Mar 16, 2009, 11:06:22 AM (14 years ago)
Author:
rwilson
Message:

Hand-merged recent changes in main trunk. Still work to be done in shallow_water.

Location:
branches/numpy/anuga
Files:
22 edited

Legend:

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

    r6481 r6517  
    319319        return self.mesh.get_boundary_polygon(*args, **kwargs)
    320320
     321    # FIXME(Ole): This doesn't seem to be required
    321322    def get_number_of_triangles_per_node(self, *args, **kwargs):
    322323        return self.mesh.get_number_of_triangles_per_node(*args, **kwargs)
     
    348349    def statistics(self, *args, **kwargs):
    349350        return self.mesh.statistics(*args, **kwargs)
     351       
     352    def get_extent(self, *args, **kwargs):
     353        return self.mesh.get_extent(*args, **kwargs)   
    350354
    351355    ##
  • branches/numpy/anuga/abstract_2d_finite_volumes/generic_boundary_conditions.py

    r6304 r6517  
    8888    # FIXME (Ole): We should rename f to function to be consistent with
    8989    # Transmissive_Momentum_Set_Stage_Boundary (cf posting by rrraman)
    90     def __init__(self, domain = None, f = None):
     90    def __init__(self, domain=None,
     91                 f=None,
     92                 default_boundary=None,
     93                 verbose=False):
    9194        Boundary.__init__(self)
     95        self.default_boundary = default_boundary
     96        self.default_boundary_invoked = False    # Flag
     97        self.domain = domain
     98        self.verbose = verbose
    9299
    93100        try:
     
    122129    def evaluate(self, vol_id=None, edge_id=None):
    123130        # FIXME (Ole): I think this should be get_time(), see ticket:306
    124         return self.f(self.domain.time)
     131        try:
     132            res = self.f(self.domain.time)
     133        except Modeltime_too_early, e:
     134            raise Modeltime_too_early, e
     135        except Modeltime_too_late, e:
     136            if self.default_boundary is None:
     137                raise Exception, e # Reraise exception
     138            else:
     139                # Pass control to default boundary
     140                res = self.default_boundary.evaluate(vol_id, edge_id)
     141               
     142                # Ensure that result cannot be manipulated
     143                # This is a real danger in case the
     144                # default_boundary is a Dirichlet type
     145                # for instance.
     146                res = res.copy()
     147               
     148                if self.default_boundary_invoked is False:
     149                    if self.verbose:               
     150                        # Issue warning the first time
     151                        msg = '%s' %str(e)
     152                        msg += 'Instead I will use the default boundary: %s\n'\
     153                            %str(self.default_boundary)
     154                        msg += 'Note: Further warnings will be supressed'
     155                        print msg
     156               
     157                    # FIXME (Ole): Replace this crude flag with
     158                    # Python's ability to print warnings only once.
     159                    # See http://docs.python.org/lib/warning-filter.html
     160                    self.default_boundary_invoked = True
     161
     162        return res
    125163
    126164
     
    299337                    if self.default_boundary_invoked is False:
    300338                        # Issue warning the first time
    301                         msg = '%s' %str(e)
    302                         msg += 'Instead I will use the default boundary: %s\n'\
    303                             %str(self.default_boundary)
    304                         msg += 'Note: Further warnings will be supressed'
    305                         warn(msg)
     339                        if self.verbose:
     340                            msg = '%s' %str(e)
     341                            msg += 'Instead I will use the default boundary: %s\n'\
     342                                %str(self.default_boundary)
     343                            msg += 'Note: Further warnings will be supressed'
     344                            print msg
    306345                   
    307346                        # FIXME (Ole): Replace this crude flag with
  • branches/numpy/anuga/abstract_2d_finite_volumes/quantity.py

    r6441 r6517  
    886886        # a crash in fittng, so disabled it until I can investigate further
    887887        # Sorry. 23 Jan 2009. Logged as ticket:314
    888         if False:
     888        #if True: # Test will fail (31 Jan 2009)
     889        if False: # Test will pass (31 Jan 2009)
    889890            # Use mesh as defined by domain
    890891            # This used to cause problems for caching due to quantities
  • branches/numpy/anuga/abstract_2d_finite_volumes/test_domain.py

    r6410 r6517  
    542542                        conserved_quantities =\
    543543                        ['stage', 'xmomentum', 'ymomentum'])
     544        domain.set_default_order(1)
    544545        domain.check_integrity()
    545546
     
    650651                        conserved_quantities =\
    651652                        ['stage', 'xmomentum', 'ymomentum'])
     653        domain.set_default_order(1)                       
    652654        domain.check_integrity()
    653655
     
    800802                             [ 11.0,  11.0,  11.0],
    801803                             [ 11.0,  11.0,  11.0]])
     804                             
     805    def test_that_mesh_methods_exist(self):
     806        """test_that_mesh_methods_exist
     807       
     808        Test that relavent mesh methods are made available in
     809        domain through composition
     810        """
     811        from mesh_factory import rectangular
     812        from shallow_water import Domain
     813
     814        # Create basic mesh
     815        points, vertices, boundary = rectangular(1, 3)
     816
     817        # Create shallow water domain
     818        domain = Domain(points, vertices, boundary)                             
     819       
     820       
     821        domain.get_centroid_coordinates()
     822        domain.get_radii()
     823        domain.get_areas()
     824        domain.get_area()
     825        domain.get_vertex_coordinates()
     826        domain.get_triangles()
     827        domain.get_nodes()
     828        domain.get_number_of_nodes()
     829        domain.get_normal(0,0)
     830        domain.get_intersecting_segments([[0.0, 0.0], [0.0, 1.0]])
     831        domain.get_disconnected_triangles()
     832        domain.get_boundary_tags()
     833        domain.get_boundary_polygon()
     834        #domain.get_number_of_triangles_per_node()
     835        domain.get_triangles_and_vertices_per_node()
     836        domain.get_interpolation_object()
     837        domain.get_tagged_elements()
     838        domain.get_lone_vertices()
     839        domain.get_unique_vertices()
     840        g = domain.get_georeference()
     841        domain.set_georeference(g)
     842        domain.build_tagged_elements_dictionary()
     843        domain.statistics()
     844        domain.get_extent()
     845
     846       
     847       
    802848
    803849#-------------------------------------------------------------
  • branches/numpy/anuga/abstract_2d_finite_volumes/test_general_mesh.py

    r6463 r6517  
    9191                  str(verts)))
    9292        self.assert_(num.allclose(num.array([nodes_absolute[1],
    93                                              nodes_absolute[0],
    94                                              nodes_absolute[2]]),
    95                                   verts), msg)
     93                                     nodes_absolute[0],
     94                                     nodes_absolute[2]]), verts))
     95        verts = domain.get_vertex_coordinates(triangle_id=0,
     96                                              absolute=True)       
     97        self.assert_(num.allclose(num.array([nodes_absolute[1],
     98                                     nodes_absolute[0],
     99                                     nodes_absolute[2]]), verts))
     100       
     101       
    96102
    97103    def test_get_vertex_coordinates_triangle_id(self):
  • branches/numpy/anuga/abstract_2d_finite_volumes/util.py

    r6441 r6517  
    925925                    verbose = False):   
    926926       
     927    # FIXME(Ole): Shouldn't print statements here be governed by verbose?
    927928    assert type(gauge_filename) == type(''), 'Gauge filename must be a string'
    928929   
     
    971972            raise msg
    972973
    973         print 'swwfile', swwfile
     974        if verbose:
     975            print 'swwfile', swwfile
    974976
    975977        # Extract parent dir name and use as label
     
    24892491                                 base_name=base,
    24902492                                 verbose=verbose)
     2493    #print 'sww files just after get_all_swwfiles()', sww_files
     2494    # fudge to get SWW files in 'correct' order, oldest on the left
     2495    sww_files.sort()
     2496
     2497    if verbose:
     2498        print 'sww files', sww_files
    24912499   
    24922500    #to make all the quantities lower case for file_function
     
    24972505
    24982506    core_quantities = ['stage', 'elevation', 'xmomentum', 'ymomentum']
    2499    
     2507    gauge_file = out_name
     2508
     2509    heading = [quantity for quantity in quantities]
     2510    heading.insert(0,'time')
     2511    heading.insert(1,'hours')
     2512
     2513    #create a list of csv writers for all the points and write header
     2514    points_writer = []
     2515    for point_i,point in enumerate(points):
     2516        points_writer.append(writer(file(dir_name + sep + gauge_file
     2517                                         + point_name[point_i] + '.csv', "wb")))
     2518        points_writer[point_i].writerow(heading)
     2519   
     2520    if verbose: print 'Writing csv files'
     2521
     2522    quake_offset_time = None
     2523
    25002524    for sww_file in sww_files:
    25012525        sww_file = join(dir_name, sww_file+'.sww')
  • branches/numpy/anuga/alpha_shape/alpha_shape.py

    r6304 r6517  
    289289                       (denom[k]< EPSILON and  denom[k] > -EPSILON)]
    290290
    291         if num.any(denom == 0.0):
    292             raise AlphaError
    293 
    294         dx = num.divide(y31*dist21 - y21*dist31, denom)
    295         dy = num.divide(x21*dist31 - x31*dist21, denom)
    296 
     291        if num.alltrue(denom != 0.0):               
     292            dx = num.divide(y31*dist21 - y21*dist31,denom)
     293            dy = num.divide(x21*dist31 - x31*dist21,denom)
     294        else:
     295            raise  AlphaError
     296           
    297297        self.triradius = 0.5*num.sqrt(dx*dx + dy*dy)
    298298        #print "triangle radii", self.triradius
  • branches/numpy/anuga/caching/test_caching.py

    r6410 r6517  
    386386
    387387        # Check that hash value of callable objects don't change
     388        # FIXME (Ole): The hash values do appear to change when OS
     389        # and/or dependencies are upgraded
    388390        if os.name == 'posix' and os.uname()[4] in ['x86_64', 'ia64']:
    389391          # 64 bit hash values
    390           f1hash = 1914027059797211698
    391           f2hash = 1914027059807087171
     392          f1hash = 7079146893884768701
     393          f2hash = -6995306676314913340
     394
     395          # Prior to cluster upgrades Feb 2009
     396          #f1hash = 1914027059797211698
     397          #f2hash = 1914027059807087171
    392398        else:
     399          # 32 bit hash values
    393400          f1hash = -758136387
    394401          f2hash = -11221564     
  • branches/numpy/anuga/compile_all.py

    r4978 r6517  
    2727#execfile('test_all.py')
    2828   
     29if sys.platform == 'win32':
     30    raw_input('Press any key')
  • branches/numpy/anuga/config.py

    r6361 r6517  
    9191       
    9292# FIXME (Ole) Maybe get rid of order altogether and use beta_w
    93 # ... and isn't it about time we make the default 2?
    94 default_order = 1
     93default_order = 2
    9594
    9695################################################################################
  • branches/numpy/anuga/coordinate_transforms/redfearn.py

    r6149 r6517  
    3737    return sign*dd, mm, ss
    3838
    39 def redfearn(lat, lon, false_easting=None, false_northing=None, zone=None):
     39def redfearn(lat, lon, false_easting=None, false_northing=None,
     40             zone=None, central_meridian=None, scale_factor=None):
    4041    """Compute UTM projection using Redfearn's formula
    4142
     
    4748    If zone is specified reproject lat and long to specified zone instead of
    4849    standard zone
     50    If meridian is specified, reproject lat and lon to that instead of zone. In this case
     51    zone will be set to -1 to indicate non-UTM projection
     52
     53    Note that zone and meridian cannot both be specifed
    4954    """
    5055
     
    5762    a = 6378137.0                       #Semi major axis
    5863    inverse_flattening = 298.257222101  #1/f
    59     K0 = 0.9996                         #Central scale factor   
     64    if scale_factor is None:
     65        K0 = 0.9996                         #Central scale factor   
     66    else:
     67        K0 = scale_factor
     68    #print 'scale', K0
    6069    zone_width = 6                      #Degrees
    6170
     
    138147    m = term1 + term2 + term3 + term4 #OK
    139148
    140     #Zone
     149    if zone is not None and central_meridian is not None:
     150        msg = 'You specified both zone and central_meridian. Provide only one of them'
     151        raise Exception, msg
     152   
     153    # Zone
    141154    if zone is None:
    142155        zone = int((lon - longitude_of_western_edge_zone0)/zone_width)
    143156
    144     central_meridian = zone*zone_width+longitude_of_central_meridian_zone0
     157    # Central meridian
     158    if central_meridian is None:
     159        central_meridian = zone*zone_width+longitude_of_central_meridian_zone0
     160    else:
     161        zone = -1
    145162
    146163    omega = (lon-central_meridian)*pi/180 #Relative longitude (radians)
  • branches/numpy/anuga/coordinate_transforms/test_geo_reference.py

    r6442 r6517  
    235235        y = 3.0
    236236       
    237         g = Geo_reference(56, x, y)
    238         points = [[3.0,34.0], [64.0,6.0]]
    239         new_points = g.get_absolute(points)
    240 
    241         self.failUnless(isinstance(new_points, list), 'failed')
    242         self.failUnless(type(new_points) == type(points), 'failed')
    243         for point, new_point in map(None, points, new_points):
    244             self.failUnless(point[0]+x == new_point[0], 'failed')
    245             self.failUnless(point[1]+y == new_point[1], 'failed')
    246 
    247         # test with no supplied offsets
     237        g = Geo_reference(56,x,y)
     238        lofl = [[3.0,34.0], [64.0,6.0]]
     239        new_lofl = g.get_absolute(lofl)
     240        #print "lofl",lofl
     241        #print "new_lofl",new_lofl
     242
     243        self.failUnless(type(new_lofl) == types.ListType, ' failed')
     244        self.failUnless(type(new_lofl) == type(lofl), ' failed')
     245        for point,new_point in map(None,lofl,new_lofl):
     246            self.failUnless(point[0]+x==new_point[0], ' failed')
     247            self.failUnless(point[1]+y==new_point[1], ' failed')
     248
     249           
    248250        g = Geo_reference()
    249         points = [[3.0,34.0], [64.0,6.0]]
    250         new_points = g.get_absolute(points)
    251 
    252         self.failUnless(isinstance(new_points, list), 'failed')
    253         self.failUnless(type(new_points) == type(points), 'failed')
    254         for point, new_point in map(None, points, new_points):
    255             self.failUnless(point[0] == new_point[0], 'failed')
    256             self.failUnless(point[1] == new_point[1], 'failed')
     251        lofl = [[3.0,34.0], [64.0,6.0]]
     252        new_lofl = g.get_absolute(lofl)
     253        #print "lofl",lofl
     254        #print "new_lofl",new_lofl
     255
     256        self.failUnless(type(new_lofl) == types.ListType, ' failed')
     257        self.failUnless(type(new_lofl) == type(lofl), ' failed')
     258        for point,new_point in map(None,lofl,new_lofl):
     259            self.failUnless(point[0]==new_point[0], ' failed')
     260            self.failUnless(point[1]==new_point[1], ' failed')
    257261           
    258262        # test that calling get_absolute twice does the right thing
  • branches/numpy/anuga/coordinate_transforms/test_redfearn.py

    r6410 r6517  
    122122
    123123    def test_UTM_6_nonstandard_projection(self):
    124         #Test 6 (Geraldton, WA)
    125 
    126         #First test native projection (zone 50)
     124        """test_UTM_6_nonstandard_projection
     125
     126        Test that projections can be forced to
     127        use other than native zone.
     128
     129        Data is from Geraldton, WA
     130        """
     131       
     132
     133        # First test native projection (zone 50)
    127134        zone, easting, northing = redfearn(-29.233299999,114.05)
    128135
     
    267274    #    #assert allclose(northing, 6181725.1724276)
    268275
     276   
     277    def Xtest_nonstandard_meridian_coinciding_with_native(self):
     278        """test_nonstandard_meridian_coinciding_with_native
     279
     280        This test will verify that redfearn can be used to project
     281        points using an arbitrary central meridian that happens to
     282        coincide with the standard meridian at the center of a UTM zone.
     283        This is a preliminary test before testing this functionality
     284        with a truly arbitrary non-standard meridian.
     285        """
     286
     287        # The file projection_test_points_z53.csv contains 10 points
     288        # which straddle the boundary between UTM zones 53 and 54.
     289        # They have been projected to zone 53 irrespective of where they
     290        # belong.
     291
     292        path = get_pathname_from_package('anuga.coordinate_transforms')
     293       
     294        for forced_zone in [53, 54]:
     295       
     296            datafile = join(path, 'projection_test_points_z%d.csv' % forced_zone)
     297            fid = open(datafile)
     298
     299            for line in fid.readlines()[1:]:
     300                fields = line.strip().split(',')
     301               
     302                lon = float(fields[1])
     303                lat = float(fields[2])
     304                x = float(fields[3])
     305                y = float(fields[4])           
     306
     307                zone, easting, northing = redfearn(lat, lon,
     308                                                   zone=forced_zone)
     309               
     310                print
     311                print 'Lat', lat
     312                print 'Lon', lon
     313                print 'Zone', zone
     314                print 'Ref x', x, 'Computed x', easting, 'Close enough:', num.allclose(x, easting)
     315                print 'Ref y', y, 'Computed y', northing, 'Close enough:', num.allclose(y, northing)
     316               
     317                # Check calculation
     318                assert zone == forced_zone
     319                print
     320                #assert num.allclose(x, easting)
     321                #assert num.allclose(y, northing)
     322
     323   
     324   
     325   
     326    def Xtest_nonstandard_meridian(self):
     327        """test_nonstandard_meridian
     328
     329        This test will verify that redfearn can be used to project
     330        points using an arbitrary central meridian.
     331        """
     332
     333        # The file projection_test_points.csv contains 10 points
     334        # which straddle the boundary between UTM zones 53 and 54.
     335        # They have been projected using a central meridian of 137.5
     336        # degrees (the boundary is 138 so it is pretty much right
     337        # in the middle of zones 53 and 54).
     338
     339        path = get_pathname_from_package('anuga.coordinate_transforms')
     340        datafile = join(path, 'projection_test_points.csv')
     341        fid = open(datafile)
     342
     343        for line in fid.readlines()[1:]:
     344            fields = line.strip().split(',')
     345
     346            lon = float(fields[1])
     347            lat = float(fields[2])
     348            x = float(fields[3])
     349            y = float(fields[4])           
     350
     351            zone, easting, northing = redfearn(lat, lon,
     352                                               central_meridian=137.5,
     353                                               scale_factor=0.9996)
     354
     355            print
     356            print 'Lat', lat
     357            print 'Lon', lon
     358            print 'Zone', zone
     359            print 'Ref x', x, 'Computed x', easting, 'Close enough:', num.allclose(x, easting)
     360            print 'Ref y', y, 'Computed y', northing, 'Close enough:', num.allclose(y, northing)
     361
     362            # Check calculation
     363            assert zone == -1 # Indicates non UTM projection
     364            print
     365            #assert num.allclose(x, easting)
     366            #assert num.allclose(y, northing)
     367
     368        # Test that zone and meridian can't both be specified
     369        try:
     370            zone, easting, northing = redfearn(lat, lon,
     371                                               zone=50,
     372                                               central_meridian=137.5)
     373        except:
     374            pass
     375        else:
     376            msg = 'Should have raised exception'
     377            raise Exception, msg
     378
     379           
    269380    def test_convert_lats_longs(self):
    270381
  • branches/numpy/anuga/culvert_flows/Test_Culvert_Flat_Water_Lev.py

    r6304 r6517  
    123123                       number_of_barrels=1,
    124124                       update_interval=2,
     125                       log_file=True,
     126                       discharge_hydrograph=True,
    125127                       verbose=True)
    126128
  • branches/numpy/anuga/culvert_flows/culvert_class.py

    r6304 r6517  
    204204            log_to_file(self.log_filename, description)
    205205            log_to_file(self.log_filename, self.culvert_type)       
     206        else:
     207            self.log_filename = None
    206208
    207209
     
    296298
    297299        # Print some diagnostics to log if requested
    298         if hasattr(self, 'log_filename'):
     300        if self.log_filename is not None:
    299301            s = 'Culvert Effective Length = %.2f m' %(self.length)
    300302            log_to_file(self.log_filename, s)
     
    324326                update = True
    325327
    326         if hasattr(self, 'log_filename'):           
     328           
     329        if self.log_filename is not None:       
    327330            s = '\nTime = %.2f, delta_t = %f' %(time, delta_t)
    328331            log_to_file(self.log_filename, s)
     
    439442            if self.verbose is True:
    440443                print msg
    441             if hasattr(self, 'log_filename'):                   
     444               
     445            if self.log_filename is not None:               
    442446                log_to_file(self.log_filename, msg)
    443447       
     
    524528                print '%.2fs - WARNING: Flow is running uphill.' % time
    525529           
    526         if hasattr(self, 'log_filename'):
     530        if self.log_filename is not None:
    527531            s = 'Time=%.2f, inlet stage = %.2f, outlet stage = %.2f'\
    528532                %(time, self.inlet.stage, self.outlet.stage)
     
    573577                    msg += 'for culvert "%s"' % self.label
    574578                   
    575                     if hasattr(self, 'log_filename'):                   
     579                    if self.log_filename is not None:                   
    576580                        log_to_file(self.log_filename, msg)
    577581            else:
    578582                # User culvert routine
    579583                Q, barrel_velocity, culvert_outlet_depth =\
    580                     self.culvert_routine(self, delta_total_energy, g)
     584                    self.culvert_routine(inlet.depth,
     585                                         outlet.depth,
     586                                         inlet.specific_energy,
     587                                         delta_total_energy,
     588                                         g,
     589                                         culvert_length=self.length,
     590                                         culvert_width=self.width,
     591                                         culvert_height=self.height,
     592                                         culvert_type=self.culvert_type,
     593                                         manning=self.manning,
     594                                         sum_loss=self.sum_loss,
     595                                         log_filename=self.log_filename)
    581596           
    582597           
     
    601616
    602617                           
    603    
     618# OBSOLETE (Except for momentum jet in Culvert_flow_energy)   
    604619class Culvert_flow_rating:
    605620    """Culvert flow - transfer water from one hole to another
  • branches/numpy/anuga/culvert_flows/culvert_routines.py

    r6128 r6517  
    88
    99#NOTE:
    10 # Inlet control:  Delta_Et > Es at the inlet
    11 # Outlet control: Delta_Et < Es at the inlet
    12 # where Et is total energy (w + 0.5*v^2/g) and
    13 # Es is the specific energy (h + 0.5*v^2/g)
    14 
    15 
    16 
    17 #   NEW DEFINED CULVERT FLOW---- Flow from INLET 1 ------> INLET 2 (Outlet)
    18 #
    19 # The First Attempt has a Simple Horizontal Circle as a Hole on the Bed
    20 # Flow Is Removed at a Rate of INFLOW
    21 # Downstream there is a similar Circular Hole on the Bed where INFLOW effectively Surcharges
    22 #
    23 # This SHould be changed to a Vertical Opening Both BOX and Circular
    24 # There will be several Culvert Routines such as:
    25 # CULVERT_Boyd_Channel
    26 # CULVERT_Orifice_and_Weir
    27 # CULVERT_Simple_FLOOR
    28 # CULVERT_Simple_WALL
    29 # CULVERT_Eqn_Floor
    30 # CULVERT_Eqn_Wall
    31 # CULVERT_Tab_Floor
    32 # CULVERT_Tab_Wall
    33 # BRIDGES.....
    34 # NOTE NEED TO DEVELOP CONCEPT 1D Model for Linked Pipe System !!!!
    35 
    36 #  COULD USE EPA SWMM Model !!!!
     10# Inlet control:  Delta_total_energy > inlet_specific_energy
     11# Outlet control: Delta_total_energy < inlet_specific_energy
     12# where total energy is (w + 0.5*v^2/g) and
     13# specific energy is (h + 0.5*v^2/g)
    3714
    3815
     
    4017
    4118
    42 def boyd_generalised_culvert_model(culvert, delta_total_energy, g):
    43 
    44     """Boyd's generalisation of the US department of transportation culvert model
    45         # == The quantity of flow passing through a culvert is controlled by many factors
    46         # == It could be that the culvert is controled by the inlet only, with it being Un submerged this is effectively equivalent to the WEIR Equation
    47         # == Else the culvert could be controlled by the inlet, with it being Submerged, this is effectively the Orifice Equation
    48         # == Else it may be controlled by Down stream conditions where depending on the down stream depth, the momentum in the culvert etc. flow is controlled
     19def boyd_generalised_culvert_model(inlet_depth,
     20                                   outlet_depth,
     21                                   inlet_specific_energy,
     22                                   delta_total_energy,
     23                                   g,
     24                                   culvert_length=0.0,
     25                                   culvert_width=0.0,
     26                                   culvert_height=0.0,
     27                                   culvert_type='box',
     28                                   manning=0.0,
     29                                   sum_loss=0.0,
     30                                   log_filename=None):
     31
     32    """Boyd's generalisation of the US department of transportation culvert
     33    model
     34     
     35    The quantity of flow passing through a culvert is controlled by many factors
     36    It could be that the culvert is controlled by the inlet only, with it
     37    being unsubmerged this is effectively equivalent to the weir Equation
     38    Else the culvert could be controlled by the inlet, with it being
     39    submerged, this is effectively the Orifice Equation
     40    Else it may be controlled by down stream conditions where depending on
     41    the down stream depth, the momentum in the culvert etc. flow is controlled
    4942    """
    5043
     
    5346    from anuga.utilities.numerical_tools import safe_acos as acos
    5447
    55     inlet = culvert.inlet
    56     outlet = culvert.outlet       
    57     Q_outlet_tailwater = 0.0
    58     inlet.rate = 0.0
    59     outlet.rate = 0.0
    60     Q_inlet_unsubmerged = 0.0
    61     Q_inlet_submerged = 0.0
    62     Q_outlet_critical_depth = 0.0
    63 
    64     if hasattr(culvert, 'log_filename'):                               
    65         log_filename = culvert.log_filename
    66 
    67     manning = culvert.manning
    68     sum_loss = culvert.sum_loss
    69     length = culvert.length
    70 
    71     if inlet.depth > 0.01:
    72 
    73         E = inlet.specific_energy
    74 
    75         if hasattr(culvert, 'log_filename'):                           
    76             s = 'Specific energy  = %f m' %E
     48    if inlet_depth > 0.01:
     49        # Water has risen above inlet
     50       
     51        if log_filename is not None:                           
     52            s = 'Specific energy  = %f m' % inlet_specific_energy
    7753            log_to_file(log_filename, s)
    7854       
    79         msg = 'Specific energy is negative'
    80         assert E >= 0.0, msg
     55        msg = 'Specific energy at inlet is negative'
     56        assert inlet_specific_energy >= 0.0, msg
    8157                     
    82        
    83         # Water has risen above inlet
    84         if culvert.culvert_type == 'circle':
    85             # Round culvert
     58        if culvert_type == 'circle':
     59            # Round culvert (use width as diameter)
     60            diameter = culvert_width           
    8661
    8762            # Calculate flows for inlet control           
    88             diameter = culvert.diameter
    89 
    90             Q_inlet_unsubmerged = 0.421*g**0.5*diameter**0.87*E**1.63    # Inlet Ctrl Inlet Unsubmerged
    91             Q_inlet_submerged = 0.530*g**0.5*diameter**1.87*E**0.63      # Inlet Ctrl Inlet Submerged
    92 
    93             if hasattr(culvert, 'log_filename'):
    94                 s = 'Q_inlet_unsubmerged = %.6f, Q_inlet_submerged = %.6f' %(Q_inlet_unsubmerged, Q_inlet_submerged)
     63            Q_inlet_unsubmerged = 0.421*g**0.5*diameter**0.87*inlet_specific_energy**1.63 # Inlet Ctrl Inlet Unsubmerged
     64            Q_inlet_submerged = 0.530*g**0.5*diameter**1.87*inlet_specific_energy**0.63   # Inlet Ctrl Inlet Submerged
     65
     66            if log_filename is not None:                                       
     67                s = 'Q_inlet_unsubmerged = %.6f, Q_inlet_submerged = %.6f' % (Q_inlet_unsubmerged, Q_inlet_submerged)
    9568                log_to_file(log_filename, s)
    9669
    97             case = ''
     70            # FIXME(Ole): Are these functions really for inlet control?                   
    9871            if Q_inlet_unsubmerged < Q_inlet_submerged:
    9972                Q = Q_inlet_unsubmerged
    10073
    101                 alpha = acos(1 - inlet.depth/diameter)
     74                alpha = acos(1 - inlet_depth/diameter)
    10275                flow_area = diameter**2 * (alpha - sin(alpha)*cos(alpha))
    103                 outlet_culvert_depth = inlet.depth
    104                 width = diameter*sin(alpha)
    105                 #perimeter = alpha*diameter               
     76                outlet_culvert_depth = inlet_depth
    10677                case = 'Inlet unsubmerged'
    10778            else:   
     
    10980                flow_area = (diameter/2)**2 * pi
    11081                outlet_culvert_depth = diameter
    111                 width = diameter
    112                 #perimeter = diameter
    11382                case = 'Inlet submerged'                   
    11483
    11584
    11685
    117             if delta_total_energy < E:
     86            if delta_total_energy < inlet_specific_energy:
    11887                # Calculate flows for outlet control
    11988                # Determine the depth at the outlet relative to the depth of flow in the Culvert
    12089
    121                 if outlet.depth > diameter:        # The Outlet is Submerged
     90                if outlet_depth > diameter:        # The Outlet is Submerged
    12291                    outlet_culvert_depth=diameter
    12392                    flow_area = (diameter/2)**2 * pi  # Cross sectional area of flow in the culvert
    12493                    perimeter = diameter * pi
    125                     width = diameter
    12694                    case = 'Outlet submerged'
    127                 elif outlet.depth==0.0:
    128                     outlet_culvert_depth=inlet.depth  # For purpose of calculation assume the outlet depth = the inlet depth
    129                     alpha = acos(1 - inlet.depth/diameter)
     95                elif outlet_depth==0.0:
     96                    outlet_culvert_depth=inlet_depth # For purpose of calculation assume the outlet depth = the inlet depth
     97                    alpha = acos(1 - inlet_depth/diameter)
    13098                    flow_area = diameter**2 * (alpha - sin(alpha)*cos(alpha))
    13199                    perimeter = alpha*diameter
    132                     width = diameter*sin(alpha)                   
    133100
    134101                    case = 'Outlet depth is zero'                       
    135102                else:   # Here really should use the Culvert Slope to calculate Actual Culvert Depth & Velocity
    136                     outlet_culvert_depth=outlet.depth   # For purpose of calculation assume the outlet depth = the inlet depth
    137                     alpha = acos(1 - outlet.depth/diameter)
     103                    outlet_culvert_depth=outlet_depth   # For purpose of calculation assume the outlet depth = the inlet depth
     104                    alpha = acos(1 - outlet_depth/diameter)
    138105                    flow_area = diameter**2 * (alpha - sin(alpha)*cos(alpha))
    139106                    perimeter = alpha*diameter
    140                     width = diameter*sin(alpha)                   
    141107                    case = 'Outlet is open channel flow'
    142108
    143109                hyd_rad = flow_area/perimeter
    144110               
    145                 if hasattr(culvert, 'log_filename'):
     111                if log_filename is not None:
    146112                    s = 'hydraulic radius at outlet = %f' %hyd_rad
    147113                    log_to_file(log_filename, s)
    148114
    149115                # Outlet control velocity using tail water
    150                 culvert_velocity = sqrt(delta_total_energy/((sum_loss/2*g)+(manning**2*length)/hyd_rad**1.33333))
     116                culvert_velocity = sqrt(delta_total_energy/((sum_loss/2*g)+(manning**2*culvert_length)/hyd_rad**1.33333))
    151117                Q_outlet_tailwater = flow_area * culvert_velocity
    152118               
    153                 if hasattr(culvert, 'log_filename'):
     119                if log_filename is not None:               
    154120                    s = 'Q_outlet_tailwater = %.6f' %Q_outlet_tailwater
    155121                    log_to_file(log_filename, s)
     
    165131
    166132            # Calculate flows for inlet control
    167             height = culvert.height
    168             width = culvert.width           
     133            height = culvert_height
     134            width = culvert_width           
    169135           
    170             Q_inlet_unsubmerged = 0.540*g**0.5*width*E**1.50 # Flow based on Inlet Ctrl Inlet Unsubmerged
    171             Q_inlet_submerged = 0.702*g**0.5*width*height**0.89*E**0.61  # Flow based on Inlet Ctrl Inlet Submerged
    172 
    173             if hasattr(culvert, 'log_filename'):
     136            Q_inlet_unsubmerged = 0.540*g**0.5*width*inlet_specific_energy**1.50 # Flow based on Inlet Ctrl Inlet Unsubmerged
     137            Q_inlet_submerged = 0.702*g**0.5*width*height**0.89*inlet_specific_energy**0.61  # Flow based on Inlet Ctrl Inlet Submerged
     138
     139            if log_filename is not None:                           
    174140                s = 'Q_inlet_unsubmerged = %.6f, Q_inlet_submerged = %.6f' %(Q_inlet_unsubmerged, Q_inlet_submerged)
    175141                log_to_file(log_filename, s)
    176142
    177             case = ''
     143            # FIXME(Ole): Are these functions really for inlet control?   
    178144            if Q_inlet_unsubmerged < Q_inlet_submerged:
    179145                Q = Q_inlet_unsubmerged
    180                 flow_area = width*inlet.depth
    181                 outlet_culvert_depth = inlet.depth
    182                 #perimeter=(width+2.0*inlet.depth)               
     146                flow_area = width*inlet_depth
     147                outlet_culvert_depth = inlet_depth
    183148                case = 'Inlet unsubmerged'
    184149            else:   
     
    189154                case = 'Inlet submerged'                   
    190155
    191             if delta_total_energy < E:
     156            if delta_total_energy < inlet_specific_energy:
    192157                # Calculate flows for outlet control
    193158                # Determine the depth at the outlet relative to the depth of flow in the Culvert
    194159
    195                 if outlet.depth > height:        # The Outlet is Submerged
     160                if outlet_depth > height:        # The Outlet is Submerged
    196161                    outlet_culvert_depth=height
    197162                    flow_area=width*height       # Cross sectional area of flow in the culvert
    198163                    perimeter=2.0*(width+height)
    199164                    case = 'Outlet submerged'
    200                 elif outlet.depth==0.0:
    201                     outlet_culvert_depth=inlet.depth   # For purpose of calculation assume the outlet depth = the inlet depth
    202                     flow_area=width*inlet.depth
    203                     perimeter=(width+2.0*inlet.depth)
     165                elif outlet_depth==0.0:
     166                    outlet_culvert_depth=inlet_depth   # For purpose of calculation assume the outlet depth = the inlet depth
     167                    flow_area=width*inlet_depth
     168                    perimeter=(width+2.0*inlet_depth)
    204169                    case = 'Outlet depth is zero'                       
    205170                else:   # Here really should use the Culvert Slope to calculate Actual Culvert Depth & Velocity
    206                     outlet_culvert_depth=outlet.depth
    207                     flow_area=width*outlet.depth
    208                     perimeter=(width+2.0*outlet.depth)
     171                    outlet_culvert_depth=outlet_depth
     172                    flow_area=width*outlet_depth
     173                    perimeter=(width+2.0*outlet_depth)
    209174                    case = 'Outlet is open channel flow'
    210175
    211176                hyd_rad = flow_area/perimeter
    212177               
    213                 if hasattr(culvert, 'log_filename'):               
    214                     s = 'hydraulic radius at outlet = %f' %hyd_rad
     178                if log_filename is not None:                               
     179                    s = 'hydraulic radius at outlet = %f' % hyd_rad
    215180                    log_to_file(log_filename, s)
    216181
    217182                # Outlet control velocity using tail water
    218                 culvert_velocity = sqrt(delta_total_energy/((sum_loss/2*g)+(manning**2*length)/hyd_rad**1.33333))
     183                culvert_velocity = sqrt(delta_total_energy/((sum_loss/2*g)+(manning**2*culvert_length)/hyd_rad**1.33333))
    219184                Q_outlet_tailwater = flow_area * culvert_velocity
    220185
    221                 if hasattr(culvert, 'log_filename'):                           
    222                     s = 'Q_outlet_tailwater = %.6f' %Q_outlet_tailwater
     186                if log_filename is not None:                           
     187                    s = 'Q_outlet_tailwater = %.6f' % Q_outlet_tailwater
    223188                    log_to_file(log_filename, s)
    224189                Q = min(Q, Q_outlet_tailwater)
     
    228193
    229194        # Common code for circle and square geometries
    230         if hasattr(culvert, 'log_filename'):                                   
    231             log_to_file(log_filename, 'Case: "%s"' %case)
     195        if log_filename is not None:
     196            log_to_file(log_filename, 'Case: "%s"' % case)
    232197           
    233         flow_rate_control=Q
    234 
    235         if hasattr(culvert, 'log_filename'):                                   
    236             s = 'Flow Rate Control = %f' %flow_rate_control
     198        if log_filename is not None:                       
     199            s = 'Flow Rate Control = %f' % Q
    237200            log_to_file(log_filename, s)
    238201
    239         inlet.rate = -flow_rate_control
    240         outlet.rate = flow_rate_control               
    241            
    242         culv_froude=sqrt(flow_rate_control**2*width/(g*flow_area**3))
    243        
    244         if hasattr(culvert, 'log_filename'):                           
    245             s = 'Froude in Culvert = %f' %culv_froude
     202        culv_froude=sqrt(Q**2*width/(g*flow_area**3))
     203       
     204        if log_filename is not None:                           
     205            s = 'Froude in Culvert = %f' % culv_froude
    246206            log_to_file(log_filename, s)
    247207
     
    250210
    251211
    252     else: #inlet.depth < 0.01:
     212    else: # inlet_depth < 0.01:
    253213        Q = barrel_velocity = outlet_culvert_depth = 0.0
    254214       
  • branches/numpy/anuga/fit_interpolate/fit.py

    r6304 r6517  
    310310                assert flag is True, msg               
    311311               
     312                # FIXME(Ole): This is the message referred to in ticket:314
    312313                minx = min(self.mesh_boundary_polygon[:,0])
    313314                maxx = max(self.mesh_boundary_polygon[:,0])               
     
    317318                msg += 'Mesh boundary extent is (%.f, %.f), (%.f, %.f)'\
    318319                    % (minx, maxx, miny, maxy)
     320                raise RuntimeError, msg
    319321               
    320322
     
    375377                self.build_fit_subset(points, z, verbose=verbose)
    376378
     379                # FIXME(Ole): I thought this test would make sense here
     380                # See test_fitting_example_that_crashed_2 in test_shallow_water_domain.py
     381                # Committed 11 March 2009
     382                msg = 'Matrix AtA was not built'
     383                assert self.AtA is not None, msg
     384               
     385                #print 'Matrix was built OK'
     386
    377387               
    378388            point_coordinates = None
     
    382392        if point_coordinates is None:
    383393            if verbose: print 'Warning: no data points in fit'
    384             #assert self.AtA != None, 'no interpolation matrix'
    385             #assert self.Atz != None
    386             assert not self.AtA is None, 'no interpolation matrix'
    387             assert not self.Atz is None
     394            msg = 'No interpolation matrix.'
     395            assert self.AtA is not None, msg
     396            assert self.Atz is not None
    388397           
    389398            # FIXME (DSG) - do  a message
  • branches/numpy/anuga/geospatial_data/geospatial_data.py

    r6441 r6517  
    18061806    points = G_small.get_data_points()
    18071807
     1808    #FIXME: Remove points outside boundary polygon
     1809#    print 'new point',len(points)
     1810#
     1811#    new_points=[]
     1812#    new_points=array([],typecode=Float)
     1813#    new_points=resize(new_points,(len(points),2))
     1814#    print "BOUNDARY", boundary_poly
     1815#    for i,point in enumerate(points):
     1816#        if is_inside_polygon(point,boundary_poly, verbose=True):
     1817#            new_points[i] = point
     1818#            print"WOW",i,new_points[i]
     1819#    points = new_points
     1820
    18081821    if verbose: print "Number of points in sample to compare: ", len(points)
    18091822
  • branches/numpy/anuga/shallow_water/__init__.py

    r6190 r6517  
    1212     Transmissive_momentum_set_stage_boundary,\
    1313     Dirichlet_discharge_boundary,\
    14      Field_boundary
     14     Field_boundary,\
     15     Transmissive_stage_zero_momentum_boundary
    1516
    1617
  • branches/numpy/anuga/shallow_water/data_manager.py

    r6481 r6517  
    58635863        # This is being used to seperate one number from a list.
    58645864        # what it is actually doing is sorting lists from numeric arrays.
    5865         if type(times) is list or isinstance(times, num.ndarray):
     5865        #if type(times) is list or isinstance(times, num.ndarray):
     5866        if isinstance(times, (list, num.ndarray)):
    58665867            number_of_times = len(times)
    58675868            times = ensure_numeric(times)
     
    59305931            #outfile.variables[q+Write_sww.RANGE][1] = -max_float # Max
    59315932
    5932         if type(times) is list or isinstance(times, num.ndarray):
    5933             outfile.variables['time'][:] = times    #Store time relative
     5933        #if type(times) is list or isinstance(times, num.ndarray):
     5934        if isinstance(times, (list, num.ndarray)):
     5935            outfile.variables['time'][:] = times    # Store time relative
    59345936
    59355937        if verbose:
     
    62656267        # This is being used to seperate one number from a list.
    62666268        # what it is actually doing is sorting lists from numeric arrays.
    6267         if type(times) is list or isinstance(times, num.ndarray):
     6269        #if type(times) is list or isinstance(times, num.ndarray):
     6270        if isinstance(times, (list, num.ndarray)):
    62686271            number_of_times = len(times)
    62696272            times = ensure_numeric(times)
     
    63126315            outfile.variables[q + Write_sts.RANGE][1] = -max_float # Max
    63136316
    6314         if type(times) is list or isinstance(times, num.ndarray):
     6317        #if type(times) is list or isinstance(times, num.ndarray):
     6318        if isinstance(times, (list, num.ndarray)):
    63156319            outfile.variables['time'][:] = times    #Store time relative
    63166320
  • branches/numpy/anuga/shallow_water/test_shallow_water_domain.py

    r6451 r6517  
    16721672        # Setup computational domain
    16731673        #-----------------------------------------------------------------
    1674         points, vertices, boundary = rectangular_cross(10, 10)    # Basic mesh
    1675         domain = Domain(points, vertices, boundary)    # Create domain
     1674        points, vertices, boundary = rectangular_cross(10, 10) # Basic mesh
     1675        domain = Domain(points, vertices, boundary) # Create domain
     1676        domain.set_default_order(1)       
    16761677        domain.set_quantities_to_be_stored(None)
    16771678        domain.set_maximum_allowed_speed(100) #FIXME (Ole): try to remove this
  • branches/numpy/anuga/test_all.py

    r6304 r6517  
    215215        #os.remove(filename)
    216216
     217   
     218    if sys.platform == 'win32':
     219        raw_input('Press any key')
Note: See TracChangeset for help on using the changeset viewer.