Changeset 6123


Ignore:
Timestamp:
Jan 8, 2009, 4:54:41 PM (10 years ago)
Author:
ole
Message:

Implented general inlet/outlet control.
Works with rating curve but not Boyd.

Location:
anuga_core/source/anuga/culvert_flows
Files:
3 edited

Legend:

Unmodified
Added
Removed
  • anuga_core/source/anuga/culvert_flows/culvert_class.py

    r6121 r6123  
    1111from anuga.config import g, epsilon
    1212from Numeric import take, sqrt
    13 from anuga.config import velocity_protection       
     13from anuga.config import minimum_allowed_height, velocity_protection       
    1414
    1515       
     
    114114   
    115115
     116   
     117
     118class Culvert_flow_general:
     119    """Culvert flow - transfer water from one hole to another
     120   
     121    This version will work with either rating curve file or with culvert
     122    routine.
     123   
     124    Input: Two points, pipe_size (either diameter or width, height),
     125    mannings_rougness,
     126    """
     127
     128    def __init__(self,
     129                 domain,
     130                 culvert_description_filename=None,
     131                 culvert_routine=None,
     132                 end_point0=None,
     133                 end_point1=None,
     134                 enquiry_point0=None,
     135                 enquiry_point1=None,
     136                 type='box',
     137                 width=None,
     138                 height=None,
     139                 length=None,
     140                 number_of_barrels=None,
     141                 trigger_depth=0.01, # Depth below which no flow happens
     142                 manning=None,          # Mannings Roughness for Culvert
     143                 sum_loss=None,
     144                 label=None,
     145                 description=None,
     146                 update_interval=None,
     147                 log_file=False,
     148                 discharge_hydrograph=False,
     149                 verbose=False):
     150       
     151
     152        self.culvert_routine = culvert_routine       
     153        self.culvert_description_filename = culvert_description_filename
     154        if culvert_description_filename is not None:
     155            label, type, width, height, length, number_of_barrels, description, rating_curve = read_culvert_description(culvert_description_filename)
     156            self.rating_curve = ensure_numeric(rating_curve)           
     157       
     158        self.domain = domain
     159        self.trigger_depth = trigger_depth       
     160               
     161        if manning is None:
     162            self.manning = 0.012   # Default roughness for pipe
     163       
     164        if sum_loss is None:
     165            self.sum_loss = 0.0
     166           
     167                       
     168        # Store culvert information
     169        self.label = label
     170        self.description = description
     171        self.culvert_type = type
     172        self.number_of_barrels = number_of_barrels
     173
     174        if label is None: label = 'culvert_flow'
     175        label += '_' + str(id(self))
     176        self.label = label
     177       
     178        # File for storing discharge_hydrograph
     179        if discharge_hydrograph is True:
     180            self.timeseries_filename = label + '_timeseries.csv'
     181            fid = open(self.timeseries_filename, 'w')
     182            fid.write('time, discharge\n')
     183            fid.close()
     184
     185        # Log file for storing general textual output
     186        if log_file is True:       
     187            self.log_filename = label + '.log'         
     188            log_to_file(self.log_filename, self.label)       
     189            log_to_file(self.log_filename, description)
     190            log_to_file(self.log_filename, self.culvert_type)       
     191
     192
     193        # Create the fundamental culvert polygons from polygon
     194        P = create_culvert_polygons(end_point0,
     195                                    end_point1,
     196                                    width=width,   
     197                                    height=height,
     198                                    number_of_barrels=number_of_barrels)
     199        self.culvert_polygons = P
     200       
     201        # Select enquiry points
     202        if enquiry_point0 is None:
     203            enquiry_point0 = P['enquiry_point0']
     204           
     205        if enquiry_point1 is None:
     206            enquiry_point1 = P['enquiry_point1']           
     207           
     208        if verbose is True:
     209            pass
     210            #plot_polygons([[end_point0, end_point1],
     211            #               P['exchange_polygon0'],
     212            #               P['exchange_polygon1'],
     213            #               [enquiry_point0, 1.005*enquiry_point0],
     214            #               [enquiry_point1, 1.005*enquiry_point1]],
     215            #              figname='culvert_polygon_output')
     216
     217           
     218           
     219        self.enquiry_points = [enquiry_point0, enquiry_point1]
     220        self.enquiry_indices = self.get_enquiry_indices()                 
     221        self.check_culvert_inside_domain()
     222       
     223                   
     224        # Create inflow object at each end of the culvert.
     225        self.openings = []
     226        self.openings.append(Inflow(domain,
     227                                    polygon=P['exchange_polygon0']))
     228        self.openings.append(Inflow(domain,
     229                                    polygon=P['exchange_polygon1']))
     230                                   
     231        # Assume two openings for now: Referred to as 0 and 1
     232        assert len(self.openings) == 2
     233
     234        # Establish initial values at each enquiry point
     235        dq = domain.quantities
     236        for i, opening in enumerate(self.openings):
     237            idx = self.enquiry_indices[i]
     238            elevation = dq['elevation'].get_values(location='centroids',
     239                                                   indices=[idx])[0]
     240            stage = dq['stage'].get_values(location='centroids',
     241                                           indices=[idx])[0]
     242            opening.elevation = elevation
     243            opening.stage = stage
     244            opening.depth = stage-elevation
     245
     246           
     247                           
     248        # Determine initial pipe direction.
     249        # This may change dynamically based on the total energy difference     
     250        # Consequently, this may be superfluous
     251        delta_z = self.openings[0].elevation - self.openings[1].elevation
     252        if delta_z > 0.0:
     253            self.inlet = self.openings[0]
     254            self.outlet = self.openings[1]
     255        else:               
     256            self.outlet = self.openings[0]
     257            self.inlet = self.openings[1]       
     258       
     259       
     260        # Store basic geometry
     261        self.end_points = [end_point0, end_point1]
     262        self.vector = P['vector']
     263        self.length = P['length']; assert self.length > 0.0
     264        if culvert_description_filename is not None:
     265            if not allclose(self.length, length, rtol=1.0e-2, atol=1.0e-2):
     266                msg = 'WARNING: barrel length specified in "%s" (%.2f m)'\
     267                    % (culvert_description_filename,
     268                       length)
     269                msg += ' does not match distance between specified'
     270                msg += ' end points (%.2f m)' %self.length
     271                print msg
     272       
     273        self.verbose = verbose
     274
     275
     276       
     277        # For use with update_interval                       
     278        self.last_update = 0.0
     279        self.update_interval = update_interval
     280       
     281
     282        # Print some diagnostics to log if requested
     283        if hasattr(self, 'log_filename'):
     284            s = 'Culvert Effective Length = %.2f m' %(self.length)
     285            log_to_file(self.log_filename, s)
     286   
     287            s = 'Culvert Direction is %s\n' %str(self.vector)
     288            log_to_file(self.log_filename, s)       
     289       
     290       
     291       
     292       
     293       
     294    def __call__(self, domain):
     295
     296        # Time stuff
     297        time = domain.get_time()
     298       
     299       
     300        update = False
     301        if self.update_interval is None:
     302            # Use next timestep as has been computed in domain.py       
     303            delta_t = domain.timestep           
     304            update = True
     305        else:   
     306            # Use update interval
     307            delta_t = self.update_interval           
     308            if time - self.last_update > self.update_interval or time == 0.0:
     309                update = True
     310
     311        if hasattr(self, 'log_filename'):           
     312            s = '\nTime = %.2f, delta_t = %f' %(time, delta_t)
     313            log_to_file(self.log_filename, s)
     314               
     315                               
     316        if update is True:
     317            self.compute_rates(delta_t)
     318           
     319   
     320        # Execute flow term for each opening
     321        # This is where Inflow objects are evaluated using the last rate
     322        # that has been calculated
     323        #
     324        # This will take place at every internal timestep and update the domain
     325        for opening in self.openings:
     326            opening(domain)
     327           
     328
     329
     330    def get_enquiry_indices(self):
     331        """Get indices for nearest centroids to self.enquiry_points
     332        """
     333       
     334        domain = self.domain
     335       
     336        enquiry_indices = []                 
     337        for point in self.enquiry_points:
     338            # Find nearest centroid
     339            N = len(domain)   
     340            points = domain.get_centroid_coordinates(absolute=True)
     341
     342            # Calculate indices in exchange area for this forcing term
     343           
     344            triangle_id = min_dist = sys.maxint
     345            for k in range(N):
     346                x, y = points[k,:] # Centroid
     347
     348                c = point                               
     349                distance = (x-c[0])**2+(y-c[1])**2
     350                if distance < min_dist:
     351                    min_dist = distance
     352                    triangle_id = k
     353
     354                   
     355            if triangle_id < sys.maxint:
     356                msg = 'found triangle with centroid (%f, %f)'\
     357                    %tuple(points[triangle_id, :])
     358                msg += ' for point (%f, %f)' %tuple(point)
     359               
     360                enquiry_indices.append(triangle_id)
     361            else:
     362                msg = 'Triangle not found for point (%f, %f)' %point
     363                raise Exception, msg
     364       
     365        return enquiry_indices
     366
     367       
     368    def check_culvert_inside_domain(self):
     369        """Check that all polygons and enquiry points lie within the mesh.
     370        """
     371        bounding_polygon = self.domain.get_boundary_polygon()
     372        P = self.culvert_polygons
     373        for key in P.keys():
     374            if key in ['exchange_polygon0',
     375                       'exchange_polygon1']:
     376                for point in list(P[key]) + self.enquiry_points:
     377                    msg = 'Point %s in polygon %s for culvert %s did not'\
     378                        %(str(point), key, self.label)
     379                    msg += 'fall within the domain boundary.'
     380                    assert is_inside_polygon(point, bounding_polygon), msg
     381           
     382
     383           
     384           
     385    def compute_rates(self, delta_t):
     386        """Compute new rates for inlet and outlet
     387        """
     388
     389        # Short hands
     390        domain = self.domain       
     391        dq = domain.quantities               
     392       
     393        # Time stuff
     394        time = domain.get_time()
     395        self.last_update = time
     396
     397           
     398        # Compute stage and energy at the
     399        # enquiry points at each end of the culvert
     400        openings = self.openings
     401        for i, opening in enumerate(openings):
     402            idx = self.enquiry_indices[i]               
     403           
     404            stage = dq['stage'].get_values(location='centroids',
     405                                               indices=[idx])[0]
     406            xmomentum = dq['xmomentum'].get_values(location='centroids',
     407                                                   indices=[idx])[0]
     408            ymomentum = dq['xmomentum'].get_values(location='centroids',
     409                                                   indices=[idx])[0]
     410           
     411            depth = h = stage-opening.elevation
     412            if h > minimum_allowed_height:
     413                u = xmomentum/(h + velocity_protection/h)
     414                v = ymomentum/(h + velocity_protection/h)
     415            else:
     416                u = v = 0.0
     417               
     418            energy_head = 0.5*(u*u + v*v)/g   
     419            opening.total_energy = energy_head + stage
     420            opening.specific_energy = energy_head + depth
     421            opening.stage = stage
     422            opening.depth = depth
     423           
     424
     425        # We now need to deal with each opening individually
     426        # Determine flow direction based on total energy difference
     427        delta_total_energy = openings[0].total_energy - openings[1].total_energy
     428        if delta_total_energy > 0:
     429            #print 'Flow U/S ---> D/S'
     430            inlet = openings[0]
     431            outlet = openings[1]
     432        else:
     433            #print 'Flow D/S ---> U/S'
     434            inlet = openings[1]
     435            outlet = openings[0]
     436           
     437            delta_total_energy = -delta_total_energy
     438
     439        self.inlet = inlet
     440        self.outlet = outlet
     441           
     442        msg = 'Total energy difference is negative'
     443        assert delta_total_energy > 0.0, msg
     444
     445        # Recompute slope and issue warning if flow is uphill
     446        # These values do not enter the computation
     447        delta_z = inlet.elevation - outlet.elevation
     448        culvert_slope = (delta_z/self.length)
     449        if culvert_slope < 0.0:
     450            # Adverse gradient - flow is running uphill
     451            # Flow will be purely controlled by uphill outlet face
     452            if self.verbose is True:
     453                print 'WARNING: Flow is running uphill. Watch Out!'
     454           
     455        if hasattr(self, 'log_filename'):
     456            s = 'Time=%.2f, inlet stage = %.2f, outlet stage = %.2f'\
     457                %(time, self.inlet.stage, self.outlet.stage)
     458            log_to_file(self.log_filename, s)
     459            s = 'Delta total energy = %.3f' %(delta_total_energy)
     460            log_to_file(log_filename, s)
     461
     462           
     463        # Determine controlling energy (driving head) for culvert
     464        if inlet.specific_energy > delta_total_energy:
     465            # Outlet control
     466            driving_head = delta_total_energy
     467        else:
     468            # Inlet control
     469            driving_head = inlet.specific_energy
     470           
     471
     472
     473        if self.inlet.depth <= self.trigger_depth:
     474            Q = 0.0
     475        else:
     476            # Calculate discharge for one barrel and
     477            # set inlet.rate and outlet.rate
     478           
     479            Q = self.compute_flow(driving_head)           
     480           
     481       
     482        # Adjust discharge for multiple barrels
     483        Q *= self.number_of_barrels
     484       
     485
     486        # Adjust Q downwards depending on available water at inlet
     487        stage = self.inlet.get_quantity_values(quantity_name='stage')
     488        elevation = self.inlet.get_quantity_values(quantity_name='elevation')
     489        depth = stage-elevation
     490       
     491       
     492        V = 0
     493        for i, d in enumerate(depth):
     494            V += d * domain.areas[i]
     495       
     496        #Vsimple = mean(depth)*self.inlet.exchange_area # Current volume in exchange area 
     497        #print 'Q', Q, 'dt', delta_t, 'Q*dt', Q*delta_t, 'V', V, 'Vsimple', Vsimple
     498
     499        dt = delta_t           
     500        if Q*dt > V:
     501       
     502            Q_reduced = 0.9*V/dt # Reduce with safety factor
     503           
     504            msg = '%.2fs: Computed extraction for this time interval (Q*dt) is ' % time
     505            msg += 'greater than current volume (V) at inlet.\n'
     506            msg += ' Q will be reduced from %.2f m^3/s to %.2f m^3/s.' % (Q, Q_reduced)
     507           
     508            #print msg
     509           
     510            if self.verbose is True:
     511                print msg
     512            if hasattr(self, 'log_filename'):                   
     513                log_to_file(self.log_filename, msg)                                       
     514           
     515            Q = Q_reduced
     516       
     517        self.inlet.rate = -Q
     518        self.outlet.rate = Q
     519
     520        # Log timeseries to file
     521        try:
     522            fid = open(self.timeseries_filename, 'a')       
     523        except:
     524            pass
     525        else:   
     526            fid.write('%.2f, %.2f\n' %(time, Q))
     527            fid.close()
     528
     529
     530    def compute_flow(self, driving_head):
     531       
     532        time = self.domain.get_time()   
     533        if self.culvert_description_filename is not None:
     534            try:
     535                Q = interpolate_linearly(driving_head,
     536                                         self.rating_curve[:,0],
     537                                         self.rating_curve[:,1])
     538            except Below_interval, e:
     539                Q = self.rating_curve[0,1]             
     540                msg = '%.2fs: ' % time
     541                msg += 'Delta head smaller than rating curve minimum: '
     542                msg += str(e)
     543                msg += '\n        '
     544                msg += 'I will use minimum discharge %.2f m^3/s ' % Q
     545                msg += 'for culvert "%s"' % self.label
     546               
     547                if hasattr(self, 'log_filename'):                   
     548                    log_to_file(self.log_filename, msg)
     549            except Above_interval, e:
     550                Q = self.rating_curve[-1,1]         
     551                msg = '%.2fs: ' % time                 
     552                msg += 'Delta head greater than rating curve maximum: '
     553                msg += str(e)
     554                msg += '\n        '
     555                msg += 'I will use maximum discharge %.2f m^3/s ' % Q
     556                msg += 'for culvert "%s"' % self.label
     557               
     558                if hasattr(self, 'log_filename'):                   
     559                    log_to_file(self.log_filename, msg)
     560        else:
     561            # User culvert routine
     562            Q, barrel_velocity, culvert_outlet_depth =\
     563                self.culvert_routine(self, driving_head, g)
     564           
     565        return Q
     566                           
     567   
    116568class Culvert_flow_rating:
    117569    """Culvert flow - transfer water from one hole to another
     
    131583                 end_point1=None,
    132584                 enquiry_point0=None,
    133                  enquiry_point1=None,                                           
     585                 enquiry_point1=None,
    134586                 update_interval=None,
    135587                 log_file=False,
     
    504956                 number_of_barrels=1,
    505957                 enquiry_point0=None,
    506                  enquiry_point1=None,                                           
     958                 enquiry_point1=None,
    507959                 update_interval=None,
    508960                 verbose=False):
     
    9401392
    9411393
    942 Culvert_flow = Culvert_flow_rating       
     1394Culvert_flow = Culvert_flow_general       
  • anuga_core/source/anuga/culvert_flows/culvert_routines.py

    r6121 r6123  
    4040
    4141
    42 def boyd_generalised_culvert_model(culvert, inlet, outlet, delta_Et, g):
     42def boyd_generalised_culvert_model(culvert, delta_Et, g):
    4343
    4444    """Boyd's generalisation of the US department of transportation culvert model
     
    5353    from anuga.utilities.numerical_tools import safe_acos as acos
    5454
    55        
     55    inlet = culvert.inlet
     56    outlet = culvert.outlet       
    5657    Q_outlet_tailwater = 0.0
    5758    inlet.rate = 0.0
  • anuga_core/source/anuga/culvert_flows/test_culvert_class.py

    r6121 r6123  
    1616    Transmissive_boundary, Time_boundary
    1717
    18 from anuga.culvert_flows.culvert_class import Culvert_flow_rating, Culvert_flow_energy
     18from anuga.culvert_flows.culvert_class import Culvert_flow, Culvert_flow_rating, Culvert_flow_energy
    1919from anuga.culvert_flows.culvert_routines import boyd_generalised_culvert_model
    2020     
     
    9696
    9797        filename=os.path.join(path, 'example_rating_curve.csv')
    98         culvert = Culvert_flow_rating(domain,
     98        culvert = Culvert_flow(domain,
    9999                               culvert_description_filename=filename,       
    100100                               end_point0=[9.0, 2.5],
     
    210210
    211211        filename = os.path.join(path, 'example_rating_curve.csv')
    212         culvert = Culvert_flow_rating(domain,
     212        culvert = Culvert_flow(domain,
    213213                               culvert_description_filename=filename,
    214214                               end_point0=[9.0, 2.5],
     
    314314
    315315        filename = os.path.join(path, 'example_rating_curve.csv')
    316         culvert = Culvert_flow_rating(domain,
    317                                       culvert_description_filename=filename,
    318                                       end_point0=[9.0, 2.5],
    319                                       end_point1=[13.0, 2.5],
    320                                       verbose=False)
     316        culvert = Culvert_flow(domain,
     317                               culvert_description_filename=filename,
     318                               end_point0=[9.0, 2.5],
     319                               end_point1=[13.0, 2.5],
     320                               verbose=False)
    321321
    322322        domain.forcing_terms.append(culvert)
     
    417417
    418418
    419         culvert = Culvert_flow_energy(domain,
    420                                       label='Culvert No. 1',
    421                                       description='This culvert is a test unit 1.2m Wide by 0.75m High',   
    422                                       end_point0=[9.0, 2.5],
    423                                       end_point1=[13.0, 2.5],
    424                                       width=1.20,height=0.75,
    425                                       culvert_routine=boyd_generalised_culvert_model,       
    426                                       number_of_barrels=1,
    427                                       update_interval=2,
    428                                       verbose=True)
     419        culvert = Culvert_flow(domain,
     420                               label='Culvert No. 1',
     421                               description='This culvert is a test unit 1.2m Wide by 0.75m High',   
     422                               end_point0=[9.0, 2.5],
     423                               end_point1=[13.0, 2.5],
     424                               width=1.20,height=0.75,
     425                               culvert_routine=boyd_generalised_culvert_model,
     426                               number_of_barrels=1,
     427                               update_interval=2,
     428                               verbose=True)
    429429       
    430430        domain.forcing_terms.append(culvert)
     
    516516
    517517
    518         culvert = Culvert_flow_energy(domain,
    519                                       label='Test culvert',
    520                                       description='4 m test culvert',   
    521                                       end_point0=[4.0, 2.5],
    522                                       end_point1=[8.0, 2.5],
    523                                       width=1.20,
    524                                       height=0.75,
    525                                       culvert_routine=boyd_generalised_culvert_model,       
    526                                       number_of_barrels=1,
    527                                       verbose=True)
     518        culvert = Culvert_flow(domain,
     519                               label='Test culvert',
     520                               description='4 m test culvert',   
     521                               end_point0=[4.0, 2.5],
     522                               end_point1=[8.0, 2.5],
     523                               width=1.20,
     524                               height=0.75,
     525                               culvert_routine=boyd_generalised_culvert_model,       
     526                               number_of_barrels=1,
     527                               verbose=True)
    528528                               
    529 
     529       
    530530        domain.forcing_terms.append(culvert)
    531531       
     
    540540#-------------------------------------------------------------
    541541if __name__ == "__main__":
    542     suite = unittest.makeSuite(Test_Culvert, 'test_that_culvert_rating_limits_flow_in_shallow_inlet_condition')
    543     #suite = unittest.makeSuite(Test_Culvert, 'test')
     542    #suite = unittest.makeSuite(Test_Culvert, 'test_that_culvert_rating_limits_flow_in_shallow_inlet_condition')
     543    suite = unittest.makeSuite(Test_Culvert, 'test')
    544544    runner = unittest.TextTestRunner()
    545545    runner.run(suite)
Note: See TracChangeset for help on using the changeset viewer.