Changeset 6533


Ignore:
Timestamp:
Mar 17, 2009, 4:02:54 PM (16 years ago)
Author:
rwilson
Message:

Revert back to 6481, prior to auto-merge of trunk and numpy branch.

Location:
branches/numpy/anuga
Files:
22 edited

Legend:

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

    r6517 r6533  
    319319        return self.mesh.get_boundary_polygon(*args, **kwargs)
    320320
    321     # FIXME(Ole): This doesn't seem to be required
    322321    def get_number_of_triangles_per_node(self, *args, **kwargs):
    323322        return self.mesh.get_number_of_triangles_per_node(*args, **kwargs)
     
    349348    def statistics(self, *args, **kwargs):
    350349        return self.mesh.statistics(*args, **kwargs)
    351        
    352     def get_extent(self, *args, **kwargs):
    353         return self.mesh.get_extent(*args, **kwargs)   
    354350
    355351    ##
  • branches/numpy/anuga/abstract_2d_finite_volumes/generic_boundary_conditions.py

    r6517 r6533  
    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,
    91                  f=None,
    92                  default_boundary=None,
    93                  verbose=False):
     90    def __init__(self, domain = None, f = None):
    9491        Boundary.__init__(self)
    95         self.default_boundary = default_boundary
    96         self.default_boundary_invoked = False    # Flag
    97         self.domain = domain
    98         self.verbose = verbose
    9992
    10093        try:
     
    129122    def evaluate(self, vol_id=None, edge_id=None):
    130123        # FIXME (Ole): I think this should be get_time(), see ticket:306
    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
     124        return self.f(self.domain.time)
    163125
    164126
     
    337299                    if self.default_boundary_invoked is False:
    338300                        # Issue warning the first time
    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
     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)
    345306                   
    346307                        # FIXME (Ole): Replace this crude flag with
  • branches/numpy/anuga/abstract_2d_finite_volumes/quantity.py

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

    r6517 r6533  
    542542                        conserved_quantities =\
    543543                        ['stage', 'xmomentum', 'ymomentum'])
    544         domain.set_default_order(1)
    545544        domain.check_integrity()
    546545
     
    651650                        conserved_quantities =\
    652651                        ['stage', 'xmomentum', 'ymomentum'])
    653         domain.set_default_order(1)                       
    654652        domain.check_integrity()
    655653
     
    802800                             [ 11.0,  11.0,  11.0],
    803801                             [ 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        
    848802
    849803#-------------------------------------------------------------
  • branches/numpy/anuga/abstract_2d_finite_volumes/test_general_mesh.py

    r6517 r6533  
    9191                  str(verts)))
    9292        self.assert_(num.allclose(num.array([nodes_absolute[1],
    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        
     93                                             nodes_absolute[0],
     94                                             nodes_absolute[2]]),
     95                                  verts), msg)
    10296
    10397    def test_get_vertex_coordinates_triangle_id(self):
  • branches/numpy/anuga/abstract_2d_finite_volumes/util.py

    r6517 r6533  
    925925                    verbose = False):   
    926926       
    927     # FIXME(Ole): Shouldn't print statements here be governed by verbose?
    928927    assert type(gauge_filename) == type(''), 'Gauge filename must be a string'
    929928   
     
    972971            raise msg
    973972
    974         if verbose:
    975             print 'swwfile', swwfile
     973        print 'swwfile', swwfile
    976974
    977975        # Extract parent dir name and use as label
     
    24912489                                 base_name=base,
    24922490                                 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
    24992491   
    25002492    #to make all the quantities lower case for file_function
     
    25052497
    25062498    core_quantities = ['stage', 'elevation', 'xmomentum', 'ymomentum']
    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 
     2499   
    25242500    for sww_file in sww_files:
    25252501        sww_file = join(dir_name, sww_file+'.sww')
  • branches/numpy/anuga/alpha_shape/alpha_shape.py

    r6517 r6533  
    289289                       (denom[k]< EPSILON and  denom[k] > -EPSILON)]
    290290
    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            
     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
    297297        self.triradius = 0.5*num.sqrt(dx*dx + dy*dy)
    298298        #print "triangle radii", self.triradius
  • branches/numpy/anuga/caching/test_caching.py

    r6517 r6533  
    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
    390388        if os.name == 'posix' and os.uname()[4] in ['x86_64', 'ia64']:
    391389          # 64 bit hash values
    392           f1hash = 7079146893884768701
    393           f2hash = -6995306676314913340
    394 
    395           # Prior to cluster upgrades Feb 2009
    396           #f1hash = 1914027059797211698
    397           #f2hash = 1914027059807087171
     390          f1hash = 1914027059797211698
     391          f2hash = 1914027059807087171
    398392        else:
    399           # 32 bit hash values
    400393          f1hash = -758136387
    401394          f2hash = -11221564     
  • branches/numpy/anuga/compile_all.py

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

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

    r6517 r6533  
    3737    return sign*dd, mm, ss
    3838
    39 def redfearn(lat, lon, false_easting=None, false_northing=None,
    40              zone=None, central_meridian=None, scale_factor=None):
     39def redfearn(lat, lon, false_easting=None, false_northing=None, zone=None):
    4140    """Compute UTM projection using Redfearn's formula
    4241
     
    4847    If zone is specified reproject lat and long to specified zone instead of
    4948    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
    5449    """
    5550
     
    6257    a = 6378137.0                       #Semi major axis
    6358    inverse_flattening = 298.257222101  #1/f
    64     if scale_factor is None:
    65         K0 = 0.9996                         #Central scale factor   
    66     else:
    67         K0 = scale_factor
    68     #print 'scale', K0
     59    K0 = 0.9996                         #Central scale factor   
    6960    zone_width = 6                      #Degrees
    7061
     
    147138    m = term1 + term2 + term3 + term4 #OK
    148139
    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
     140    #Zone
    154141    if zone is None:
    155142        zone = int((lon - longitude_of_western_edge_zone0)/zone_width)
    156143
    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
     144    central_meridian = zone*zone_width+longitude_of_central_meridian_zone0
    162145
    163146    omega = (lon-central_meridian)*pi/180 #Relative longitude (radians)
  • branches/numpy/anuga/coordinate_transforms/test_geo_reference.py

    r6517 r6533  
    235235        y = 3.0
    236236       
    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            
     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
    250248        g = Geo_reference()
    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')
     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')
    261257           
    262258        # test that calling get_absolute twice does the right thing
  • branches/numpy/anuga/coordinate_transforms/test_redfearn.py

    r6517 r6533  
    122122
    123123    def test_UTM_6_nonstandard_projection(self):
    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)
     124        #Test 6 (Geraldton, WA)
     125
     126        #First test native projection (zone 50)
    134127        zone, easting, northing = redfearn(-29.233299999,114.05)
    135128
     
    274267    #    #assert allclose(northing, 6181725.1724276)
    275268
    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            
    380269    def test_convert_lats_longs(self):
    381270
  • branches/numpy/anuga/culvert_flows/Test_Culvert_Flat_Water_Lev.py

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

    r6517 r6533  
    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
    208206
    209207
     
    298296
    299297        # Print some diagnostics to log if requested
    300         if self.log_filename is not None:
     298        if hasattr(self, 'log_filename'):
    301299            s = 'Culvert Effective Length = %.2f m' %(self.length)
    302300            log_to_file(self.log_filename, s)
     
    326324                update = True
    327325
    328            
    329         if self.log_filename is not None:       
     326        if hasattr(self, 'log_filename'):           
    330327            s = '\nTime = %.2f, delta_t = %f' %(time, delta_t)
    331328            log_to_file(self.log_filename, s)
     
    442439            if self.verbose is True:
    443440                print msg
    444                
    445             if self.log_filename is not None:               
     441            if hasattr(self, 'log_filename'):                   
    446442                log_to_file(self.log_filename, msg)
    447443       
     
    528524                print '%.2fs - WARNING: Flow is running uphill.' % time
    529525           
    530         if self.log_filename is not None:
     526        if hasattr(self, 'log_filename'):
    531527            s = 'Time=%.2f, inlet stage = %.2f, outlet stage = %.2f'\
    532528                %(time, self.inlet.stage, self.outlet.stage)
     
    577573                    msg += 'for culvert "%s"' % self.label
    578574                   
    579                     if self.log_filename is not None:                   
     575                    if hasattr(self, 'log_filename'):                   
    580576                        log_to_file(self.log_filename, msg)
    581577            else:
    582578                # User culvert routine
    583579                Q, barrel_velocity, culvert_outlet_depth =\
    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)
     580                    self.culvert_routine(self, delta_total_energy, g)
    596581           
    597582           
     
    616601
    617602                           
    618 # OBSOLETE (Except for momentum jet in Culvert_flow_energy)   
     603   
    619604class Culvert_flow_rating:
    620605    """Culvert flow - transfer water from one hole to another
  • branches/numpy/anuga/culvert_flows/culvert_routines.py

    r6517 r6533  
    88
    99#NOTE:
    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)
     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 !!!!
    1437
    1538
     
    1740
    1841
    19 def 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
     42def 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
    4249    """
    4350
     
    4653    from anuga.utilities.numerical_tools import safe_acos as acos
    4754
    48     if inlet_depth > 0.01:
     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
     77            log_to_file(log_filename, s)
     78       
     79        msg = 'Specific energy is negative'
     80        assert E >= 0.0, msg
     81                     
     82       
    4983        # Water has risen above inlet
    50        
    51         if log_filename is not None:                           
    52             s = 'Specific energy  = %f m' % inlet_specific_energy
    53             log_to_file(log_filename, s)
    54        
    55         msg = 'Specific energy at inlet is negative'
    56         assert inlet_specific_energy >= 0.0, msg
    57                      
    58         if culvert_type == 'circle':
    59             # Round culvert (use width as diameter)
    60             diameter = culvert_width           
     84        if culvert.culvert_type == 'circle':
     85            # Round culvert
    6186
    6287            # Calculate flows for inlet control           
    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)
     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)
    6895                log_to_file(log_filename, s)
    6996
    70             # FIXME(Ole): Are these functions really for inlet control?                   
     97            case = ''
    7198            if Q_inlet_unsubmerged < Q_inlet_submerged:
    7299                Q = Q_inlet_unsubmerged
    73100
    74                 alpha = acos(1 - inlet_depth/diameter)
     101                alpha = acos(1 - inlet.depth/diameter)
    75102                flow_area = diameter**2 * (alpha - sin(alpha)*cos(alpha))
    76                 outlet_culvert_depth = inlet_depth
     103                outlet_culvert_depth = inlet.depth
     104                width = diameter*sin(alpha)
     105                #perimeter = alpha*diameter               
    77106                case = 'Inlet unsubmerged'
    78107            else:   
     
    80109                flow_area = (diameter/2)**2 * pi
    81110                outlet_culvert_depth = diameter
     111                width = diameter
     112                #perimeter = diameter
    82113                case = 'Inlet submerged'                   
    83114
    84115
    85116
    86             if delta_total_energy < inlet_specific_energy:
     117            if delta_total_energy < E:
    87118                # Calculate flows for outlet control
    88119                # Determine the depth at the outlet relative to the depth of flow in the Culvert
    89120
    90                 if outlet_depth > diameter:        # The Outlet is Submerged
     121                if outlet.depth > diameter:        # The Outlet is Submerged
    91122                    outlet_culvert_depth=diameter
    92123                    flow_area = (diameter/2)**2 * pi  # Cross sectional area of flow in the culvert
    93124                    perimeter = diameter * pi
     125                    width = diameter
    94126                    case = 'Outlet submerged'
    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)
     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)
    98130                    flow_area = diameter**2 * (alpha - sin(alpha)*cos(alpha))
    99131                    perimeter = alpha*diameter
     132                    width = diameter*sin(alpha)                   
    100133
    101134                    case = 'Outlet depth is zero'                       
    102135                else:   # Here really should use the Culvert Slope to calculate Actual Culvert Depth & Velocity
    103                     outlet_culvert_depth=outlet_depth   # For purpose of calculation assume the outlet depth = the inlet depth
    104                     alpha = acos(1 - outlet_depth/diameter)
     136                    outlet_culvert_depth=outlet.depth   # For purpose of calculation assume the outlet depth = the inlet depth
     137                    alpha = acos(1 - outlet.depth/diameter)
    105138                    flow_area = diameter**2 * (alpha - sin(alpha)*cos(alpha))
    106139                    perimeter = alpha*diameter
     140                    width = diameter*sin(alpha)                   
    107141                    case = 'Outlet is open channel flow'
    108142
    109143                hyd_rad = flow_area/perimeter
    110144               
    111                 if log_filename is not None:
     145                if hasattr(culvert, 'log_filename'):
    112146                    s = 'hydraulic radius at outlet = %f' %hyd_rad
    113147                    log_to_file(log_filename, s)
    114148
    115149                # Outlet control velocity using tail water
    116                 culvert_velocity = sqrt(delta_total_energy/((sum_loss/2*g)+(manning**2*culvert_length)/hyd_rad**1.33333))
     150                culvert_velocity = sqrt(delta_total_energy/((sum_loss/2*g)+(manning**2*length)/hyd_rad**1.33333))
    117151                Q_outlet_tailwater = flow_area * culvert_velocity
    118152               
    119                 if log_filename is not None:               
     153                if hasattr(culvert, 'log_filename'):
    120154                    s = 'Q_outlet_tailwater = %.6f' %Q_outlet_tailwater
    121155                    log_to_file(log_filename, s)
     
    131165
    132166            # Calculate flows for inlet control
    133             height = culvert_height
    134             width = culvert_width           
     167            height = culvert.height
     168            width = culvert.width           
    135169           
    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:                           
     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'):
    140174                s = 'Q_inlet_unsubmerged = %.6f, Q_inlet_submerged = %.6f' %(Q_inlet_unsubmerged, Q_inlet_submerged)
    141175                log_to_file(log_filename, s)
    142176
    143             # FIXME(Ole): Are these functions really for inlet control?   
     177            case = ''
    144178            if Q_inlet_unsubmerged < Q_inlet_submerged:
    145179                Q = Q_inlet_unsubmerged
    146                 flow_area = width*inlet_depth
    147                 outlet_culvert_depth = inlet_depth
     180                flow_area = width*inlet.depth
     181                outlet_culvert_depth = inlet.depth
     182                #perimeter=(width+2.0*inlet.depth)               
    148183                case = 'Inlet unsubmerged'
    149184            else:   
     
    154189                case = 'Inlet submerged'                   
    155190
    156             if delta_total_energy < inlet_specific_energy:
     191            if delta_total_energy < E:
    157192                # Calculate flows for outlet control
    158193                # Determine the depth at the outlet relative to the depth of flow in the Culvert
    159194
    160                 if outlet_depth > height:        # The Outlet is Submerged
     195                if outlet.depth > height:        # The Outlet is Submerged
    161196                    outlet_culvert_depth=height
    162197                    flow_area=width*height       # Cross sectional area of flow in the culvert
    163198                    perimeter=2.0*(width+height)
    164199                    case = 'Outlet submerged'
    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)
     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)
    169204                    case = 'Outlet depth is zero'                       
    170205                else:   # Here really should use the Culvert Slope to calculate Actual Culvert Depth & Velocity
    171                     outlet_culvert_depth=outlet_depth
    172                     flow_area=width*outlet_depth
    173                     perimeter=(width+2.0*outlet_depth)
     206                    outlet_culvert_depth=outlet.depth
     207                    flow_area=width*outlet.depth
     208                    perimeter=(width+2.0*outlet.depth)
    174209                    case = 'Outlet is open channel flow'
    175210
    176211                hyd_rad = flow_area/perimeter
    177212               
    178                 if log_filename is not None:                               
    179                     s = 'hydraulic radius at outlet = %f' % hyd_rad
     213                if hasattr(culvert, 'log_filename'):               
     214                    s = 'hydraulic radius at outlet = %f' %hyd_rad
    180215                    log_to_file(log_filename, s)
    181216
    182217                # Outlet control velocity using tail water
    183                 culvert_velocity = sqrt(delta_total_energy/((sum_loss/2*g)+(manning**2*culvert_length)/hyd_rad**1.33333))
     218                culvert_velocity = sqrt(delta_total_energy/((sum_loss/2*g)+(manning**2*length)/hyd_rad**1.33333))
    184219                Q_outlet_tailwater = flow_area * culvert_velocity
    185220
    186                 if log_filename is not None:                           
    187                     s = 'Q_outlet_tailwater = %.6f' % Q_outlet_tailwater
     221                if hasattr(culvert, 'log_filename'):                           
     222                    s = 'Q_outlet_tailwater = %.6f' %Q_outlet_tailwater
    188223                    log_to_file(log_filename, s)
    189224                Q = min(Q, Q_outlet_tailwater)
     
    193228
    194229        # Common code for circle and square geometries
    195         if log_filename is not None:
    196             log_to_file(log_filename, 'Case: "%s"' % case)
     230        if hasattr(culvert, 'log_filename'):                                   
     231            log_to_file(log_filename, 'Case: "%s"' %case)
    197232           
    198         if log_filename is not None:                       
    199             s = 'Flow Rate Control = %f' % Q
     233        flow_rate_control=Q
     234
     235        if hasattr(culvert, 'log_filename'):                                   
     236            s = 'Flow Rate Control = %f' %flow_rate_control
    200237            log_to_file(log_filename, s)
    201238
    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
     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
    206246            log_to_file(log_filename, s)
    207247
     
    210250
    211251
    212     else: # inlet_depth < 0.01:
     252    else: #inlet.depth < 0.01:
    213253        Q = barrel_velocity = outlet_culvert_depth = 0.0
    214254       
  • branches/numpy/anuga/fit_interpolate/fit.py

    r6517 r6533  
    310310                assert flag is True, msg               
    311311               
    312                 # FIXME(Ole): This is the message referred to in ticket:314
    313312                minx = min(self.mesh_boundary_polygon[:,0])
    314313                maxx = max(self.mesh_boundary_polygon[:,0])               
     
    318317                msg += 'Mesh boundary extent is (%.f, %.f), (%.f, %.f)'\
    319318                    % (minx, maxx, miny, maxy)
    320                 raise RuntimeError, msg
    321319               
    322320
     
    377375                self.build_fit_subset(points, z, verbose=verbose)
    378376
    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 
    387377               
    388378            point_coordinates = None
     
    392382        if point_coordinates is None:
    393383            if verbose: print 'Warning: no data points in fit'
    394             msg = 'No interpolation matrix.'
    395             assert self.AtA is not None, msg
    396             assert self.Atz is not None
     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
    397388           
    398389            # FIXME (DSG) - do  a message
  • branches/numpy/anuga/geospatial_data/geospatial_data.py

    r6517 r6533  
    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 
    18211808    if verbose: print "Number of points in sample to compare: ", len(points)
    18221809
  • branches/numpy/anuga/shallow_water/__init__.py

    r6517 r6533  
    1212     Transmissive_momentum_set_stage_boundary,\
    1313     Dirichlet_discharge_boundary,\
    14      Field_boundary,\
    15      Transmissive_stage_zero_momentum_boundary
     14     Field_boundary
    1615
    1716
  • branches/numpy/anuga/shallow_water/data_manager.py

    r6517 r6533  
    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):
    5866         if isinstance(times, (list, num.ndarray)):
     5865        if type(times) is list or isinstance(times, num.ndarray):
    58675866            number_of_times = len(times)
    58685867            times = ensure_numeric(times)
     
    59315930            #outfile.variables[q+Write_sww.RANGE][1] = -max_float # Max
    59325931
    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
     5932        if type(times) is list or isinstance(times, num.ndarray):
     5933            outfile.variables['time'][:] = times    #Store time relative
    59365934
    59375935        if verbose:
     
    62676265        # This is being used to seperate one number from a list.
    62686266        # what it is actually doing is sorting lists from numeric arrays.
    6269         #if type(times) is list or isinstance(times, num.ndarray):
    6270         if isinstance(times, (list, num.ndarray)):
     6267        if type(times) is list or isinstance(times, num.ndarray):
    62716268            number_of_times = len(times)
    62726269            times = ensure_numeric(times)
     
    63156312            outfile.variables[q + Write_sts.RANGE][1] = -max_float # Max
    63166313
    6317         #if type(times) is list or isinstance(times, num.ndarray):
    6318         if isinstance(times, (list, num.ndarray)):
     6314        if type(times) is list or isinstance(times, num.ndarray):
    63196315            outfile.variables['time'][:] = times    #Store time relative
    63206316
  • branches/numpy/anuga/shallow_water/test_shallow_water_domain.py

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

    r6517 r6533  
    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.