Changeset 7962


Ignore:
Timestamp:
Aug 20, 2010, 5:27:03 PM (9 years ago)
Author:
steve
Message:

Got a flow from one region to another.

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

Legend:

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

    r7958 r7962  
    4646        self.domain = domain
    4747
     48        self.domain.set_fractional_step_operator(self)
     49
    4850        self.end_points= [end_point0, end_point1]
    4951        self.enquiry_gap_factor = enquiry_gap_factor
     
    6264        self.compute_enquiry_indices()
    6365        self.check_culvert_inside_domain()
    64 
    65         # Establish initial values at each enquiry point
    66         self.enquiry_quantity_values = []
    67         dq = domain.quantities
    68         for i in [0, 1]:
    69             idx = self.enquiry_indices[i]
    70             elevation = dq['elevation'].get_values(location='centroids',
    71                                                    indices=[idx])[0]
    72             stage = dq['stage'].get_values(location='centroids',
    73                                            indices=[idx])[0]
    74             depth = stage - elevation
    75            
    76             quantity_values = {'stage' : stage, 'elevation' : elevation, 'depth' : depth }
    77             self.enquiry_quantity_values.append(quantity_values)
     66        self.compute_exchange_triangle_indices()
     67
     68
     69        #self.print_stats()
     70
     71
     72
     73    def __call__(self):
     74
     75
     76        # Time stuff
     77        time     = self.domain.get_time()
     78        timestep = self.domain.get_timestep()
     79
     80
     81        inlet_indices  = self.exchange_triangle_indices[0]
     82        outlet_indices = self.exchange_triangle_indices[1]
     83
     84        areas      = self.domain.areas
     85        stage     = self.domain.quantities['stage'].centroid_values
     86        elevation = self.domain.quantities['elevation'].centroid_values
     87       
     88        xmom = self.domain.quantities['xmomentum'].centroid_values
     89        ymom = self.domain.quantities['ymomentum'].centroid_values
     90
     91        # Inlet averages
     92        inlet_heights  = stage[inlet_indices]-elevation[inlet_indices]
     93        inlet_areas = areas[inlet_indices]
     94
     95        inlet_water = num.sum(inlet_heights*inlet_areas)
     96
     97        average_inlet_water = inlet_water/self.exchange_areas[0]
     98
     99        # Outlet averages
     100        outlet_heights  = stage[outlet_indices]-elevation[outlet_indices]
     101        outlet_areas = areas[outlet_indices]
     102
     103        outlet_water = num.sum(outlet_heights*outlet_areas)
     104
     105        average_outlet_water = outlet_water/self.exchange_areas[1]
     106
     107
     108        # Transfer
     109        transfer_water = timestep*inlet_water
     110
     111        stage[inlet_indices] = elevation[inlet_indices] + average_inlet_water - transfer_water
     112        xmom[inlet_indices] = 0.0
     113        ymom[inlet_indices] = 0.0
     114
     115
     116        stage[outlet_indices] = elevation[outlet_indices] + average_outlet_water + transfer_water
     117        xmom[outlet_indices] = 0.0
     118        ymom[outlet_indices] = 0.0
     119
     120
     121    def print_stats(self):
     122
     123        print '====================================='
     124        print 'Generic Culvert Operator'
     125        print '====================================='
     126        print "enquiry_gap_factor"
     127        print self.enquiry_gap_factor
     128       
     129        for i in [0,1]:
     130            print '-------------------------------------'
     131            print 'exchange_region %i' % i
     132            print '-------------------------------------'
     133
     134            print 'exchange triangle indices and centres'
     135            print self.exchange_triangle_indices[i]
     136            print self.domain.get_centroid_coordinates()[self.exchange_triangle_indices[i]]
     137       
     138            print 'end_point'
     139            print self.end_points[i]
     140
     141
     142            print 'exchange_polygon'
     143            print self.exchange_polygons[i]
     144
     145            print 'enquiry_point'
     146            print self.enquiry_points[i]
     147
     148        print '====================================='
     149
     150
     151 
     152
     153
     154    def compute_exchange_triangle_indices(self):
     155
     156        # Get boundary (in absolute coordinates)
     157        domain = self.domain
     158        bounding_polygon = domain.get_boundary_polygon()
     159        centroids = domain.get_centroid_coordinates(absolute=True)
     160        self.exchange_triangle_indices = []
     161        self.exchange_areas = []
     162
     163        for i in [0,1]:
     164            exchange_polygon = self.exchange_polygons[i]
     165
     166            # Check that polygon lies within the mesh.
     167            for point in exchange_polygon:
     168                msg = 'Point %s in polygon for forcing term' % str(point)
     169                msg += ' did not fall within the domain boundary.'
     170                assert is_inside_polygon(point, bounding_polygon), msg
     171
     172           
     173            exchange_triangle_indices = inside_polygon(centroids, exchange_polygon)
     174
     175            if len(exchange_triangle_indices) == 0:
     176                region = 'polygon=%s' % (exchange_polygon)
     177                msg = 'No triangles have been identified in '
     178                msg += 'specified region: %s' % region
     179                raise Exception, msg
     180
     181            # Compute exchange area as the sum of areas of triangles identified
     182            # by polygon
     183            exchange_area = 0.0
     184            for j in exchange_triangle_indices:
     185                exchange_area += domain.areas[j]
     186
     187
     188            msg = 'Exchange area %f in culvert' % i
     189            msg += ' has area = %f' % exchange_area
     190            assert exchange_area > 0.0
     191
     192            self.exchange_triangle_indices.append(exchange_triangle_indices)
     193            self.exchange_areas.append(exchange_area)
     194
     195
    78196
    79197    def set_store_hydrograph_discharge(self, filename=None):
     
    123241        # Build exchange polygon and enquiry points 0 and 1
    124242        for i in [0, 1]:
     243            i0 = (2*i-1)
    125244            p0 = self.end_points[i] + w
    126245            p1 = self.end_points[i] - w
    127             p2 = p1 - h
    128             p3 = p0 - h
     246            p2 = p1 + i0*h
     247            p3 = p0 + i0*h
    129248            self.exchange_polygons.append(num.array([p0, p1, p2, p3]))
    130             self.enquiry_points.append(self.end_points[i] + (2*i-1)*gap)
    131 
    132         self.polygon_areas = []
     249            self.enquiry_points.append(self.end_points[i] + i0*gap)
     250
     251
     252
    133253
    134254        # Check that enquiry points are outside exchange polygons
    135255        for i in [0,1]:
    136256            polygon = self.exchange_polygons[i]
    137             # FIXME(SR) should use area of triangles associated with polygon
     257            # FIXME (SR) Probably should calculate the area of all the triangles
     258            # associated with this polygon, as there is likely to be some
     259            # inconsistency between triangles and ploygon
    138260            area = polygon_area(polygon)
    139             self.polygon_areas.append(area)
     261           
    140262
    141263            msg = 'Polygon %s ' %(polygon)
     
    150272
    151273       
    152     def __call__(self, domain):
    153 
    154         # Time stuff
    155         time = domain.get_time()
    156        
    157        
    158         update = False
    159         if self.update_interval is None:
    160             # Use next timestep as has been computed in domain.py       
    161             delta_t = domain.timestep           
    162             update = True
    163         else:   
    164             # Use update interval
    165             delta_t = self.update_interval           
    166             if time - self.last_update > self.update_interval or time == 0.0:
    167                 update = True
    168            
    169         if self.log_filename is not None:       
    170             s = '\nTime = %.2f, delta_t = %f' %(time, delta_t)
    171             log_to_file(self.log_filename, s)
    172                
    173                                
    174         if update is True:
    175             self.compute_rates(delta_t)
    176            
    177    
    178         # Execute flow term for each opening
    179         # This is where Inflow objects are evaluated using the last rate
    180         # that has been calculated
    181         #
    182         # This will take place at every internal timestep and update the domain
    183         for opening in self.openings:
    184             opening(domain)
    185            
    186274
    187275
  • trunk/anuga_core/source/anuga/structures/testing_culvert.py

    r7960 r7962  
    3333width = 5.
    3434
    35 dx = dy = 1           # Resolution: Length of subdivisions on both axes
     35dx = dy = 0.5          # Resolution: Length of subdivisions on both axes
    3636
    3737points, vertices, boundary = rectangular_cross(int(length/dx),
     
    8585
    8686filename=os.path.join(path, 'example_rating_curve.csv')
    87 culvert = Generic_box_culvert(domain,
     87culvert1 = Generic_box_culvert(domain,
    8888                              end_point0=[9.0, 2.5],
    8989                              end_point1=[13.0, 2.5],
     
    9191                              verbose=False)
    9292
    93 #domain.forcing_terms.append(culvert)
    9493
     94#culvert2 = Generic_box_culvert(domain,
     95#                              end_point0=[19.0, 2.5],
     96#                              end_point1=[25.0, 2.5],
     97#                              width=1.00,
     98#                              verbose=False)
     99
     100
     101
     102
     103#print domain.fractional_step_operators
     104
     105#domain.apply_fractional_steps()
    95106
    96107##-----------------------------------------------------------------------
     
    99110
    100111## Inflow based on Flow Depth and Approaching Momentum
    101 #Bi = anuga.Dirichlet_boundary([0.0, 0.0, 0.0])
    102 #Br = anuga.Reflective_boundary(domain)              # Solid reflective wall
     112Bi = anuga.Dirichlet_boundary([1.0, 0.0, 0.0])
     113Br = anuga.Reflective_boundary(domain)              # Solid reflective wall
    103114#Bo = anuga.Dirichlet_boundary([-5, 0, 0])           # Outflow
    104115
     
    111122            #lambda t: [-5*(num.cos(2*pi*(t-4)/20)), 0.0, 0.0])
    112123#domain.set_boundary({'left': Btus, 'right': Btds, 'top': Br, 'bottom': Br})
     124domain.set_boundary({'left': Bi, 'right': Br, 'top': Br, 'bottom': Br})
    113125
    114126
     
    119131#min_delta_w = sys.maxint
    120132#max_delta_w = -min_delta_w
    121 #for t in domain.evolve(yieldstep = 1, finaltime = 25):
     133for t in domain.evolve(yieldstep = 1.0, finaltime = 300):
     134    domain.write_time()
    122135    #delta_w = culvert.inlet.stage - culvert.outlet.stage
    123136   
     
    125138    #if delta_w < min_delta_w: min_delta_w = delta_w           
    126139   
    127     #pass
     140    pass
    128141
    129142## Check that extreme values in rating curve have been exceeded
Note: See TracChangeset for help on using the changeset viewer.