Changeset 7962
 Timestamp:
 Aug 20, 2010, 5:27:03 PM (12 years ago)
 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 46 46 self.domain = domain 47 47 48 self.domain.set_fractional_step_operator(self) 49 48 50 self.end_points= [end_point0, end_point1] 49 51 self.enquiry_gap_factor = enquiry_gap_factor … … 62 64 self.compute_enquiry_indices() 63 65 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 78 196 79 197 def set_store_hydrograph_discharge(self, filename=None): … … 123 241 # Build exchange polygon and enquiry points 0 and 1 124 242 for i in [0, 1]: 243 i0 = (2*i1) 125 244 p0 = self.end_points[i] + w 126 245 p1 = self.end_points[i]  w 127 p2 = p1 h128 p3 = p0 h246 p2 = p1 + i0*h 247 p3 = p0 + i0*h 129 248 self.exchange_polygons.append(num.array([p0, p1, p2, p3])) 130 self.enquiry_points.append(self.end_points[i] + (2*i1)*gap) 131 132 self.polygon_areas = [] 249 self.enquiry_points.append(self.end_points[i] + i0*gap) 250 251 252 133 253 134 254 # Check that enquiry points are outside exchange polygons 135 255 for i in [0,1]: 136 256 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 138 260 area = polygon_area(polygon) 139 self.polygon_areas.append(area)261 140 262 141 263 msg = 'Polygon %s ' %(polygon) … … 150 272 151 273 152 def __call__(self, domain):153 154 # Time stuff155 time = domain.get_time()156 157 158 update = False159 if self.update_interval is None:160 # Use next timestep as has been computed in domain.py161 delta_t = domain.timestep162 update = True163 else:164 # Use update interval165 delta_t = self.update_interval166 if time  self.last_update > self.update_interval or time == 0.0:167 update = True168 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 opening179 # This is where Inflow objects are evaluated using the last rate180 # that has been calculated181 #182 # This will take place at every internal timestep and update the domain183 for opening in self.openings:184 opening(domain)185 186 274 187 275 
trunk/anuga_core/source/anuga/structures/testing_culvert.py
r7960 r7962 33 33 width = 5. 34 34 35 dx = dy = 1# Resolution: Length of subdivisions on both axes35 dx = dy = 0.5 # Resolution: Length of subdivisions on both axes 36 36 37 37 points, vertices, boundary = rectangular_cross(int(length/dx), … … 85 85 86 86 filename=os.path.join(path, 'example_rating_curve.csv') 87 culvert = Generic_box_culvert(domain,87 culvert1 = Generic_box_culvert(domain, 88 88 end_point0=[9.0, 2.5], 89 89 end_point1=[13.0, 2.5], … … 91 91 verbose=False) 92 92 93 #domain.forcing_terms.append(culvert)94 93 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() 95 106 96 107 ## … … 99 110 100 111 ## 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 wall112 Bi = anuga.Dirichlet_boundary([1.0, 0.0, 0.0]) 113 Br = anuga.Reflective_boundary(domain) # Solid reflective wall 103 114 #Bo = anuga.Dirichlet_boundary([5, 0, 0]) # Outflow 104 115 … … 111 122 #lambda t: [5*(num.cos(2*pi*(t4)/20)), 0.0, 0.0]) 112 123 #domain.set_boundary({'left': Btus, 'right': Btds, 'top': Br, 'bottom': Br}) 124 domain.set_boundary({'left': Bi, 'right': Br, 'top': Br, 'bottom': Br}) 113 125 114 126 … … 119 131 #min_delta_w = sys.maxint 120 132 #max_delta_w = min_delta_w 121 #for t in domain.evolve(yieldstep = 1, finaltime = 25): 133 for t in domain.evolve(yieldstep = 1.0, finaltime = 300): 134 domain.write_time() 122 135 #delta_w = culvert.inlet.stage  culvert.outlet.stage 123 136 … … 125 138 #if delta_w < min_delta_w: min_delta_w = delta_w 126 139 127 #pass140 pass 128 141 129 142 ## Check that extreme values in rating curve have been exceeded
Note: See TracChangeset
for help on using the changeset viewer.