Changeset 7984


Ignore:
Timestamp:
Sep 1, 2010, 6:11:25 PM (14 years ago)
Author:
steve
Message:

Added enquiry points to culverts

Location:
trunk/anuga_core/source/anuga/structures
Files:
8 edited

Legend:

Unmodified
Added
Removed

  • trunk/anuga_core/source/anuga/structures/boyd_box_culvert.py

    r7980 r7984  
    1515                 domain,
    1616                 end_points,
    17                  width=None,
    18                  height=None,
     17                 width,
     18                 height,
     19                 apron,
     20                 enquiry_gap,
    1921                 verbose=False):
    2022
    21         culvert.Culvert.__init__(self, domain, end_points, width, height, verbose)
     23        culvert.Culvert.__init__(self, domain, end_points, width, height, apron, enquiry_gap, verbose)
    2224
    2325        self.routine = boyd_box_routine.Boyd_box_routine(self)

  • trunk/anuga_core/source/anuga/structures/boyd_box_culvert_routine.py

    r7982 r7984  
    7070        local_debug ='false'
    7171       
    72         if self.inflow.get_average_height() > 0.01: #this value was 0.01:
     72        if self.inflow.get_enquiry_height() > 0.01: #this value was 0.01:
    7373            if local_debug =='true':
    7474                log.critical('Specific E & Deltat Tot E = %s, %s'
    75                              % (str(self.inflow.get_average_specific_energy()),
     75                             % (str(self.inflow.get_enquiry_specific_energy()),
    7676                                str(self.delta_total_energy)))
    7777                log.critical('culvert type = %s' % str(culvert_type))
    … …  
    7979
    8080            if self.log_filename is not None:
    81                 s = 'Specific energy  = %f m' % self.inflow.get_average_specific_energy()
     81                s = 'Specific energy  = %f m' % self.inflow.get_enquiry_specific_energy()
    8282                log_to_file(self.log_filename, s)
    8383
    8484            msg = 'Specific energy at inlet is negative'
    85             assert self.inflow.get_average_specific_energy() >= 0.0, msg
     85            assert self.inflow.get_enquiry_specific_energy() >= 0.0, msg
    8686
    8787            height = self.culvert_height
    … …  
    8989            flow_width = self.culvert_width
    9090
    91             Q_inlet_unsubmerged = 0.540*g**0.5*width*self.inflow.get_average_specific_energy()**1.50 # Flow based on Inlet Ctrl Inlet Unsubmerged
    92             Q_inlet_submerged = 0.702*g**0.5*width*height**0.89*self.inflow.get_average_specific_energy()**0.61  # Flow based on Inlet Ctrl Inlet Submerged
     91            Q_inlet_unsubmerged = 0.540*g**0.5*width*self.inflow.get_enquiry_specific_energy()**1.50 # Flow based on Inlet Ctrl Inlet Unsubmerged
     92            Q_inlet_submerged = 0.702*g**0.5*width*height**0.89*self.inflow.get_enquiry_specific_energy()**0.61  # Flow based on Inlet Ctrl Inlet Submerged
    9393
    9494            # FIXME(Ole): Are these functions really for inlet control?
    … …  
    120120                case = 'INLET CTRL Culvert is open channel flow we will for now assume critical depth'
    121121
    122             if self.delta_total_energy < self.inflow.get_average_specific_energy():
     122            if self.delta_total_energy < self.inflow.get_enquiry_specific_energy():
    123123                # Calculate flows for outlet control
    124124
    125125                # Determine the depth at the outlet relative to the depth of flow in the Culvert
    126                 if self.outflow.get_average_height() > height:        # The Outlet is Submerged
     126                if self.outflow.get_enquiry_height() > height:        # The Outlet is Submerged
    127127                    outlet_culvert_depth=height
    128128                    flow_area=width*height       # Cross sectional area of flow in the culvert
    … …  
    172172        # END CODE BLOCK for DEPTH  > Required depth for CULVERT Flow
    173173
    174         else: # self.inflow.get_average_height() < 0.01:
     174        else: # self.inflow.get_enquiry_height() < 0.01:
    175175            Q = barrel_velocity = outlet_culvert_depth = 0.0
    176176

  • trunk/anuga_core/source/anuga/structures/boyd_box_routine.py

    r7980 r7984  
    7070        local_debug ='false'
    7171       
    72         if self.inflow.get_average_height() > 0.01: #this value was 0.01:
     72        if self.inflow.get_enquiry_height() > 0.01: #this value was 0.01:
    7373            if local_debug =='true':
    7474                log.critical('Specific E & Deltat Tot E = %s, %s'
    75                              % (str(self.inflow.get_average_specific_energy()),
     75                             % (str(self.inflow.get_enquiry_specific_energy()),
    7676                                str(self.delta_total_energy)))
    7777                log.critical('culvert type = %s' % str(culvert_type))
    … …  
    7979
    8080            if self.log_filename is not None:
    81                 s = 'Specific energy  = %f m' % self.inflow.get_average_specific_energy()
     81                s = 'Specific energy  = %f m' % self.inflow.get_enquiry_specific_energy()
    8282                log_to_file(self.log_filename, s)
    8383
    8484            msg = 'Specific energy at inlet is negative'
    85             assert self.inflow.get_average_specific_energy() >= 0.0, msg
     85            assert self.inflow.get_enquiry_specific_energy() >= 0.0, msg
    8686
    8787            height = self.culvert_height
    … …  
    8989            flow_width = self.culvert_width
    9090
    91             Q_inlet_unsubmerged = 0.540*g**0.5*width*self.inflow.get_average_specific_energy()**1.50 # Flow based on Inlet Ctrl Inlet Unsubmerged
    92             Q_inlet_submerged = 0.702*g**0.5*width*height**0.89*self.inflow.get_average_specific_energy()**0.61  # Flow based on Inlet Ctrl Inlet Submerged
     91            Q_inlet_unsubmerged = 0.540*g**0.5*width*self.inflow.get_enquiry_specific_energy()**1.50 # Flow based on Inlet Ctrl Inlet Unsubmerged
     92            Q_inlet_submerged = 0.702*g**0.5*width*height**0.89*self.inflow.get_enquiry_specific_energy()**0.61  # Flow based on Inlet Ctrl Inlet Submerged
    9393
    9494            # FIXME(Ole): Are these functions really for inlet control?
    … …  
    120120                case = 'INLET CTRL Culvert is open channel flow we will for now assume critical depth'
    121121
    122             if self.delta_total_energy < self.inflow.get_average_specific_energy():
     122            if self.delta_total_energy < self.inflow.get_enquiry_specific_energy():
    123123                # Calculate flows for outlet control
    124124
    125125                # Determine the depth at the outlet relative to the depth of flow in the Culvert
    126                 if self.outflow.get_average_height() > height:        # The Outlet is Submerged
     126                if self.outflow.get_enquiry_height() > height:        # The Outlet is Submerged
    127127                    outlet_culvert_depth=height
    128128                    flow_area=width*height       # Cross sectional area of flow in the culvert
    … …  
    172172        # END CODE BLOCK for DEPTH  > Required depth for CULVERT Flow
    173173
    174         else: # self.inflow.get_average_height() < 0.01:
     174        else: # self.inflow.get_enquiry_height() < 0.01:
    175175            Q = barrel_velocity = outlet_culvert_depth = 0.0
    176176

  • trunk/anuga_core/source/anuga/structures/culvert.py

    r7980 r7984  
    2020                 domain,
    2121                 end_points,
    22                  width=None,
    23                  height=None,
    24                  verbose=False):
     22                 width,
     23                 height,
     24                 apron,
     25                 enquiry_gap,
     26                 verbose):
    2527       
    2628        # Input check
    2729       
    2830        self.domain = domain
    29 
    3031        self.end_points = end_points
    31  
    32         self.width = width
     32        self.width  = width
    3333        self.height = height
    34        
     34        self.apron  = apron
     35        self.enquiry_gap = enquiry_gap
    3536        self.verbose=verbose
    3637       
    … …  
    4243        self.inlets = []
    4344        polygon0 = self.inlet_polygons[0]
     45        ep0 = self.inlet_equiry_pts[0]
    4446        outward_vector0 = self.culvert_vector
    45         self.inlets.append(inlet.Inlet(self.domain, polygon0, outward_vector0))
     47        self.inlets.append(inlet.Inlet(self.domain, polygon0, ep0, outward_vector0))
    4648
    4749        polygon1 = self.inlet_polygons[1]
     50        ep1 = self.inlet_equiry_pts[1]
    4851        outward_vector1  = - self.culvert_vector
    49         self.inlets.append(inlet.Inlet(self.domain, polygon1, outward_vector1))
     52        self.inlets.append(inlet.Inlet(self.domain, polygon1, ep1, outward_vector1))
    5053
    5154
    … …  
    7881        # Short hands
    7982        w = 0.5*self.width*self.culvert_normal # Perpendicular vector of 1/2 width
    80         h = 0.5*self.height*self.culvert_vector    # Vector of length=height in the
     83        h = self.apron*self.culvert_vector    # Vector of length=height in the
    8184                             # direction of the culvert
    8285
     86        gap = (1 + self.enquiry_gap)*h
     87
    8388        self.inlet_polygons = []
     89        self.inlet_equiry_pts = []
    8490
    85         # Build exchange polygon and enquiry points 0 and 1
     91        # Build exchange polygon and enquiry point
    8692        for i in [0, 1]:
    8793            i0 = (2*i-1)
    … …  
    9096            p2 = p1 + i0*h
    9197            p3 = p0 + i0*h
     98            ep = self.end_points[i] + i0*gap
     99
    92100            self.inlet_polygons.append(num.array([p0, p1, p2, p3]))
     101            self.inlet_equiry_pts.append(ep)
    93102
    94103        # Check that enquiry points are outside inlet polygons
    95104        for i in [0,1]:
    96105            polygon = self.inlet_polygons[i]
    97             # FIXME (SR) Probably should calculate the area of all the triangles
    98             # associated with this polygon, as there is likely to be some
    99             # inconsistency between triangles and ploygon
     106            ep = self.inlet_equiry_pts[i]
     107           
    100108            area = polygon_area(polygon)
    101109           
    … …  
    103111            msg += ' has area = %f' % area
    104112            assert area > 0.0, msg
    105            
     113
     114            msg = 'Enquiry point falls inside an exchange polygon.'
     115            assert not inside_polygon(ep, polygon), msg
    106116   
    107117    def get_inlets(self):
    … …  
    123133   
    124134        return self.height
     135
     136
     137    def get_culvert_apron(self):
     138
     139        return self.apron
    125140       

  • trunk/anuga_core/source/anuga/structures/culvert_operator.py

    r7981 r7984  
    2222                 width,
    2323                 height=None,
     24                 apron=None,
     25                 enquiry_gap=0.2,
    2426                 verbose=False):
    2527       
    2628        self.domain = domain
    2729        self.domain.set_fractional_step_operator(self)
    28         end_points = [end_point0, end_point1]
     30        self.end_points = [end_point0, end_point1]
    2931       
    3032        if height is None:
    3133            height = width
    3234
    33         self.width = width
     35        if apron is None:
     36            apron = width
     37
     38        self.width  = width
    3439        self.height = height
     40        self.apron  = apron
     41        self.enquiry_gap = enquiry_gap
     42        self.verbose = verbose
    3543       
    36         self.culvert = Boyd_box_culvert(self.domain, end_points, self.width, self.height)
    37         print self.culvert
     44        self.culvert = Boyd_box_culvert(self.domain,
     45                                        self.end_points,
     46                                        self.width,
     47                                        self.height,
     48                                        self.apron,
     49                                        self.enquiry_gap,
     50                                        self.verbose)
     51       
    3852        self.routine = self.culvert.routine
    39         print self.routine
    40         self.inlets = self.culvert.get_inlets()
    41    
    42         self.print_stats()
     53
     54        self.inlets  = self.culvert.get_inlets()
     55
     56        if self.verbose:
     57            self.print_stats()
    4358
    4459
    … …  
    4863       
    4964        Q, barrel_speed, outlet_depth = self.routine()
     65
    5066
    5167        inflow  = self.routine.get_inflow()
    … …  
    5470
    5571        old_inflow_height = inflow.get_average_height()
    56                 old_inflow_xmom = inflow.get_average_xmom()
    57                 old_inflow_ymom = inflow.get_average_ymom()
    58                
    59                 if old_inflow_height > 0.0 :
    60                         Qstar = Q/old_inflow_height/inflow.get_area()
    61                 else:
    62                         Qstar = 0.0
    63 
    64                 factor = 1.0/(1.0 + Qstar*timestep)
    65 
    66                
    67                
    68                 new_inflow_height = old_inflow_height*factor
    69                 new_inflow_xmom = old_inflow_xmom*factor
    70                 new_inflow_ymom = old_inflow_ymom*factor
     72        old_inflow_xmom = inflow.get_average_xmom()
     73        old_inflow_ymom = inflow.get_average_ymom()
     74               
     75        if old_inflow_height > 0.0 :
     76            Qstar = Q/old_inflow_height
     77        else:
     78            Qstar = 0.0
     79
     80        factor = 1.0/(1.0 + Qstar*timestep/inflow.get_area())
     81
     82               
     83               
     84        new_inflow_height = old_inflow_height*factor
     85        new_inflow_xmom = old_inflow_xmom*factor
     86        new_inflow_ymom = old_inflow_ymom*factor
    7187               
    7288
    7389        inflow.set_heights(new_inflow_height)
     90
     91        #inflow.set_xmoms(Q/inflow.get_area())
     92        #inflow.set_ymoms(0.0)
     93
     94
    7495        inflow.set_xmoms(new_inflow_xmom)
    7596        inflow.set_ymoms(new_inflow_ymom)
    7697
    77                
    78                 # set outflow
    79                 if old_inflow_height > 0.0 :
    80                         timestep_star = timestep*new_inflow_height/old_inflow_height
    81                 else:
    82                         timestep_star = 0.0
    83                
    84                 print Q, barrel_speed, outlet_depth, Qstar, factor, timestep_star
    85                
     98
     99        loss = (old_inflow_height - new_inflow_height)*inflow.get_area()
     100
     101               
     102        # set outflow
     103        if old_inflow_height > 0.0 :
     104            timestep_star = timestep*new_inflow_height/old_inflow_height
     105        else:
     106            timestep_star = 0.0
     107
    86108               
    87109        outflow_extra_height = Q*timestep_star/outflow.get_area()
    … …  
    90112               
    91113
    92         outflow.set_heights(outflow.get_average_height() + outflow_extra_height)
    93         outflow.set_xmoms(outflow.get_average_xmom() + outflow_extra_momentum[0] )
    94         outflow.set_ymoms(outflow.get_average_ymom() + outflow_extra_momentum[1] )
     114        gain = outflow_extra_height*outflow.get_area()
     115       
     116        #print Q, Q*timestep, barrel_speed, outlet_depth, Qstar, factor, timestep_star
     117        #print '  ', loss, gain
     118
     119
     120        new_outflow_height = outflow.get_average_height() + outflow_extra_height
     121        new_outflow_xmom = outflow.get_average_xmom() + outflow_extra_momentum[0]
     122        new_outflow_ymom = outflow.get_average_ymom() + outflow_extra_momentum[1]
     123
     124        outflow.set_heights(new_outflow_height)
     125
     126        outflow.set_xmoms(barrel_speed*new_outflow_height*outflow_direction[0])
     127        outflow.set_ymoms(barrel_speed*new_outflow_height*outflow_direction[1])
     128
     129        #outflow.set_xmoms(new_outflow_xmom)
     130        #outflow.set_ymoms(new_outflow_ymom)
     131       
     132        #print '   outflow volume ',outflow.get_total_water_volume()
    95133
    96134    def print_stats(self):
    … …  
    99137        print 'Generic Culvert Operator'
    100138        print '====================================='
     139
     140        print 'Culvert'
     141        print self.culvert
     142
     143        print 'Culvert Routine'
     144        print self.routine
    101145       
    102146        for i, inlet in enumerate(self.inlets):

  • trunk/anuga_core/source/anuga/structures/culvert_routine.py

    r7980 r7984  
    5050        # Determine flow direction based on total energy difference
    5151
    52         self.delta_total_energy = self.inlets[0].get_average_total_energy() - self.inlets[1].get_average_total_energy()
     52        self.delta_total_energy = self.inlets[0].get_enquiry_total_energy() - self.inlets[1].get_enquiry_total_energy()
    5353
    5454        self.inflow  = self.inlets[0]

  • trunk/anuga_core/source/anuga/structures/inlet.py

    r7980 r7984  
    99    """
    1010
    11     def __init__(self, domain, polygon, outward_culvert_vector=None):
     11    def __init__(self, domain, polygon, enquiry_pt,  outward_culvert_vector=None, verbose=False):
    1212
    1313        self.domain = domain
    1414        self.domain_bounding_polygon = self.domain.get_boundary_polygon()
    1515        self.polygon = polygon
     16        self.enquiry_pt = enquiry_pt
    1617        self.outward_culvert_vector = outward_culvert_vector
     18        self.verbose = verbose
    1719
    1820        # FIXME (SR) Using get_triangle_containing_point which needs to be sped up
    1921
    20         self.compute_triangle_indices()
     22        self.compute_indices()
    2123        self.compute_area()
    2224
    2325
    24     def compute_triangle_indices(self):
     26    def compute_indices(self):
    2527
    2628        # Get boundary (in absolute coordinates)
    … …  
    2830        domain_centroids = self.domain.get_centroid_coordinates(absolute=True)
    2931
    30         self.inlet_polygon = self.polygon
    3132
    3233        # Check that polygon lies within the mesh.
    33         for point in self.inlet_polygon:
    34                 msg = 'Point %s in polygon for forcing term' %  str(point)
     34        for point in self.polygon + self.enquiry_pt:
     35                msg = 'Point %s ' %  str(point)
    3536                msg += ' did not fall within the domain boundary.'
    3637                assert is_inside_polygon(point, bounding_polygon), msg
    3738
    3839
    39         self.triangle_indices = inside_polygon(domain_centroids, self.inlet_polygon, verbose=True)
     40        self.triangle_indices = inside_polygon(domain_centroids, self.polygon, verbose=self.verbose)
    4041
    4142        if len(self.triangle_indices) == 0:
    42             region = 'Inlet polygon=%s' % (self.inlet_polygon)
    43             msg = 'No triangles have been identified in '
    44             msg += 'specified region: %s' % region
     43            region = 'Inlet polygon=%s' % (self.polygon)
     44            msg = 'No triangles have been identified in region '
    4545            raise Exception, msg
    4646
     47        self.enquiry_index = self.domain.get_triangle_containing_point(self.enquiry_pt)
    4748
    4849
    … …  
    5455        if len(self.triangle_indices) == 0:
    5556            region = 'Inlet polygon=%s' % (self.inlet_polygon)
    56             msg = 'No triangles have been identified in '
    57             msg += 'specified region: %s' % region
     57            msg = 'No triangles have been identified in region '
    5858            raise Exception, msg
    5959       
    … …  
    115115        return num.sum(self.get_ymoms()*self.get_areas())/self.area
    116116   
     117
    117118    def get_heights(self):
    118119   
    … …  
    175176
    176177
     178    def get_enquiry_stage(self):
     179
     180        return self.domain.quantities['stage'].centroid_values[self.enquiry_index]
     181
     182
     183    def get_enquiry_xmom(self):
     184
     185        return self.domain.quantities['xmomentum'].centroid_values[self.enquiry_index]
     186
     187    def get_enquiry_ymom(self):
     188
     189        return self.domain.quantities['ymomentum'].centroid_values[self.enquiry_index]
     190
     191
     192    def get_enquiry_elevation(self):
     193
     194        return self.domain.quantities['elevation'].centroid_values[self.enquiry_index]
     195
     196    def get_enquiry_height(self):
     197
     198        return self.get_enquiry_stage() - self.get_enquiry_elevation()
     199
     200
     201    def get_enquiry_velocity(self):
     202
     203            height = self.get_enquiry_height()
     204            u = self.get_enquiry_xmom()/(height + velocity_protection/height)
     205            v = self.get_enquiry_ymom()/(height + velocity_protection/height)
     206
     207            return u, v
     208
     209
     210    def get_enquiry_xvelocity(self):
     211
     212            height = self.get_enquiry_height()
     213            return self.get_enquiry_xmom()/(height + velocity_protection/height)
     214
     215    def get_enquiry_yvelocity(self):
     216
     217            height = self.get_enquiry_height()
     218            return self.get_enquiry_ymom()/(height + velocity_protection/height)
     219
     220
     221    def get_enquiry_speed(self):
     222
     223            u, v = self.get_enquiry_velocity()
     224
     225            return math.sqrt(u**2 + v**2)
     226
     227
     228    def get_enquiry_velocity_head(self):
     229
     230        return 0.5*self.get_enquiry_speed()**2/g
     231
     232
     233    def get_enquiry_total_energy(self):
     234
     235        return self.get_enquiry_velocity_head() + self.get_enquiry_stage()
     236
     237
     238    def get_enquiry_specific_energy(self):
     239
     240        return self.get_enquiry_velocity_head() + self.get_enquiry_height()
     241
     242
    177243    def set_heights(self,height):
    178244
    … …  
    187253    def set_xmoms(self,xmom):
    188254
    189         self.xmoms=self.domain.quantities['xmomentum'].centroid_values.put(self.triangle_indices, xmom)
     255        self.domain.quantities['xmomentum'].centroid_values.put(self.triangle_indices, xmom)
    190256
    191257
    192258    def set_ymoms(self,ymom):
    193259
    194         self.xmoms=self.domain.quantities['ymomentum'].centroid_values.put(self.triangle_indices, ymom)
     260        self.domain.quantities['ymomentum'].centroid_values.put(self.triangle_indices, ymom)
    195261
    196262

  • trunk/anuga_core/source/anuga/structures/testing_wide_bridge.py

    r7980 r7984  
    4646
    4747points, vertices, boundary = rectangular_cross(int(length/dx), int(width/dy),
    48                                                len1=length, len2=width)
     48                                                    len1=length, len2=width)
    4949domain = Domain(points, vertices, boundary)   
    5050domain.set_name('Test_WIDE_BRIDGE')                 # Output name
    … …  
    168168"""
    169169from anuga.structures.culvert_operator import Culvert_operator
     170#culvert0 = Culvert_operator(domain,
     171#                            end_point0=[40.0, 75.0],
     172#                            end_point1=[50.0, 75.0],
     173#                            width=50.0,
     174#                            height=10.0,
     175#                            apron=5.0,
     176#                            verbose=False)
     177
     178
    170179culvert1 = Culvert_operator(domain,
    171                             end_point0=[40.0, 75.0],
    172                             end_point1=[50.0, 75.0],
    173                             width=50.0,
     180                            end_point0=[40.0, 87.5],
     181                            end_point1=[50.0, 87.5],
     182                            width=25.0,
    174183                            height=10.0,
     184                            apron=5.0,
     185                            verbose=False)
     186
     187
     188culvert2 = Culvert_operator(domain,
     189                            end_point0=[40.0, 62.5],
     190                            end_point1=[50.0, 62.5],
     191                            width=25.0,
     192                            height=10.0,
     193                            apron=5.0,
    175194                            verbose=False)
    176195
    … …  
    220239for t in domain.evolve(yieldstep = 1, finaltime = 100):
    221240    print domain.timestepping_statistics()
     241    print domain.volumetric_balance_statistics()
    222242   
    223243
Note: See TracChangeset for help on using the changeset viewer.