Changeset 6055


Ignore:
Timestamp:
Dec 10, 2008, 5:07:51 PM (16 years ago)
Author:
ole
Message:

Culvert work with Petar Milevski

Location:
anuga_core
Files:
5 edited

Legend:

Unmodified
Added
Removed
  • anuga_core/documentation/user_manual/version.tex

    r4788 r6055  
    77% release version; this is used to define the
    88% \version macro
    9 \release{1.0beta\_0000}
     9\release{1.0beta\_6051}
  • anuga_core/source/anuga/culvert_flows/Test_Culvert_Flat_Water_Lev.py

    r6051 r6055  
    6161    z=-x/1000
    6262   
    63         # Changing Slopes in the X- Direction
    64     # N = len(x)
    65     # for i in range(N):
    66     #     if 0 <x[i] < 5.1:
    67     #         z[i] -= -x[i]/10
    68     #     if 5 <x[i] < 10.1:
    69     #         z[i] -= -x[i]/100                             
    70     #     if 10 <x[i]:
    71     #         z[i] -= -x[i]/20                           
    72     # return z
    73 
    7463    #       NOW Add bits and Pieces to topography
    7564    N = len(x)
     
    9180                           
    9281       
    93         # Constriction
    94         #if 27 < x[i] < 29 and y[i] > 3:
    95         #    z[i] += 2       
    96        
    97         # Pole
    98         #if (x[i] - 34)**2 + (y[i] - 2)**2 < 0.4**2:
    99         #    z[i] += 2
    100 
    101         # HOLE For Pit at Opening[0]
    102         #if (x[i]-4)**2 + (y[i]-1.5)**2<0.75**2:
    103          #  z[i]-=1             
    104 
    105            # HOLE For Pit at Opening[1]
    106         #if (x[i]-20)**2 + (y[i]-3.5)**2<0.5**2:
    107          #  z[i]-=1             
    10882               
    10983    return z
     
    131105
    132106culvert = Culvert_flow(domain,
    133                        label='Culvert No. 1',
    134                        description='This culvert is a test unit 1.2m Wide by 0.75m High',   
     107                       culvert_description_filename='example_rating_curve.csv',
    135108                       end_point0=[9.0, 2.5],
    136109                       end_point1=[13.0, 2.5],
    137                        width=1.20,height=0.75,
    138                        culvert_routine=boyd_generalised_culvert_model,       
    139                        number_of_barrels=1,
    140                        update_interval=2,
    141110                       verbose=True)
     111
     112
     113#culvert = Culvert_flow(domain,
     114#                       label='Culvert No. 1',
     115#                       description='This culvert is a test unit 1.2m Wide by 0.75m High',   
     116#                       end_point0=[9.0, 2.5],
     117#                       end_point1=[13.0, 2.5],
     118#                       width=1.20,height=0.75,
     119#                       culvert_routine=boyd_generalised_culvert_model,       
     120#                       number_of_barrels=1,
     121#                       update_interval=2,
     122#                       verbose=True)
    142123
    143124domain.forcing_terms.append(culvert)
  • anuga_core/source/anuga/culvert_flows/culvert_class.py

    r6054 r6055  
    66from anuga.utilities.polygon import plot_polygons
    77
    8 
     8from anuga.utilities.numerical_tools import ensure_numeric
     9from Numeric import allclose
     10
     11import sys
     12
     13class Below_interval(Exception): pass
     14class Above_interval(Exception): pass
     15
     16
     17
     18# FIXME(Ole): Write in C and reuse this function by similar code in interpolate.py
     19def interpolate_linearly(x, xvec, yvec):
     20         
     21    # Find appropriate slot           
     22    i = 0
     23    while x > xvec[i]: i += 1
     24
     25    if i == 0:
     26        msg = 'Value provided = %.2f, interpolation minimum = %.2f.' %(x, xvec[0])
     27        raise Below_interval, msg
     28       
     29    if i == len(xvec):
     30        msg = 'Value provided = %.2f, interpolation maximum = %.2f.' %(x, xvec[-1])
     31        raise Above_interval, msg       
     32       
     33   
     34    x0 = xvec[i-1]
     35    x1 = xvec[i]           
     36    alpha = (x - x0)/(x1 - x0)
     37           
     38    y0 = yvec[i-1]
     39    y1 = yvec[i]                       
     40    y = alpha*y1 + (1-alpha)*y0
     41   
     42    return y
     43           
     44
     45           
     46def read_culvert_description(culvert_description_filename):           
     47   
     48    # Read description file
     49    fid = open(culvert_description_filename)
     50   
     51    read_rating_curve_data = False       
     52    rating_curve = []
     53    for i, line in enumerate(fid.readlines()):
     54       
     55        if read_rating_curve_data is True:
     56            fields = line.split(',')
     57            head_difference = float(fields[0].strip())
     58            flow_rate = float(fields[1].strip())               
     59            barrel_velocity = float(fields[2].strip())
     60           
     61            rating_curve.append( [head_difference, flow_rate, barrel_velocity] )
     62       
     63        if i == 0:
     64            # Header
     65            continue
     66        if i == 1:
     67            # Metadata
     68            fields = line.split(',')
     69            label=fields[0].strip()
     70            type=fields[1].strip().lower()
     71            assert type in ['box', 'pipe']
     72           
     73            width=float(fields[2].strip())
     74            height=float(fields[3].strip())
     75            length=float(fields[4].strip())
     76            number_of_barrels=int(fields[5].strip())
     77            #fields[6] refers to losses
     78            description=fields[7].strip()               
     79               
     80        if line.strip() == '': continue # Skip blanks
     81
     82        if line.startswith('Rating'):
     83            read_rating_curve_data = True
     84            # Flow data follows
     85               
     86    fid.close()
     87   
     88    return label, type, width, height, length, number_of_barrels, description, rating_curve
     89   
    990
    1091class Culvert_flow:
     
    21102    def __init__(self,
    22103                 domain,
    23                  
     104                 culvert_description_filename=None,
     105                 end_point0=None,
     106                 end_point1=None,
     107                 update_interval=None,
     108                 log_file=False,
     109                 discharge_hydrograph=False,
     110                 verbose=False):
     111       
     112        from Numeric import sqrt, sum
     113
     114       
     115        label, type, width, height, length, number_of_barrels, description, rating_curve = read_culvert_description(culvert_description_filename)
     116       
     117               
     118        # Store culvert information
     119        self.label = label
     120        self.description = description
     121        self.culvert_type = type
     122        self.rating_curve = ensure_numeric(rating_curve)
     123        self.number_of_barrels = number_of_barrels
     124
     125        if label is None: label = 'culvert_flow'
     126        label += '_' + str(id(self))
     127        self.label = label
     128       
     129        # File for storing discharge_hydrograph
     130        if discharge_hydrograph is True:
     131            self.timeseries_filename = label + '_timeseries.csv'
     132            fid = open(self.timeseries_filename, 'w')
     133            fid.write('time, discharge\n')
     134            fid.close()
     135
     136        # Log file for storing general textual output
     137        if log_file is True:       
     138            self.log_filename = label + '.log'         
     139            log_to_file(self.log_filename, self.label)       
     140            log_to_file(self.log_filename, description)
     141            log_to_file(self.log_filename, self.culvert_type)       
     142
     143
     144        # Create the fundamental culvert polygons from POLYGON
     145        #if self.culvert_type == 'circle':
     146        #    # Redefine width and height for use with create_culvert_polygons
     147        #    width = height = diameter
     148       
     149        P = create_culvert_polygons(end_point0,
     150                                    end_point1,
     151                                    width=width,   
     152                                    height=height,
     153                                    number_of_barrels=number_of_barrels)
     154       
     155        if verbose is True:
     156            pass
     157            #plot_polygons([[end_point0, end_point1],
     158            #               P['exchange_polygon0'],
     159            #               P['exchange_polygon1'],
     160            #               P['enquiry_polygon0'],
     161            #               P['enquiry_polygon1']],
     162            #              figname='culvert_polygon_output')
     163
     164                         
     165        # Compute the average point for enquiry
     166        enquiry_point0 = sum(P['enquiry_polygon0'][:2])/2
     167        enquiry_point1 = sum(P['enquiry_polygon1'][:2])/2         
     168       
     169        self.enquiry_points = [enquiry_point0, enquiry_point1]                           
     170        self.enquiry_indices = []                 
     171        for point in self.enquiry_points:
     172            # Find nearest centroid
     173            N = len(domain)   
     174            points = domain.get_centroid_coordinates(absolute=True)
     175
     176            # Calculate indices in exchange area for this forcing term
     177           
     178            triangle_id = min_dist = sys.maxint
     179            for k in range(N):
     180                x, y = points[k,:] # Centroid
     181
     182                c = point                               
     183                distance = (x-c[0])**2+(y-c[1])**2
     184                if distance < min_dist:
     185                    min_dist = distance
     186                    triangle_id = k
     187
     188                   
     189            if triangle_id < sys.maxint:
     190                msg = 'found triangle with centroid (%f, %f)'\
     191                    %tuple(points[triangle_id, :])
     192                msg += ' for point (%f, %f)' %tuple(point)
     193               
     194                self.enquiry_indices.append(triangle_id)
     195            else:
     196                msg = 'Triangle not found for point (%f, %f)' %point
     197                raise Exception, msg
     198       
     199                         
     200
     201        # Check that all polygons lie within the mesh.
     202        bounding_polygon = domain.get_boundary_polygon()
     203        for key in P.keys():
     204            if key in ['exchange_polygon0',
     205                       'exchange_polygon1']:
     206                for point in list(P[key]) + self.enquiry_points:
     207                    msg = 'Point %s in polygon %s for culvert %s did not'\
     208                        %(str(point), key, self.label)
     209                    msg += 'fall within the domain boundary.'
     210                    assert is_inside_polygon(point, bounding_polygon), msg
     211       
     212                   
     213        # Create inflow object at each end of the culvert.
     214        self.openings = []
     215        self.openings.append(Inflow(domain,
     216                                    polygon=P['exchange_polygon0']))
     217
     218        self.openings.append(Inflow(domain,
     219                                    polygon=P['exchange_polygon1']))                                   
     220
     221
     222
     223        dq = domain.quantities                                           
     224        for i, opening in enumerate(self.openings):                           
     225            #elevation = dq['elevation'].get_values(location='centroids',
     226            #                                      interpolation_points=[self.enquiry_points[i]])
     227           
     228            elevation = dq['elevation'].get_values(location='centroids',
     229                                                   indices=[self.enquiry_indices[i]])           
     230            opening.elevation = elevation
     231
     232        # Assume two openings for now: Referred to as 0 and 1
     233        assert len(self.openings) == 2
     234                           
     235        # Determine pipe direction     
     236        self.delta_z = delta_z = self.openings[0].elevation - self.openings[1].elevation
     237        if delta_z > 0.0:
     238            self.inlet = self.openings[0]
     239            self.outlet = self.openings[1]
     240        else:               
     241            self.outlet = self.openings[0]
     242            self.inlet = self.openings[1]       
     243       
     244       
     245        # Store basic geometry
     246        self.end_points = [end_point0, end_point1]
     247        self.vector = P['vector']
     248        self.length = P['length']; assert self.length > 0.0
     249        if not allclose(self.length, length, rtol=1.0e-2, atol=1.0e-2):
     250            msg = 'WARNING: barrel length specified in "%s" (%.2f m)' %(culvert_description_filename, length)
     251            msg += ' does not match distance between specified'
     252            msg += ' end points (%.2f m)' %self.length
     253            print msg
     254       
     255        self.verbose = verbose
     256        self.last_update = 0.0 # For use with update_interval       
     257        self.update_interval = update_interval
     258       
     259
     260        # Print something
     261        if hasattr(self, 'log_filename'):
     262            s = 'Culvert Effective Length = %.2f m' %(self.length)
     263            log_to_file(self.log_filename, s)
     264   
     265            s = 'Culvert Direction is %s\n' %str(self.vector)
     266            log_to_file(self.log_filename, s)       
     267       
     268       
     269       
     270       
     271       
     272    def __call__(self, domain):
     273        from anuga.utilities.numerical_tools import mean
     274       
     275        from anuga.config import g, epsilon
     276        from Numeric import take, sqrt
     277        from anuga.config import velocity_protection       
     278
     279
     280        # Time stuff
     281        time = domain.get_time()
     282       
     283       
     284        update = False
     285        if self.update_interval is None:
     286            update = True
     287        else:   
     288            if time - self.last_update > self.update_interval or time == 0.0:
     289                update = True
     290
     291                               
     292        if update is True:
     293            self.last_update = time
     294            dq = domain.quantities
     295                       
     296            # Get average water depths at each opening       
     297            openings = self.openings   # There are two Opening [0] and [1]
     298            for i, opening in enumerate(openings):
     299               
     300                # Compute mean values of selected quantitites in the
     301                # enquiry area in front of the culvert
     302                # Stage and velocity comes from enquiry area
     303                # and elevation from exchange area
     304               
     305                stage = dq['stage'].get_values(location='centroids',
     306                                               indices=[self.enquiry_indices[i]])
     307
     308                   
     309                # Store current average stage and depth with each opening object
     310                opening.stage = stage
     311
     312               
     313
     314            #################  End of the FOR loop ################################################
     315
     316            # We now need to deal with each opening individually
     317               
     318            # Determine flow direction based on total energy difference
     319
     320            delta_w = self.inlet.stage - self.outlet.stage
     321           
     322            if hasattr(self, 'log_filename'):
     323                s = 'Time=%.2f, inlet stage = %.2f, outlet stage = %.2f' %(time,
     324                                                                           self.inlet.stage,
     325                                                                           self.outlet.stage)
     326                log_to_file(self.log_filename, s)
     327
     328
     329            # Calculate discharge for one barrel and set inlet.rate and outlet.rate
     330            try:
     331                Q = interpolate_linearly(delta_w, self.rating_curve[:,0], self.rating_curve[:,1])
     332            except Below_interval, e:
     333                Q = self.rating_curve[0,1]             
     334                msg = '%.2fs: Delta head smaller than rating curve minimum: ' %time
     335                msg += str(e)
     336                msg += '\n        I will use minimum discharge %.2f m^3/s for culvert "%s"'\
     337                    %(Q, self.label)
     338                if hasattr(self, 'log_filename'):                   
     339                    log_to_file(self.log_filename, msg)
     340            except Above_interval, e:
     341                Q = self.rating_curve[-1,1]             
     342                msg = '%.2fs: Delta head greater than rating curve maximum: ' %time
     343                msg += str(e)
     344                msg += '\n        I will use maximum discharge %.2f m^3/s for culvert "%s"'\
     345                    %(Q, self.label)
     346                if hasattr(self, 'log_filename'):                   
     347                    log_to_file(self.log_filename, msg)                   
     348
     349               
     350               
     351           
     352            # Adjust discharge for multiple barrels
     353            Q *= self.number_of_barrels
     354           
     355           
     356            self.inlet.rate = -Q
     357            self.outlet.rate = Q
     358
     359            # Log timeseries to file
     360            try:
     361                fid = open(self.timeseries_filename, 'a')       
     362            except:
     363                pass
     364            else:   
     365                fid.write('%.2f, %.2f\n' %(time, Q))
     366                fid.close()
     367
     368
     369
     370           
     371   
     372        # Execute flow term for each opening
     373        # This is where Inflow objects are evaluated using the last rate that has been calculated
     374        #
     375        # This will take place at every internal timestep and update the domain
     376        for opening in self.openings:
     377            opening(domain)
     378
     379
     380
     381       
     382       
     383
     384class Culvert_flow_energy:
     385    """Culvert flow - transfer water from one hole to another
     386   
     387    Using Momentum as Calculated by Culvert Flow !!
     388    Could be Several Methods Investigated to do This !!!
     389
     390    2008_May_08
     391    To Ole:
     392    OK so here we need to get the Polygon Creating code to create
     393    polygons for the culvert Based on
     394    the two input Points (X0,Y0) and (X1,Y1) - need to be passed
     395    to create polygon
     396
     397    The two centers are now passed on to create_polygon.
     398   
     399
     400    Input: Two points, pipe_size (either diameter or width, height),
     401    mannings_rougness,
     402    inlet/outlet energy_loss_coefficients, internal_bend_coefficent,
     403    top-down_blockage_factor and bottom_up_blockage_factor
     404   
     405   
     406    And the Delta H enquiry should be change from Openings in line 412
     407    to the enquiry Polygons infront of the culvert
     408    At the moment this script uses only Depth, later we can change it to
     409    Energy...
     410
     411    Once we have Delta H can calculate a Flow Rate and from Flow Rate
     412    an Outlet Velocity
     413    The Outlet Velocity x Outlet Depth = Momentum to be applied at the Outlet...
     414       
     415    Invert levels are optional. If left out they will default to the
     416    elevation at the opening.
     417       
     418    """
     419
     420    def __init__(self,
     421                 domain,
    24422                 label=None,
    25423                 description=None,
    26424                 end_point0=None,
    27425                 end_point1=None,
    28                  
    29                  
    30 
     426                 width=None,
     427                 height=None,
     428                 diameter=None,
    31429                 manning=None,          # Mannings Roughness for Culvert
    32430                 invert_level0=None,    # Invert level at opening 0
     
    428826
    429827
    430 
    431        
    432        
    433 
    434 class Culvert_flow_energy:
    435     """Culvert flow - transfer water from one hole to another
    436    
    437     Using Momentum as Calculated by Culvert Flow !!
    438     Could be Several Methods Investigated to do This !!!
    439 
    440     2008_May_08
    441     To Ole:
    442     OK so here we need to get the Polygon Creating code to create
    443     polygons for the culvert Based on
    444     the two input Points (X0,Y0) and (X1,Y1) - need to be passed
    445     to create polygon
    446 
    447     The two centers are now passed on to create_polygon.
    448    
    449 
    450     Input: Two points, pipe_size (either diameter or width, height),
    451     mannings_rougness,
    452     inlet/outlet energy_loss_coefficients, internal_bend_coefficent,
    453     top-down_blockage_factor and bottom_up_blockage_factor
    454    
    455    
    456     And the Delta H enquiry should be change from Openings in line 412
    457     to the enquiry Polygons infront of the culvert
    458     At the moment this script uses only Depth, later we can change it to
    459     Energy...
    460 
    461     Once we have Delta H can calculate a Flow Rate and from Flow Rate
    462     an Outlet Velocity
    463     The Outlet Velocity x Outlet Depth = Momentum to be applied at the Outlet...
    464        
    465     Invert levels are optional. If left out they will default to the
    466     elevation at the opening.
    467        
    468     """
    469 
    470     def __init__(self,
    471                  domain,
    472                  label=None,
    473                  description=None,
    474                  end_point0=None,
    475                  end_point1=None,
    476                  width=None,
    477                  height=None,
    478                  diameter=None,
    479                  manning=None,          # Mannings Roughness for Culvert
    480                  invert_level0=None,    # Invert level at opening 0
    481                  invert_level1=None,    # Invert level at opening 1
    482                  loss_exit=None,
    483                  loss_entry=None,
    484                  loss_bend=None,
    485                  loss_special=None,
    486                  blockage_topdwn=None,
    487                  blockage_bottup=None,
    488                  culvert_routine=None,
    489                  number_of_barrels=1,
    490                  update_interval=None,
    491                  verbose=False):
    492        
    493         from Numeric import sqrt, sum
    494 
    495         # Input check
    496         if diameter is not None:
    497             self.culvert_type = 'circle'
    498             self.diameter = diameter
    499             if height is not None or width is not None:
    500                 msg = 'Either diameter or width&height must be specified, '
    501                 msg += 'but not both.'
    502                 raise Exception, msg
    503         else:
    504             if height is not None:
    505                 if width is None:
    506                     self.culvert_type = 'square'                               
    507                     width = height
    508                 else:
    509                     self.culvert_type = 'rectangle'
    510             elif width is not None:
    511                 if height is None:
    512                     self.culvert_type = 'square'                               
    513                     height = width
    514             else:
    515                 msg = 'Either diameter or width&height must be specified.'
    516                 raise Exception, msg               
    517                
    518             if height == width:
    519                 self.culvert_type = 'square'                                               
    520                
    521             self.height = height
    522             self.width = width
    523 
    524            
    525         assert self.culvert_type in ['circle', 'square', 'rectangle']
    526        
    527         assert number_of_barrels >= 1
    528         self.number_of_barrels = number_of_barrels
    529        
    530        
    531         # Set defaults
    532         if manning is None: manning = 0.012   # Default roughness for pipe
    533         if loss_exit is None: loss_exit = 1.00
    534         if loss_entry is None: loss_entry = 0.50
    535         if loss_bend is None: loss_bend=0.00
    536         if loss_special is None: loss_special=0.00
    537         if blockage_topdwn is None: blockage_topdwn=0.00
    538         if blockage_bottup is None: blockage_bottup=0.00
    539         if culvert_routine is None:
    540             culvert_routine=boyd_generalised_culvert_model
    541            
    542         if label is None: label = 'culvert_flow'
    543         label += '_' + str(id(self))
    544         self.label = label
    545        
    546         # File for storing culvert quantities
    547         self.timeseries_filename = label + '_timeseries.csv'
    548         fid = open(self.timeseries_filename, 'w')
    549         fid.write('time, E0, E1, Velocity, Discharge\n')
    550         fid.close()
    551 
    552         # Log file for storing general textual output
    553         self.log_filename = label + '.log'         
    554         log_to_file(self.log_filename, self.label)       
    555         log_to_file(self.log_filename, description)
    556         log_to_file(self.log_filename, self.culvert_type)       
    557 
    558 
    559         # Create the fundamental culvert polygons from POLYGON
    560         if self.culvert_type == 'circle':
    561             # Redefine width and height for use with create_culvert_polygons
    562             width = height = diameter
    563        
    564         P = create_culvert_polygons(end_point0,
    565                                     end_point1,
    566                                     width=width,   
    567                                     height=height,
    568                                     number_of_barrels=number_of_barrels)
    569        
    570         if verbose is True:
    571             pass
    572             #plot_polygons([[end_point0, end_point1],
    573             #               P['exchange_polygon0'],
    574             #               P['exchange_polygon1'],
    575             #               P['enquiry_polygon0'],
    576             #               P['enquiry_polygon1']],
    577             #              figname='culvert_polygon_output')
    578             #import sys; sys.exit()                           
    579 
    580 
    581         # Check that all polygons lie within the mesh.
    582         bounding_polygon = domain.get_boundary_polygon()
    583         for key in P.keys():
    584             if key in ['exchange_polygon0',
    585                        'exchange_polygon1',
    586                        'enquiry_polygon0',
    587                        'enquiry_polygon1']:
    588                 for point in P[key]:
    589                
    590                     msg = 'Point %s in polygon %s for culvert %s did not'\
    591                         %(str(point), key, self.label)
    592                     msg += 'fall within the domain boundary.'
    593                     assert is_inside_polygon(point, bounding_polygon), msg
    594        
    595 
    596         # Create inflow object at each end of the culvert.
    597         self.openings = []
    598         self.openings.append(Inflow(domain,
    599                                     polygon=P['exchange_polygon0']))
    600 
    601         self.openings.append(Inflow(domain,
    602                                     polygon=P['exchange_polygon1']))                                   
    603 
    604 
    605         # Assume two openings for now: Referred to as 0 and 1
    606         assert len(self.openings) == 2
    607        
    608         # Store basic geometry
    609         self.end_points = [end_point0, end_point1]
    610         self.invert_levels = [invert_level0, invert_level1]               
    611         #self.enquiry_polygons = [P['enquiry_polygon0'], P['enquiry_polygon1']]
    612         self.enquiry_polylines = [P['enquiry_polygon0'][:2],
    613                                   P['enquiry_polygon1'][:2]]
    614         self.vector = P['vector']
    615         self.length = P['length']; assert self.length > 0.0
    616         self.verbose = verbose
    617         self.last_time = 0.0       
    618         self.last_update = 0.0 # For use with update_interval       
    619         self.update_interval = update_interval
    620        
    621 
    622         # Store hydraulic parameters
    623         self.manning = manning
    624         self.loss_exit = loss_exit
    625         self.loss_entry = loss_entry
    626         self.loss_bend = loss_bend
    627         self.loss_special = loss_special
    628         self.sum_loss = loss_exit + loss_entry + loss_bend + loss_special
    629         self.blockage_topdwn = blockage_topdwn
    630         self.blockage_bottup = blockage_bottup
    631 
    632         # Store culvert routine
    633         self.culvert_routine = culvert_routine
    634 
    635        
    636         # Create objects to update momentum (a bit crude at this stage)
    637 
    638        
    639         xmom0 = General_forcing(domain, 'xmomentum',
    640                                 polygon=P['exchange_polygon0'])
    641 
    642         xmom1 = General_forcing(domain, 'xmomentum',
    643                                 polygon=P['exchange_polygon1'])
    644 
    645         ymom0 = General_forcing(domain, 'ymomentum',
    646                                 polygon=P['exchange_polygon0'])
    647 
    648         ymom1 = General_forcing(domain, 'ymomentum',
    649                                 polygon=P['exchange_polygon1'])
    650 
    651         self.opening_momentum = [ [xmom0, ymom0], [xmom1, ymom1] ]
    652        
    653 
    654         # Print something
    655         s = 'Culvert Effective Length = %.2f m' %(self.length)
    656         log_to_file(self.log_filename, s)
    657    
    658         s = 'Culvert Direction is %s\n' %str(self.vector)
    659         log_to_file(self.log_filename, s)       
    660        
    661        
    662     def __call__(self, domain):
    663         from anuga.utilities.numerical_tools import mean
    664        
    665         from anuga.config import g, epsilon
    666         from Numeric import take, sqrt
    667         from anuga.config import velocity_protection       
    668 
    669 
    670         log_filename = self.log_filename
    671          
    672         # Time stuff
    673         time = domain.get_time()
    674        
    675        
    676         update = False
    677         if self.update_interval is None:
    678             update = True
    679         else:   
    680             if time - self.last_update > self.update_interval or time == 0.0:
    681                 update = True
    682 
    683         #print 'call', time, time - self.last_update
    684                
    685                                
    686         if update is True:
    687             #print 'Updating', time, time - self.last_update
    688             self.last_update = time
    689        
    690             delta_t = time-self.last_time
    691             s = '\nTime = %.2f, delta_t = %f' %(time, delta_t)
    692             log_to_file(log_filename, s)
    693        
    694             msg = 'Time did not advance'
    695             if time > 0.0: assert delta_t > 0.0, msg
    696 
    697 
    698             # Get average water depths at each opening       
    699             openings = self.openings   # There are two Opening [0] and [1]
    700             for i, opening in enumerate(openings):
    701                 dq = domain.quantities
    702                
    703                 # Compute mean values of selected quantitites in the
    704                 # exchange area in front of the culvert
    705                 # Stage and velocity comes from enquiry area
    706                 # and elevation from exchange area
    707                
    708                 stage = dq['stage'].get_values(location='centroids',
    709                                                indices=opening.exchange_indices)           
    710                 w = mean(stage) # Average stage
    711 
    712                 # Use invert level instead of elevation if specified
    713                 invert_level = self.invert_levels[i]
    714                 if invert_level is not None:
    715                     z = invert_level
    716                 else:
    717                     elevation = dq['elevation'].get_values(location='centroids',
    718                                                            indices=opening.exchange_indices)
    719                     z = mean(elevation) # Average elevation
    720 
    721                 # Estimated depth above the culvert inlet
    722                 d = w - z  # Used for calculations involving depth
    723                 if d < 0.0:
    724                     # This is possible since w and z are taken at different locations
    725                     #msg = 'D < 0.0: %f' %d
    726                     #raise msg
    727                     d = 0.0
    728                
    729 
    730                 # Ratio of depth to culvert height.
    731                 # If ratio > 1 then culvert is running full
    732                 if self.culvert_type == 'circle':
    733                     ratio = d/self.diameter
    734                 else:   
    735                     ratio = d/self.height 
    736                 opening.ratio = ratio
    737                    
    738                    
    739                 # Average measures of energy in front of this opening
    740                 polyline = self.enquiry_polylines[i]
    741                 #print 't = %.4f, opening=%d,' %(domain.time, i),
    742                 opening.total_energy = domain.get_energy_through_cross_section(polyline,
    743                                                                                kind='total')           
    744                 #print 'Et = %.3f m' %opening.total_energy
    745 
    746                 # Store current average stage and depth with each opening object
    747                 opening.depth = d
    748                 opening.depth_trigger = d # Use this for now
    749                 opening.stage = w
    750                 opening.elevation = z
    751                
    752 
    753             #################  End of the FOR loop ################################################
    754 
    755             # We now need to deal with each opening individually
    756                
    757             # Determine flow direction based on total energy difference
    758             delta_Et = openings[0].total_energy - openings[1].total_energy
    759 
    760             if delta_Et > 0:
    761                 #print 'Flow U/S ---> D/S'
    762                 inlet = openings[0]
    763                 outlet = openings[1]
    764 
    765                 inlet.momentum = self.opening_momentum[0]
    766                 outlet.momentum = self.opening_momentum[1]
    767 
    768             else:
    769                 #print 'Flow D/S ---> U/S'
    770                 inlet = openings[1]
    771                 outlet = openings[0]
    772 
    773                 inlet.momentum = self.opening_momentum[1]
    774                 outlet.momentum = self.opening_momentum[0]
    775                
    776                 delta_Et = -delta_Et
    777 
    778             self.inlet = inlet
    779             self.outlet = outlet
    780                
    781             msg = 'Total energy difference is negative'
    782             assert delta_Et > 0.0, msg
    783 
    784             delta_h = inlet.stage - outlet.stage
    785             delta_z = inlet.elevation - outlet.elevation
    786             culvert_slope = (delta_z/self.length)
    787 
    788             if culvert_slope < 0.0:
    789                 # Adverse gradient - flow is running uphill
    790                 # Flow will be purely controlled by uphill outlet face
    791                 if self.verbose is True:
    792                     print 'WARNING: Flow is running uphill. Watch Out!', inlet.elevation, outlet.elevation
    793 
    794 
    795             s = 'Delta total energy = %.3f' %(delta_Et)
    796             log_to_file(log_filename, s)
    797 
    798 
    799             # Calculate discharge for one barrel and set inlet.rate and outlet.rate
    800             Q, barrel_velocity, culvert_outlet_depth = self.culvert_routine(self, inlet, outlet, delta_Et, g)
    801        
    802             # Adjust discharge for multiple barrels
    803             Q *= self.number_of_barrels
    804 
    805             # Compute barrel momentum
    806             barrel_momentum = barrel_velocity*culvert_outlet_depth
    807                    
    808             s = 'Barrel velocity = %f' %barrel_velocity
    809             log_to_file(log_filename, s)
    810 
    811             # Compute momentum vector at outlet
    812             outlet_mom_x, outlet_mom_y = self.vector * barrel_momentum
    813                
    814             s = 'Directional momentum = (%f, %f)' %(outlet_mom_x, outlet_mom_y)
    815             log_to_file(log_filename, s)
    816 
    817             # Log timeseries to file
    818             fid = open(self.timeseries_filename, 'a')       
    819             fid.write('%f, %f, %f, %f, %f\n'\
    820                           %(time,
    821                             openings[0].total_energy,
    822                             openings[1].total_energy,
    823                             barrel_velocity,
    824                             Q))
    825             fid.close()
    826 
    827             # Update momentum       
    828             delta_t = time - self.last_time
    829             if delta_t > 0.0:
    830                 xmomentum_rate = outlet_mom_x - outlet.momentum[0].value
    831                 xmomentum_rate /= delta_t
    832                    
    833                 ymomentum_rate = outlet_mom_y - outlet.momentum[1].value
    834                 ymomentum_rate /= delta_t
    835                        
    836                 s = 'X Y MOM_RATE = (%f, %f) ' %(xmomentum_rate, ymomentum_rate)
    837                 log_to_file(log_filename, s)                   
    838             else:
    839                 xmomentum_rate = ymomentum_rate = 0.0
    840 
    841 
    842             # Set momentum rates for outlet jet
    843             outlet.momentum[0].rate = xmomentum_rate
    844             outlet.momentum[1].rate = ymomentum_rate
    845 
    846             # Remember this value for next step (IMPORTANT)
    847             outlet.momentum[0].value = outlet_mom_x
    848             outlet.momentum[1].value = outlet_mom_y                   
    849 
    850             if int(domain.time*100) % 100 == 0:
    851                 s = 'T=%.5f, Culvert Discharge = %.3f f'\
    852                     %(time, Q)
    853                 s += ' Depth= %0.3f  Momentum = (%0.3f, %0.3f)'\
    854                      %(culvert_outlet_depth, outlet_mom_x,outlet_mom_y)
    855                 s += ' Momentum rate: (%.4f, %.4f)'\
    856                      %(xmomentum_rate, ymomentum_rate)                   
    857                 s+='Outlet Vel= %.3f'\
    858                     %(barrel_velocity)
    859                 log_to_file(log_filename, s)
    860            
    861                
    862 
    863 
    864         # Execute flow term for each opening
    865         # This is where Inflow objects are evaluated and update the domain
    866         for opening in self.openings:
    867             opening(domain)
    868 
    869         # Execute momentum terms
    870         # This is where Inflow objects are evaluated and update the domain
    871         self.outlet.momentum[0](domain)
    872         self.outlet.momentum[1](domain)       
    873            
    874         # Store value of time #FIXME(Ole): Maybe only every time we update   
    875         self.last_time = time
    876 
    877 
    878        
     828       
  • anuga_core/source/anuga/culvert_flows/example_rating_curve.csv

    r6054 r6055  
    1 Name, type, width or diameter, height, length, losses, description     
    2 Test Culvert, box, 3, 1.8, 5, 1, 50% blocked test case
     1Name, type, width or diameter, height, length, number of barrels, losses, description   
     2Test Culvert, box, 3, 1.8, 4, 1, 1, 50% blocked single barrel test case
    33
    44
    5 Rating Curve
    6 Delta W (m), Q (m3/s), Velocity (m/s)
     5Rating Curve: Delta W (m), Q (m3/s), Velocity (m/s)
    760,0,0
    870.1,0.720243203,0.800270226
  • anuga_core/source/anuga/shallow_water/shallow_water_domain.py

    r6045 r6055  
    19821982    #------------------------------------------------------------------------
    19831983    # This is the new element implemented by Ole and Rudy to allow direct
    1984     # input of Inflow in mm/s
     1984    # input of Rainfall in mm/s
    19851985
    19861986    catchmentrainfall = Rainfall(rain=file_function('Q100_2hr_Rain.tms')) 
Note: See TracChangeset for help on using the changeset viewer.