Changeset 7971


Ignore:
Timestamp:
Aug 24, 2010, 11:06:01 PM (8 years ago)
Author:
habili
Message:

Refactoring. Moving Inlet class to Inlet.py

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/anuga_core/source/anuga/structures/culvert_operator.py

    r7968 r7971  
    66from anuga.geometry.polygon import plot_polygons, polygon_area
    77
    8 
    98from anuga.utilities.numerical_tools import mean
    109from anuga.utilities.numerical_tools import ensure_numeric, sign
     
    1413import anuga.utilities.log as log
    1514
     15import Inlet
     16
    1617import numpy as num
    17 from math import sqrt
     18import math
    1819
    1920class Below_interval(Exception): pass
    2021class Above_interval(Exception): pass
    21 
    22 
    23 
    24 class Inlet:
    25     """Contains information associated with each inlet
    26     """
    27 
    28     def __init__(self,domain,polygon,enquiry_point,inlet_vector):
    29 
    30         self.domain = domain
    31         self.domain_bounding_polygon = self.domain.get_boundary_polygon()
    32         self.polygon = polygon
    33         self.enquiry_point = enquiry_point
    34         self.inlet_vector = inlet_vector
    35 
    36         # FIXME (SR) Using get_triangle_containing_point which needs to be sped up
    37         self.enquiry_index = self.domain.get_triangle_containing_point(self.enquiry_point)
    38 
    39         self.compute_inlet_triangle_indices()
    40 
    41 
    42     def compute_inlet_averages(self):
    43 
    44 
    45 
    46         self.cell_indices = self.exchange_triangle_indices[0]
    47         self.areas    = areas[cell_indices]
    48 
    49 
    50         # Inlet Averages
    51         self.heights            = stage[cell_indices]-elevation[cell_indices]
    52         self.total_water_volume = num.sum(self.heights*self.areas)
    53         self.average_height     = self.total_water_volume/self.total_inlet_area
    54 
    55 
    56 
    57 
    58     def compute_inlet_triangle_indices(self):
    59 
    60         # Get boundary (in absolute coordinates)
    61         bounding_polygon = self.domain_bounding_polygon
    62         centroids = self.domain.get_centroid_coordinates(absolute=True)
    63 
    64         inlet_polygon = self.polygon
    65 
    66         # Check that polygon lies within the mesh.
    67         for point in inlet_polygon:
    68                 msg = 'Point %s in polygon for forcing term' % str(point)
    69                 msg += ' did not fall within the domain boundary.'
    70                 assert is_inside_polygon(point, bounding_polygon), msg
    71 
    72 
    73         self.triangle_indices = inside_polygon(centroids, inlet_polygon)
    74 
    75         if len(self.triangle_indices) == 0:
    76             region = 'Inlet polygon=%s' % (inlet_polygon)
    77             msg = 'No triangles have been identified in '
    78             msg += 'specified region: %s' % region
    79             raise Exception, msg
    80 
    81         # Compute exchange area as the sum of areas of triangles identified
    82         # by polygon
    83         self.area = 0.0
    84         for j in self.triangle_indices:
    85             self.area += self.domain.areas[j]
    86 
    87 
    88         msg = 'Inlet exchange area has area = %f' % self.area
    89         assert self.area > 0.0
    90 
    9122
    9223
     
    11748        self.domain.set_fractional_step_operator(self)
    11849
    119         self.end_points= [end_point0, end_point1]
     50        self.end_points = [end_point0, end_point1]
    12051        self.enquiry_gap_factor = enquiry_gap_factor
    12152       
     
    13465        #FIXME (SR) Put this into a foe loop to deal with more inlets
    13566        self.inlets = []
    136         polygon0 = self.exchange_polygons[0]
     67        polygon0 = self.inlet_polygons[0]
    13768        enquiry_pt0 = self.enquiry_points[0]
    13869        inlet0_vector = self.culvert_vector
    139 
    140         self.inlets.append(Inlet(self.domain,polygon0,enquiry_pt0,inlet0_vector))
    141 
    142         polygon1 = self.exchange_polygons[1]
     70        self.inlets.append(Inlet.Inlet(self.domain, polygon0, enquiry_pt0, inlet0_vector))
     71
     72        polygon1 = self.inlet_polygons[1]
    14373        enquiry_pt1 = self.enquiry_points[1]
    14474        inlet1_vector = - self.culvert_vector
    145 
    146         self.inlets.append(Inlet(self.domain,polygon1,enquiry_pt1, inlet1_vector))
    147 
    148 
    149         # aliases to quantity centroid values and cell areas
    150         self.areas      = self.domain.areas
    151         self.stage      = self.domain.quantities['stage'].centroid_values
    152         self.elevation  = self.domain.quantities['elevation'].centroid_values
    153         self.xmom       = self.domain.quantities['xmomentum'].centroid_values
    154         self.ymom       = self.domain.quantities['ymomentum'].centroid_values
    155 
     75        self.inlets.append(Inlet.Inlet(self.domain, polygon1, enquiry_pt1, inlet1_vector))
    15676 
     77        self.areas = self.domain.areas
     78        self.stages = self.domain.quantities['stage'].centroid_values
     79        self.elevations = self.domain.quantities['elevation'].centroid_values
     80        self.xmoms = self.domain.quantities['xmomentum'].centroid_values
     81        self.ymoms = self.domain.quantities['ymomentum'].centroid_values
     82   
    15783        self.print_stats()
    15884
    15985
    160 
    16186    def __call__(self):
    162 
    16387
    16488        # Time stuff
     
    16690        timestep = self.domain.get_timestep()
    16791
    168 
    16992        inlet0 = self.inlets[0]
    17093        inlet1 = self.inlets[1]
    171 
    17294
    17395        # Aliases to cell indices
     
    17597        inlet1_indices = inlet1.triangle_indices
    17698
    177 
    17899        # Inlet0 averages
    179         inlet0_heights  = self.stage[inlet0_indices]-self.elevation[inlet0_indices]
    180         inlet0_areas    = self.areas[inlet0_indices]
    181 
     100        inlet0_heights = inlet0.get_heights()
     101        inlet0_areas = inlet0.get_areas()
    182102        inlet0_water_volume = num.sum(inlet0_heights*inlet0_areas)
    183 
    184103        average_inlet0_height = inlet0_water_volume/inlet0.area
    185104
    186105        # Inlet1 averages
    187         inlet1_heights  = self.stage[inlet1_indices]-self.elevation[inlet1_indices]
    188         inlet1_areas    = self.areas[inlet1_indices]
    189 
     106        inlet1_heights = inlet1.get_heights()
     107        inlet1_areas = inlet1.get_areas()
    190108        inlet1_water_volume = num.sum(inlet1_heights*inlet1_areas)
    191 
    192109        average_inlet1_height = inlet1_water_volume/inlet1.area
    193 
    194110
    195111        # Transfer
    196112        transfer_water = timestep*inlet0_water_volume
    197 
    198         self.stage[inlet0_indices] = self.elevation[inlet0_indices] + average_inlet0_height - transfer_water
    199         self.xmom[inlet0_indices]  = 0.0
    200         self.ymom[inlet0_indices]  = 0.0
    201 
    202 
    203         self.stage[inlet1_indices] = self.elevation[inlet1_indices] + average_inlet1_height + transfer_water
    204         self.xmom[inlet1_indices]  = 0.0
    205         self.ymom[inlet1_indices]  = 0.0
     113       
     114        self.stages.put(inlet0_indices, inlet0.get_elevations() + average_inlet0_height - transfer_water)
     115        self.xmoms.put(inlet0_indices, 0.0)
     116        self.ymoms.put(inlet0_indices, 0.0)
     117
     118        self.stages.put(inlet1_indices, inlet1.get_elevations() + average_inlet1_height + transfer_water)
     119        self.xmoms.put(inlet1_indices, 0.0)
     120        self.ymoms.put(inlet1_indices, 0.0)
    206121
    207122
     
    232147
    233148
    234  
    235 
    236 
    237 
    238 
    239149    def set_store_hydrograph_discharge(self, filename=None):
    240150
     
    251161        fid.close()
    252162
     163
    253164    def create_culvert_polygons(self):
     165
    254166        """Create polygons at the end of a culvert inlet and outlet.
    255167        At either end two polygons will be created; one for the actual flow to pass through and one a little further away
     
    265177
    266178        self.culvert_vector = num.array([dx, dy])
    267         self.culvert_length = sqrt(num.sum(self.culvert_vector**2))
     179        self.culvert_length = math.sqrt(num.sum(self.culvert_vector**2))
    268180        assert self.culvert_length > 0.0, 'The length of culvert is less than 0'
    269181
     
    278190        gap = (1 + self.enquiry_gap_factor)*h
    279191
    280         self.exchange_polygons = []
     192        self.inlet_polygons = []
    281193        self.enquiry_points = []
    282194
     
    288200            p2 = p1 + i0*h
    289201            p3 = p0 + i0*h
    290             self.exchange_polygons.append(num.array([p0, p1, p2, p3]))
     202            self.inlet_polygons.append(num.array([p0, p1, p2, p3]))
    291203            self.enquiry_points.append(self.end_points[i] + i0*gap)
    292204
    293 
    294 
    295 
    296         # Check that enquiry points are outside exchange polygons
     205        # Check that enquiry points are outside inlet polygons
    297206        for i in [0,1]:
    298             polygon = self.exchange_polygons[i]
     207            polygon = self.inlet_polygons[i]
    299208            # FIXME (SR) Probably should calculate the area of all the triangles
    300209            # associated with this polygon, as there is likely to be some
     
    313222                assert not inside_polygon(point, polygon), msg
    314223
    315        
    316 
    317 
    318 
    319        
    320 
    321            
    322 
     224   
    323225    def adjust_flow_for_available_water_at_inlet(self, Q, delta_t):
    324226        """Adjust Q downwards depending on available water at inlet
     
    463365        """Compute new rates for inlet and outlet
    464366        """
    465 
    466         # Short hands
    467         domain = self.domain       
    468         dq = domain.quantities               
    469        
    470367        # Time stuff
    471         time = domain.get_time()
     368        time = self.domain.get_time()
    472369        self.last_update = time
    473370
    474            
    475371        if hasattr(self, 'log_filename'):
    476372            log_filename = self.log_filename
     
    478374        # Compute stage, energy and velocity at the
    479375        # enquiry points at each end of the culvert
    480         openings = self.openings
    481         for i, opening in enumerate(openings):
    482             idx = self.enquiry_indices[i]               
    483            
    484             stage = dq['stage'].get_values(location='centroids',
    485                                            indices=[idx])[0]
    486             depth = h = stage-opening.elevation
    487                                                            
    488            
    489             # Get velocity                                 
    490             xmomentum = dq['xmomentum'].get_values(location='centroids',
    491                                                    indices=[idx])[0]
    492             ymomentum = dq['xmomentum'].get_values(location='centroids',
    493                                                    indices=[idx])[0]
    494 
    495             if h > minimum_allowed_height:
    496                 u = xmomentum/(h + velocity_protection/h)
    497                 v = ymomentum/(h + velocity_protection/h)
     376        total_energies = []
     377        specific_energies = []
     378        velocities = []
     379       
     380        for i, inlet in enumerate(self.inlets):
     381           
     382            stage = inlet.get_average_stage()
     383            depth = stage - inlet.get_average_elevation()
     384                                                                       
     385            if depth > minimum_allowed_height:
     386                u, v = inlet.get_average_velocities
    498387            else:
    499                 u = v = 0.0
     388                u = 0.0
     389                v = 0.0
    500390               
    501             v_squared = u*u + v*v
    502            
     391            v_squared =  u**2 + v**2   
     392               
    503393            if self.use_velocity_head is True:
    504394                velocity_head = 0.5*v_squared/g   
     
    506396                velocity_head = 0.0
    507397           
    508             opening.total_energy = velocity_head + stage
    509             opening.specific_energy = velocity_head + depth
    510             opening.stage = stage
    511             opening.depth = depth
    512             opening.velocity = sqrt(v_squared)
    513            
     398            total_energies[i] = velocity_head + stage
     399            specific_energies[i] = velocity_head + depth
     400            velocities[i] = math.sqrt(v_squared)           
    514401
    515402        # We now need to deal with each opening individually
    516403        # Determine flow direction based on total energy difference
    517         delta_total_energy = openings[0].total_energy - openings[1].total_energy
    518         if delta_total_energy > 0:
    519             inlet = openings[0]
    520             outlet = openings[1]
    521 
    522             # FIXME: I think this whole momentum jet thing could be a bit more elegant
    523             inlet.momentum = self.opening_momentum[0]
    524             outlet.momentum = self.opening_momentum[1]
     404        delta_total_energy = total_energies[0] - total_energies[1]
     405        if delta_total_energy >= 0:
     406            inlet = inlets[0]
     407            outlet = inlets[1]
    525408        else:
    526             inlet = openings[1]
    527             outlet = openings[0]
    528            
    529             inlet.momentum = self.opening_momentum[1]
    530             outlet.momentum = self.opening_momentum[0]
    531 
    532             delta_total_energy = -delta_total_energy
    533 
    534         self.inlet = inlet
    535         self.outlet = outlet
    536            
    537         msg = 'Total energy difference is negative'
    538         assert delta_total_energy >= 0.0, msg
     409            msg = 'Total energy difference is negative'
     410            assert delta_total_energy >= 0.0, msg
    539411
    540412        # Recompute slope and issue warning if flow is uphill
    541413        # These values do not enter the computation
    542         delta_z = inlet.elevation - outlet.elevation
    543         culvert_slope = (delta_z/self.length)
     414        delta_z = inlet.get_average_elevation() - outlet.get_average_elevation()
     415        culvert_slope = (delta_z/self.culvert_length)
    544416        if culvert_slope < 0.0:
    545417            # Adverse gradient - flow is running uphill
     
    554426            s = 'Delta total energy = %.3f' %(delta_total_energy)
    555427            log_to_file(log_filename, s)
    556 
    557428           
    558429        # Determine controlling energy (driving head) for culvert
     
    564435            driving_head = inlet.specific_energy
    565436           
    566 
    567 
    568437        if self.inlet.depth <= self.trigger_depth:
    569438            Q = 0.0
     
    617486                                         log_filename=self.log_filename)
    618487           
    619            
    620        
    621488        # Adjust discharge for multiple barrels
    622489        Q *= self.number_of_barrels
     
    628495        self.outlet.rate = Q
    629496
    630 
    631497        # Momentum jet stuff
    632498        if self.use_momentum_jet is True:
    633 
    634499
    635500            # Compute barrel momentum
     
    646511                s = 'Directional momentum = (%f, %f)' %(outlet_mom_x, outlet_mom_y)
    647512                log_to_file(self.log_filename, s)
    648 
    649513
    650514            # Update momentum       
Note: See TracChangeset for help on using the changeset viewer.