[8004] | 1 | import anuga |
---|
[7993] | 2 | import numpy as num |
---|
| 3 | import math |
---|
| 4 | import inlet |
---|
| 5 | |
---|
[8018] | 6 | from anuga.utilities.system_tools import log_to_file |
---|
| 7 | |
---|
| 8 | |
---|
[7993] | 9 | class Structure_operator: |
---|
[7995] | 10 | """Structure Operator - transfer water from one rectangular box to another. |
---|
[7993] | 11 | Sets up the geometry of problem |
---|
| 12 | |
---|
[8008] | 13 | This is the base class for structures (culverts, pipes, bridges etc). Inherit from this class (and overwrite |
---|
| 14 | discharge_routine method for specific subclasses) |
---|
[7993] | 15 | |
---|
| 16 | Input: Two points, pipe_size (either diameter or width, height), |
---|
| 17 | mannings_rougness, |
---|
| 18 | """ |
---|
| 19 | |
---|
[8018] | 20 | counter = 0 |
---|
| 21 | |
---|
[7993] | 22 | def __init__(self, |
---|
| 23 | domain, |
---|
| 24 | end_point0, |
---|
| 25 | end_point1, |
---|
| 26 | width, |
---|
| 27 | height, |
---|
| 28 | apron, |
---|
| 29 | manning, |
---|
| 30 | enquiry_gap, |
---|
[7996] | 31 | description, |
---|
[8018] | 32 | label, |
---|
| 33 | structure_type, |
---|
| 34 | logging, |
---|
[7993] | 35 | verbose): |
---|
| 36 | |
---|
| 37 | self.domain = domain |
---|
| 38 | self.domain.set_fractional_step_operator(self) |
---|
| 39 | self.end_points = [end_point0, end_point1] |
---|
[8018] | 40 | |
---|
[7993] | 41 | |
---|
[8018] | 42 | |
---|
[7993] | 43 | if height is None: |
---|
| 44 | height = width |
---|
| 45 | |
---|
| 46 | if apron is None: |
---|
| 47 | apron = width |
---|
| 48 | |
---|
| 49 | self.width = width |
---|
| 50 | self.height = height |
---|
| 51 | self.apron = apron |
---|
| 52 | self.manning = manning |
---|
| 53 | self.enquiry_gap = enquiry_gap |
---|
[8018] | 54 | |
---|
| 55 | if description == None: |
---|
| 56 | self.description = ' ' |
---|
| 57 | else: |
---|
| 58 | self.description = description |
---|
| 59 | |
---|
| 60 | |
---|
| 61 | if label == None: |
---|
| 62 | self.label = "structure_%g" % Structure_operator.counter |
---|
| 63 | else: |
---|
| 64 | self.label = label |
---|
| 65 | print label |
---|
| 66 | |
---|
| 67 | if structure_type == None: |
---|
| 68 | self.structure_type = 'generic structure' |
---|
| 69 | else: |
---|
| 70 | self.structure_type = structure_type |
---|
| 71 | |
---|
[7993] | 72 | self.verbose = verbose |
---|
[7995] | 73 | |
---|
[8018] | 74 | |
---|
| 75 | |
---|
| 76 | # Keep count of structures |
---|
| 77 | Structure_operator.counter += 1 |
---|
| 78 | |
---|
| 79 | # Slots for recording current statistics |
---|
[7995] | 80 | self.discharge = 0.0 |
---|
| 81 | self.velocity = 0.0 |
---|
[7996] | 82 | self.delta_total_energy = 0.0 |
---|
| 83 | self.driving_energy = 0.0 |
---|
[7993] | 84 | |
---|
| 85 | self.__create_exchange_polygons() |
---|
| 86 | |
---|
| 87 | self.inlets = [] |
---|
| 88 | polygon0 = self.inlet_polygons[0] |
---|
| 89 | enquiry_point0 = self.inlet_equiry_points[0] |
---|
| 90 | outward_vector0 = self.culvert_vector |
---|
| 91 | self.inlets.append(inlet.Inlet(self.domain, polygon0, enquiry_point0, outward_vector0)) |
---|
| 92 | |
---|
| 93 | polygon1 = self.inlet_polygons[1] |
---|
| 94 | exchange_polygon1 = self.inlet_equiry_points[1] |
---|
| 95 | outward_vector1 = - self.culvert_vector |
---|
| 96 | self.inlets.append(inlet.Inlet(self.domain, polygon1, exchange_polygon1, outward_vector1)) |
---|
| 97 | |
---|
[8018] | 98 | self.set_logging(logging) |
---|
[8008] | 99 | |
---|
[7993] | 100 | def __call__(self): |
---|
| 101 | |
---|
[8008] | 102 | timestep = self.domain.get_timestep() |
---|
| 103 | |
---|
| 104 | self.__determine_inflow_outflow() |
---|
| 105 | |
---|
| 106 | Q, barrel_speed, outlet_depth = self.discharge_routine() |
---|
[7993] | 107 | |
---|
[8008] | 108 | old_inflow_height = self.inflow.get_average_height() |
---|
| 109 | old_inflow_xmom = self.inflow.get_average_xmom() |
---|
| 110 | old_inflow_ymom = self.inflow.get_average_ymom() |
---|
| 111 | |
---|
| 112 | if old_inflow_height > 0.0 : |
---|
| 113 | Q_star = Q/old_inflow_height |
---|
| 114 | else: |
---|
| 115 | Q_star = 0.0 |
---|
| 116 | |
---|
| 117 | factor = 1.0/(1.0 + Q_star*timestep/self.inflow.get_area()) |
---|
| 118 | |
---|
| 119 | new_inflow_height = old_inflow_height*factor |
---|
| 120 | new_inflow_xmom = old_inflow_xmom*factor |
---|
| 121 | new_inflow_ymom = old_inflow_ymom*factor |
---|
| 122 | |
---|
| 123 | self.inflow.set_heights(new_inflow_height) |
---|
| 124 | |
---|
| 125 | #inflow.set_xmoms(Q/inflow.get_area()) |
---|
| 126 | #inflow.set_ymoms(0.0) |
---|
| 127 | |
---|
| 128 | self.inflow.set_xmoms(new_inflow_xmom) |
---|
| 129 | self.inflow.set_ymoms(new_inflow_ymom) |
---|
| 130 | |
---|
| 131 | loss = (old_inflow_height - new_inflow_height)*self.inflow.get_area() |
---|
| 132 | |
---|
| 133 | # set outflow |
---|
| 134 | if old_inflow_height > 0.0 : |
---|
| 135 | timestep_star = timestep*new_inflow_height/old_inflow_height |
---|
| 136 | else: |
---|
| 137 | timestep_star = 0.0 |
---|
| 138 | |
---|
| 139 | outflow_extra_height = Q*timestep_star/self.outflow.get_area() |
---|
| 140 | outflow_direction = - self.outflow.outward_culvert_vector |
---|
| 141 | outflow_extra_momentum = outflow_extra_height*barrel_speed*outflow_direction |
---|
| 142 | |
---|
| 143 | gain = outflow_extra_height*self.outflow.get_area() |
---|
| 144 | |
---|
| 145 | #print Q, Q*timestep, barrel_speed, outlet_depth, Qstar, factor, timestep_star |
---|
| 146 | #print ' ', loss, gain |
---|
| 147 | |
---|
| 148 | # Stats |
---|
| 149 | self.discharge = Q#outflow_extra_height*self.outflow.get_area()/timestep |
---|
| 150 | self.velocity = barrel_speed#self.discharge/outlet_depth/self.width |
---|
| 151 | |
---|
| 152 | new_outflow_height = self.outflow.get_average_height() + outflow_extra_height |
---|
| 153 | |
---|
| 154 | if self.use_momentum_jet : |
---|
| 155 | # FIXME (SR) Review momentum to account for possible hydraulic jumps at outlet |
---|
| 156 | #new_outflow_xmom = outflow.get_average_xmom() + outflow_extra_momentum[0] |
---|
| 157 | #new_outflow_ymom = outflow.get_average_ymom() + outflow_extra_momentum[1] |
---|
| 158 | |
---|
| 159 | new_outflow_xmom = barrel_speed*new_outflow_height*outflow_direction[0] |
---|
| 160 | new_outflow_ymom = barrel_speed*new_outflow_height*outflow_direction[1] |
---|
| 161 | |
---|
| 162 | else: |
---|
| 163 | #new_outflow_xmom = outflow.get_average_xmom() |
---|
| 164 | #new_outflow_ymom = outflow.get_average_ymom() |
---|
| 165 | |
---|
| 166 | new_outflow_xmom = 0.0 |
---|
| 167 | new_outflow_ymom = 0.0 |
---|
| 168 | |
---|
| 169 | self.outflow.set_heights(new_outflow_height) |
---|
| 170 | self.outflow.set_xmoms(new_outflow_xmom) |
---|
| 171 | self.outflow.set_ymoms(new_outflow_ymom) |
---|
| 172 | |
---|
| 173 | |
---|
| 174 | def __determine_inflow_outflow(self): |
---|
| 175 | # Determine flow direction based on total energy difference |
---|
| 176 | |
---|
| 177 | if self.use_velocity_head: |
---|
| 178 | self.delta_total_energy = self.inlets[0].get_enquiry_total_energy() - self.inlets[1].get_enquiry_total_energy() |
---|
| 179 | else: |
---|
| 180 | self.delta_total_energy = self.inlets[0].get_enquiry_stage() - self.inlets[1].get_enquiry_stage() |
---|
| 181 | |
---|
| 182 | |
---|
| 183 | self.inflow = self.inlets[0] |
---|
| 184 | self.outflow = self.inlets[1] |
---|
| 185 | |
---|
| 186 | |
---|
| 187 | if self.delta_total_energy < 0: |
---|
| 188 | self.inflow = self.inlets[1] |
---|
| 189 | self.outflow = self.inlets[0] |
---|
| 190 | self.delta_total_energy = -self.delta_total_energy |
---|
| 191 | |
---|
| 192 | |
---|
[7993] | 193 | def __create_exchange_polygons(self): |
---|
| 194 | |
---|
| 195 | """Create polygons at the end of a culvert inlet and outlet. |
---|
| 196 | At either end two polygons will be created; one for the actual flow to pass through and one a little further away |
---|
| 197 | for enquiring the total energy at both ends of the culvert and transferring flow. |
---|
| 198 | """ |
---|
| 199 | |
---|
| 200 | # Calculate geometry |
---|
| 201 | x0, y0 = self.end_points[0] |
---|
| 202 | x1, y1 = self.end_points[1] |
---|
| 203 | |
---|
| 204 | dx = x1 - x0 |
---|
| 205 | dy = y1 - y0 |
---|
| 206 | |
---|
| 207 | self.culvert_vector = num.array([dx, dy]) |
---|
| 208 | self.culvert_length = math.sqrt(num.sum(self.culvert_vector**2)) |
---|
| 209 | assert self.culvert_length > 0.0, 'The length of culvert is less than 0' |
---|
| 210 | |
---|
| 211 | # Unit direction vector and normal |
---|
| 212 | self.culvert_vector /= self.culvert_length # Unit vector in culvert direction |
---|
| 213 | self.culvert_normal = num.array([-dy, dx])/self.culvert_length # Normal vector |
---|
| 214 | |
---|
| 215 | # Short hands |
---|
| 216 | w = 0.5*self.width*self.culvert_normal # Perpendicular vector of 1/2 width |
---|
| 217 | h = self.apron*self.culvert_vector # Vector of length=height in the |
---|
| 218 | # direction of the culvert |
---|
| 219 | |
---|
| 220 | gap = (1 + self.enquiry_gap)*h |
---|
| 221 | |
---|
| 222 | self.inlet_polygons = [] |
---|
| 223 | self.inlet_equiry_points = [] |
---|
| 224 | |
---|
| 225 | # Build exchange polygon and enquiry point |
---|
| 226 | for i in [0, 1]: |
---|
| 227 | i0 = (2*i-1) |
---|
| 228 | p0 = self.end_points[i] + w |
---|
| 229 | p1 = self.end_points[i] - w |
---|
| 230 | p2 = p1 + i0*h |
---|
| 231 | p3 = p0 + i0*h |
---|
| 232 | ep = self.end_points[i] + i0*gap |
---|
| 233 | |
---|
| 234 | self.inlet_polygons.append(num.array([p0, p1, p2, p3])) |
---|
| 235 | self.inlet_equiry_points.append(ep) |
---|
| 236 | |
---|
| 237 | # Check that enquiry points are outside inlet polygons |
---|
| 238 | for i in [0,1]: |
---|
| 239 | polygon = self.inlet_polygons[i] |
---|
| 240 | ep = self.inlet_equiry_points[i] |
---|
| 241 | |
---|
[8004] | 242 | area = anuga.polygon_area(polygon) |
---|
[7993] | 243 | |
---|
| 244 | msg = 'Polygon %s ' %(polygon) |
---|
| 245 | msg += ' has area = %f' % area |
---|
| 246 | assert area > 0.0, msg |
---|
| 247 | |
---|
| 248 | msg = 'Enquiry point falls inside an exchange polygon.' |
---|
[8004] | 249 | assert not anuga.inside_polygon(ep, polygon), msg |
---|
[8008] | 250 | |
---|
[7993] | 251 | |
---|
[8008] | 252 | def discharge_routine(self): |
---|
| 253 | |
---|
| 254 | pass |
---|
[7993] | 255 | |
---|
| 256 | |
---|
[8018] | 257 | def structure_statistics(self): |
---|
[7993] | 258 | |
---|
| 259 | |
---|
[8018] | 260 | message = '=====================================\n' |
---|
| 261 | message += 'Structure Operator: %s\n' % self.label |
---|
| 262 | message += '=====================================\n' |
---|
[7993] | 263 | |
---|
[8018] | 264 | message += 'Structure Type: %s\n' % self.structure_type |
---|
| 265 | |
---|
| 266 | message += 'Description\n' |
---|
| 267 | message += '%s' % self.description |
---|
| 268 | message += '\n' |
---|
[7993] | 269 | |
---|
| 270 | for i, inlet in enumerate(self.inlets): |
---|
[8018] | 271 | message += '-------------------------------------\n' |
---|
| 272 | message += 'Inlet %i\n' % i |
---|
| 273 | message += '-------------------------------------\n' |
---|
[7993] | 274 | |
---|
[8018] | 275 | message += 'inlet triangle indices and centres\n' |
---|
| 276 | message += '%s' % inlet.triangle_indices |
---|
| 277 | message += '\n' |
---|
| 278 | |
---|
| 279 | message += '%s' % self.domain.get_centroid_coordinates()[inlet.triangle_indices] |
---|
| 280 | message += '\n' |
---|
[7993] | 281 | |
---|
[8018] | 282 | message += 'polygon\n' |
---|
| 283 | message += '%s' % inlet.polygon |
---|
| 284 | message += '\n' |
---|
[7993] | 285 | |
---|
[8018] | 286 | message += '=====================================\n' |
---|
[7995] | 287 | |
---|
[8018] | 288 | return message |
---|
[7995] | 289 | |
---|
[8018] | 290 | |
---|
| 291 | def print_structure_statistics(self): |
---|
| 292 | |
---|
| 293 | print self.structure_statistics() |
---|
| 294 | |
---|
| 295 | |
---|
| 296 | def print_timestepping_statistics(self): |
---|
| 297 | |
---|
[7995] | 298 | message = '---------------------------\n' |
---|
[8018] | 299 | message += 'Structure report for %s:\n' % self.label |
---|
[7995] | 300 | message += '--------------------------\n' |
---|
[8018] | 301 | message += 'Type: %s\n' % self.structure_type |
---|
[7995] | 302 | message += 'Discharge [m^3/s]: %.2f\n' % self.discharge |
---|
| 303 | message += 'Velocity [m/s]: %.2f\n' % self.velocity |
---|
[7996] | 304 | message += 'Inlet Driving Energy %.2f\n' % self.driving_energy |
---|
| 305 | message += 'delta total energy %.2f\n' % self.delta_total_energy |
---|
[7995] | 306 | |
---|
[8018] | 307 | print message |
---|
| 308 | |
---|
| 309 | |
---|
| 310 | def set_logging(self, flag=True): |
---|
| 311 | |
---|
| 312 | self.logging = flag |
---|
| 313 | |
---|
| 314 | # If flag is true open file with mode = "w" to form a clean file for logging |
---|
| 315 | if self.logging: |
---|
| 316 | self.log_filename = self.label + '.log' |
---|
| 317 | log_to_file(self.log_filename, self.structure_statistics(), mode='w') |
---|
| 318 | log_to_file(self.log_filename, 'time,discharge,velocity,driving_energy,delta_total_energy') |
---|
| 319 | |
---|
| 320 | #log_to_file(self.log_filename, self.culvert_type) |
---|
| 321 | |
---|
| 322 | |
---|
| 323 | def timestepping_statistics(self): |
---|
| 324 | |
---|
| 325 | message = '%.5f, ' % self.domain.get_time() |
---|
| 326 | message += '%.5f, ' % self.discharge |
---|
| 327 | message += '%.5f, ' % self.velocity |
---|
| 328 | message += '%.5f, ' % self.driving_energy |
---|
| 329 | message += '%.5f' % self.delta_total_energy |
---|
| 330 | |
---|
[7995] | 331 | return message |
---|
| 332 | |
---|
[8018] | 333 | def log_timestepping_statistics(self): |
---|
[8008] | 334 | |
---|
[8018] | 335 | if self.logging: |
---|
| 336 | log_to_file(self.log_filename, self.timestepping_statistics()) |
---|
| 337 | |
---|
| 338 | |
---|
[7993] | 339 | def get_inlets(self): |
---|
| 340 | |
---|
| 341 | return self.inlets |
---|
| 342 | |
---|
| 343 | |
---|
| 344 | def get_culvert_length(self): |
---|
| 345 | |
---|
| 346 | return self.culvert_length |
---|
| 347 | |
---|
| 348 | |
---|
| 349 | def get_culvert_width(self): |
---|
| 350 | |
---|
| 351 | return self.width |
---|
| 352 | |
---|
| 353 | |
---|
[7998] | 354 | def get_culvert_diameter(self): |
---|
| 355 | |
---|
| 356 | return self.width |
---|
| 357 | |
---|
| 358 | |
---|
[7993] | 359 | def get_culvert_height(self): |
---|
| 360 | |
---|
| 361 | return self.height |
---|
| 362 | |
---|
| 363 | |
---|
| 364 | def get_culvert_apron(self): |
---|
| 365 | |
---|
| 366 | return self.apron |
---|