[1280] | 1 | """Class Domain - 2D triangular domains for finite-volume computations of |
---|
| 2 | the shallow water wave equation |
---|
| 3 | |
---|
| 4 | |
---|
| 5 | Copyright 2004 |
---|
| 6 | Ole Nielsen, Stephen Roberts, Duncan Gray, Christopher Zoppou |
---|
| 7 | Geoscience Australia |
---|
| 8 | """ |
---|
| 9 | |
---|
| 10 | from mesh import Mesh |
---|
| 11 | from generic_boundary_conditions import * |
---|
| 12 | import types |
---|
| 13 | |
---|
| 14 | class Domain(Mesh): |
---|
| 15 | |
---|
| 16 | def __init__(self, coordinates, vertices, boundary = None, |
---|
| 17 | conserved_quantities = None, other_quantities = None, |
---|
[1360] | 18 | tagged_elements = None, geo_reference = None, use_inscribed_circle=False): |
---|
[1280] | 19 | |
---|
| 20 | Mesh.__init__(self, coordinates, vertices, boundary, |
---|
[1360] | 21 | tagged_elements, geo_reference, use_inscribed_circle) |
---|
[1280] | 22 | |
---|
| 23 | from Numeric import zeros, Float |
---|
| 24 | from quantity import Quantity, Conserved_quantity |
---|
| 25 | |
---|
| 26 | #List of quantity names entering |
---|
| 27 | #the conservation equations |
---|
| 28 | #(Must be a subset of quantities) |
---|
| 29 | if conserved_quantities is None: |
---|
| 30 | self.conserved_quantities = [] |
---|
| 31 | else: |
---|
| 32 | self.conserved_quantities = conserved_quantities |
---|
| 33 | |
---|
| 34 | if other_quantities is None: |
---|
| 35 | self.other_quantities = [] |
---|
| 36 | else: |
---|
| 37 | self.other_quantities = other_quantities |
---|
| 38 | |
---|
| 39 | |
---|
| 40 | #Build dictionary of Quantity instances keyed by quantity names |
---|
| 41 | self.quantities = {} |
---|
| 42 | |
---|
| 43 | #FIXME: remove later - maybe OK, though.... |
---|
| 44 | for name in self.conserved_quantities: |
---|
| 45 | self.quantities[name] = Conserved_quantity(self) |
---|
| 46 | for name in self.other_quantities: |
---|
| 47 | self.quantities[name] = Quantity(self) |
---|
| 48 | |
---|
| 49 | #Create an empty list for explicit forcing terms |
---|
| 50 | self.forcing_terms = [] |
---|
| 51 | |
---|
| 52 | |
---|
| 53 | #Defaults |
---|
| 54 | from config import max_smallsteps, beta_w, beta_h, epsilon, CFL |
---|
| 55 | self.beta_w = beta_w |
---|
| 56 | self.beta_h = beta_h |
---|
| 57 | self.epsilon = epsilon |
---|
| 58 | |
---|
| 59 | #FIXME: Maybe have separate orders for h-limiter and w-limiter? |
---|
| 60 | #Or maybe get rid of order altogether and use beta_w and beta_h |
---|
| 61 | self.default_order = 1 |
---|
| 62 | self.order = self.default_order |
---|
| 63 | self.smallsteps = 0 |
---|
| 64 | self.max_smallsteps = max_smallsteps |
---|
| 65 | self.number_of_steps = 0 |
---|
| 66 | self.number_of_first_order_steps = 0 |
---|
| 67 | self.CFL = CFL |
---|
| 68 | |
---|
| 69 | #Model time |
---|
| 70 | self.time = 0.0 |
---|
| 71 | self.finaltime = None |
---|
| 72 | self.min_timestep = self.max_timestep = 0.0 |
---|
| 73 | self.starttime = 0 #Physical starttime if any (0 is 1 Jan 1970 00:00:00) |
---|
| 74 | #Origin in UTM coordinates |
---|
| 75 | #FIXME: This should be set if read by a msh file |
---|
| 76 | #self.zone = zone |
---|
| 77 | #self.xllcorner = xllcorner |
---|
| 78 | #self.yllcorner = yllcorner |
---|
| 79 | |
---|
| 80 | |
---|
| 81 | #Checkpointing and storage |
---|
| 82 | from config import default_datadir |
---|
| 83 | self.datadir = default_datadir |
---|
| 84 | self.filename = 'domain' |
---|
| 85 | self.checkpoint = False |
---|
| 86 | |
---|
| 87 | |
---|
| 88 | #Public interface to Domain |
---|
| 89 | def get_conserved_quantities(self, vol_id, vertex=None, edge=None): |
---|
| 90 | """Get conserved quantities at volume vol_id |
---|
| 91 | |
---|
| 92 | If vertex is specified use it as index for vertex values |
---|
| 93 | If edge is specified use it as index for edge values |
---|
| 94 | If neither are specified use centroid values |
---|
| 95 | If both are specified an exeception is raised |
---|
| 96 | |
---|
| 97 | Return value: Vector of length == number_of_conserved quantities |
---|
| 98 | |
---|
| 99 | """ |
---|
| 100 | |
---|
| 101 | from Numeric import zeros, Float |
---|
| 102 | |
---|
| 103 | if not (vertex is None or edge is None): |
---|
| 104 | msg = 'Values for both vertex and edge was specified.' |
---|
| 105 | msg += 'Only one (or none) is allowed.' |
---|
| 106 | raise msg |
---|
| 107 | |
---|
| 108 | q = zeros( len(self.conserved_quantities), Float) |
---|
| 109 | |
---|
| 110 | for i, name in enumerate(self.conserved_quantities): |
---|
| 111 | Q = self.quantities[name] |
---|
| 112 | if vertex is not None: |
---|
| 113 | q[i] = Q.vertex_values[vol_id, vertex] |
---|
| 114 | elif edge is not None: |
---|
| 115 | q[i] = Q.edge_values[vol_id, edge] |
---|
| 116 | else: |
---|
| 117 | q[i] = Q.centroid_values[vol_id] |
---|
| 118 | |
---|
| 119 | return q |
---|
| 120 | |
---|
| 121 | |
---|
| 122 | def set_quantity_vertices_dict(self, quantity_dict): |
---|
| 123 | """Set values for named quantities. |
---|
| 124 | The index is the quantity |
---|
| 125 | |
---|
| 126 | name: Name of quantity |
---|
| 127 | X: Compatible list, Numeric array, const or function (see below) |
---|
| 128 | |
---|
| 129 | The values will be stored in elements following their |
---|
| 130 | internal ordering. |
---|
| 131 | |
---|
| 132 | """ |
---|
| 133 | for key in quantity_dict.keys(): |
---|
| 134 | self.set_quantity(key, quantity_dict[key], location='vertices') |
---|
| 135 | |
---|
| 136 | |
---|
| 137 | def set_quantity(self, name, X, location='vertices', indexes = None): |
---|
| 138 | """Set values for named quantity |
---|
| 139 | |
---|
| 140 | name: Name of quantity |
---|
| 141 | X: Compatible list, Numeric array, const or function (see below) |
---|
| 142 | location: Where values are to be stored. |
---|
[1290] | 143 | Permissible options are: vertices, edges, centroids |
---|
[1280] | 144 | |
---|
[1290] | 145 | In case of location == 'centroids' the dimension values must |
---|
[1280] | 146 | be a list of a Numerical array of length N, N being the number |
---|
| 147 | of elements. Otherwise it must be of dimension Nx3. |
---|
| 148 | |
---|
| 149 | Indexes is the set of element ids that the operation applies to |
---|
| 150 | |
---|
| 151 | The values will be stored in elements following their |
---|
| 152 | internal ordering. |
---|
| 153 | |
---|
[1290] | 154 | #FIXME (Ole): I suggest the following interface |
---|
| 155 | set_quantity(name, X, location, region) |
---|
| 156 | where |
---|
[1280] | 157 | name: Name of quantity |
---|
| 158 | X: |
---|
| 159 | -Compatible list, |
---|
| 160 | -Numeric array, |
---|
| 161 | -const or function (see below) |
---|
| 162 | -another quantity Q or an expression of the form |
---|
| 163 | a*Q+b, where a is a scalar or a compatible array or a quantity |
---|
| 164 | Q is a quantity |
---|
| 165 | b is either a scalar, a quantity or a compatible array |
---|
| 166 | location: Where values are to be stored. |
---|
| 167 | Permissible options are: vertices, edges, centroid |
---|
| 168 | region: Identify subset of triangles. Permissible values are |
---|
| 169 | - tag name (refers to tagged region) |
---|
| 170 | - indices (refers to specific triangles) |
---|
| 171 | - polygon (identifies region) |
---|
| 172 | |
---|
| 173 | This should work for all values of X |
---|
| 174 | |
---|
| 175 | |
---|
| 176 | |
---|
| 177 | """ |
---|
| 178 | |
---|
| 179 | #from quantity import Quantity, Conserved_quantity |
---|
| 180 | |
---|
| 181 | #Create appropriate quantity object |
---|
| 182 | ##if name in self.conserved_quantities: |
---|
| 183 | ## self.quantities[name] = Conserved_quantity(self) |
---|
| 184 | ##else: |
---|
| 185 | ## self.quantities[name] = Quantity(self) |
---|
| 186 | |
---|
| 187 | #Set value |
---|
| 188 | self.quantities[name].set_values(X, location, indexes = indexes) |
---|
| 189 | |
---|
| 190 | |
---|
| 191 | def get_quantity(self, name, location='vertices', indexes = None): |
---|
| 192 | """Get values for named quantity |
---|
| 193 | |
---|
| 194 | name: Name of quantity |
---|
| 195 | |
---|
[1290] | 196 | In case of location == 'centroids' the dimension values must |
---|
[1280] | 197 | be a list of a Numerical array of length N, N being the number |
---|
| 198 | of elements. Otherwise it must be of dimension Nx3. |
---|
| 199 | |
---|
| 200 | Indexes is the set of element ids that the operation applies to. |
---|
| 201 | |
---|
| 202 | The values will be stored in elements following their |
---|
| 203 | internal ordering. |
---|
| 204 | """ |
---|
| 205 | |
---|
| 206 | return self.quantities[name].get_values( location, indexes = indexes) |
---|
| 207 | |
---|
| 208 | |
---|
| 209 | def set_boundary(self, boundary_map): |
---|
| 210 | """Associate boundary objects with tagged boundary segments. |
---|
| 211 | |
---|
| 212 | Input boundary_map is a dictionary of boundary objects keyed |
---|
| 213 | by symbolic tags to matched against tags in the internal dictionary |
---|
| 214 | self.boundary. |
---|
| 215 | |
---|
| 216 | As result one pointer to a boundary object is stored for each vertex |
---|
| 217 | in the list self.boundary_objects. |
---|
| 218 | More entries may point to the same boundary object |
---|
| 219 | |
---|
| 220 | Schematically the mapping is from two dictionaries to one list |
---|
| 221 | where the index is used as pointer to the boundary_values arrays |
---|
| 222 | within each quantity. |
---|
| 223 | |
---|
| 224 | self.boundary: (vol_id, edge_id): tag |
---|
| 225 | boundary_map (input): tag: boundary_object |
---|
| 226 | ---------------------------------------------- |
---|
| 227 | self.boundary_objects: ((vol_id, edge_id), boundary_object) |
---|
| 228 | |
---|
| 229 | |
---|
| 230 | Pre-condition: |
---|
| 231 | self.boundary has been built. |
---|
| 232 | |
---|
| 233 | Post-condition: |
---|
| 234 | self.boundary_objects is built |
---|
| 235 | |
---|
| 236 | If a tag from the domain doesn't appear in the input dictionary an |
---|
| 237 | exception is raised. |
---|
| 238 | However, if a tag is not used to the domain, no error is thrown. |
---|
| 239 | FIXME: This would lead to implementation of a |
---|
| 240 | default boundary condition |
---|
| 241 | |
---|
| 242 | Note: If a segment is listed in the boundary dictionary and if it is |
---|
| 243 | not None, it *will* become a boundary - |
---|
| 244 | even if there is a neighbouring triangle. |
---|
| 245 | This would be the case for internal boundaries |
---|
| 246 | |
---|
| 247 | Boundary objects that are None will be skipped. |
---|
| 248 | |
---|
| 249 | FIXME: If set_boundary is called multiple times and if Boundary |
---|
| 250 | object is changed into None, the neighbour structure will not be |
---|
| 251 | restored!!! |
---|
| 252 | |
---|
| 253 | |
---|
| 254 | """ |
---|
| 255 | |
---|
| 256 | self.boundary_objects = [] |
---|
| 257 | self.boundary_map = boundary_map #Store for use with eg. boundary_stats. |
---|
| 258 | |
---|
| 259 | #FIXME: Try to remove the sorting and fix test_mesh.py |
---|
| 260 | x = self.boundary.keys() |
---|
| 261 | x.sort() |
---|
[1387] | 262 | |
---|
[1280] | 263 | for k, (vol_id, edge_id) in enumerate(x): |
---|
| 264 | tag = self.boundary[ (vol_id, edge_id) ] |
---|
| 265 | |
---|
| 266 | if boundary_map.has_key(tag): |
---|
| 267 | B = boundary_map[tag] |
---|
| 268 | |
---|
| 269 | if B is not None: |
---|
| 270 | self.boundary_objects.append( ((vol_id, edge_id), B) ) |
---|
| 271 | self.neighbours[vol_id, edge_id] = -len(self.boundary_objects) |
---|
| 272 | else: |
---|
| 273 | pass |
---|
| 274 | #FIXME: Check and perhaps fix neighbour structure |
---|
| 275 | |
---|
| 276 | |
---|
| 277 | else: |
---|
| 278 | msg = 'ERROR (domain.py): Tag "%s" has not been ' %tag |
---|
| 279 | msg += 'bound to a boundary object.\n' |
---|
| 280 | msg += 'All boundary tags defined in domain must appear ' |
---|
| 281 | msg += 'in the supplied dictionary.\n' |
---|
| 282 | msg += 'The tags are: %s' %self.get_boundary_tags() |
---|
| 283 | raise msg |
---|
| 284 | |
---|
| 285 | |
---|
| 286 | def set_region(self, functions): |
---|
| 287 | # The order of functions in the list is used. |
---|
| 288 | if type(functions) not in [types.ListType,types.TupleType]: |
---|
| 289 | functions = [functions] |
---|
| 290 | for function in functions: |
---|
| 291 | for tag in self.tagged_elements.keys(): |
---|
| 292 | function(tag, self.tagged_elements[tag], self) |
---|
| 293 | |
---|
| 294 | #Do we need to do this sort of thing? |
---|
| 295 | #self = function(tag, self.tagged_elements[tag], self) |
---|
| 296 | |
---|
| 297 | #MISC |
---|
| 298 | def check_integrity(self): |
---|
| 299 | Mesh.check_integrity(self) |
---|
| 300 | |
---|
| 301 | for quantity in self.conserved_quantities: |
---|
| 302 | msg = 'Conserved quantities must be a subset of all quantities' |
---|
| 303 | assert quantity in self.quantities, msg |
---|
| 304 | |
---|
| 305 | ##assert hasattr(self, 'boundary_objects') |
---|
| 306 | |
---|
| 307 | def write_time(self): |
---|
| 308 | if self.min_timestep == self.max_timestep: |
---|
| 309 | print 'Time = %.4f, delta t = %.8f, steps=%d (%d)'\ |
---|
| 310 | %(self.time, self.min_timestep, self.number_of_steps, |
---|
| 311 | self.number_of_first_order_steps) |
---|
| 312 | elif self.min_timestep > self.max_timestep: |
---|
| 313 | print 'Time = %.4f, steps=%d (%d)'\ |
---|
| 314 | %(self.time, self.number_of_steps, |
---|
| 315 | self.number_of_first_order_steps) |
---|
| 316 | else: |
---|
| 317 | print 'Time = %.4f, delta t in [%.8f, %.8f], steps=%d (%d)'\ |
---|
| 318 | %(self.time, self.min_timestep, |
---|
| 319 | self.max_timestep, self.number_of_steps, |
---|
| 320 | self.number_of_first_order_steps) |
---|
| 321 | |
---|
| 322 | def boundary_stats(self, quantities = None, tag = None): |
---|
| 323 | """Output statistics about boundary forcing |
---|
| 324 | |
---|
| 325 | |
---|
| 326 | """ |
---|
| 327 | |
---|
| 328 | if quantities is None: |
---|
| 329 | quantities = self.conserved_quantities |
---|
| 330 | |
---|
| 331 | |
---|
| 332 | print 'Boundary values at time %.4f:' %self.time |
---|
| 333 | for name in quantities: |
---|
| 334 | q = self.quantities[name] |
---|
| 335 | |
---|
| 336 | if tag is None: |
---|
| 337 | #Take entire boundary |
---|
| 338 | print ' Quantity %s: min = %12.8f, max = %12.8f'\ |
---|
| 339 | %(name, min(q.boundary_values), max(q.boundary_values)) |
---|
| 340 | else: |
---|
| 341 | #Take only boundary associated with tag |
---|
| 342 | maxval = minval = None |
---|
| 343 | for i, ((vol_id, edge_id), B) in enumerate(self.boundary_objects): |
---|
| 344 | if self.boundary[(vol_id, edge_id)] == tag: |
---|
| 345 | v = q.boundary_values[i] |
---|
| 346 | if minval is None or v < minval: minval = v |
---|
| 347 | if maxval is None or v > maxval: maxval = v |
---|
| 348 | |
---|
| 349 | if minval is None or maxval is None: |
---|
| 350 | print 'Sorry no information about tag %s' %tag |
---|
| 351 | else: |
---|
| 352 | print ' Quantity %s, tag %s: min = %12.8f, max = %12.8f'\ |
---|
| 353 | %(name, tag, minval, maxval) |
---|
| 354 | |
---|
| 355 | |
---|
| 356 | |
---|
| 357 | |
---|
| 358 | |
---|
| 359 | |
---|
| 360 | |
---|
| 361 | |
---|
| 362 | |
---|
| 363 | |
---|
| 364 | def get_name(self): |
---|
| 365 | return self.filename |
---|
| 366 | |
---|
| 367 | def set_name(self, name): |
---|
| 368 | self.filename = name |
---|
| 369 | |
---|
| 370 | def get_datadir(self): |
---|
| 371 | return self.datadir |
---|
| 372 | |
---|
| 373 | def set_datadir(self, name): |
---|
| 374 | self.datadir = name |
---|
| 375 | |
---|
| 376 | |
---|
| 377 | |
---|
| 378 | #def set_defaults(self): |
---|
| 379 | # """Set default values for uninitialised quantities. |
---|
| 380 | # Should be overridden or specialised by specific modules |
---|
| 381 | # """# |
---|
| 382 | # |
---|
| 383 | # for name in self.conserved_quantities + self.other_quantities: |
---|
| 384 | # self.set_quantity(name, 0.0) |
---|
| 385 | |
---|
| 386 | |
---|
| 387 | ########################### |
---|
| 388 | #Main components of evolve |
---|
| 389 | |
---|
| 390 | def evolve(self, yieldstep = None, finaltime = None): |
---|
| 391 | """Evolve model from time=0.0 to finaltime yielding results |
---|
| 392 | every yieldstep. |
---|
| 393 | |
---|
| 394 | Internally, smaller timesteps may be taken. |
---|
| 395 | |
---|
| 396 | Evolve is implemented as a generator and is to be called as such, e.g. |
---|
| 397 | |
---|
| 398 | for t in domain.evolve(timestep, yieldstep, finaltime): |
---|
| 399 | <Do something with domain and t> |
---|
| 400 | |
---|
| 401 | """ |
---|
| 402 | |
---|
| 403 | from config import min_timestep, max_timestep, epsilon |
---|
| 404 | |
---|
| 405 | #FIXME: Maybe lump into a larger check prior to evolving |
---|
| 406 | msg = 'Boundary tags must be bound to boundary objects before evolving system, ' |
---|
| 407 | msg += 'e.g. using the method set_boundary.\n' |
---|
| 408 | msg += 'This system has the boundary tags %s ' %self.get_boundary_tags() |
---|
| 409 | assert hasattr(self, 'boundary_objects'), msg |
---|
| 410 | |
---|
| 411 | ##self.set_defaults() |
---|
| 412 | |
---|
| 413 | if yieldstep is None: |
---|
| 414 | yieldstep = max_timestep |
---|
| 415 | else: |
---|
| 416 | yieldstep = float(yieldstep) |
---|
| 417 | |
---|
| 418 | self.order = self.default_order |
---|
| 419 | |
---|
| 420 | |
---|
| 421 | self.yieldtime = 0.0 #Time between 'yields' |
---|
| 422 | |
---|
| 423 | #Initialise interval of timestep sizes (for reporting only) |
---|
| 424 | self.min_timestep = max_timestep |
---|
| 425 | self.max_timestep = min_timestep |
---|
| 426 | self.finaltime = finaltime |
---|
| 427 | self.number_of_steps = 0 |
---|
| 428 | self.number_of_first_order_steps = 0 |
---|
| 429 | |
---|
| 430 | #Initial update of vertex and edge values |
---|
| 431 | self.distribute_to_vertices_and_edges() |
---|
| 432 | |
---|
[1402] | 433 | #Initial update boundary values |
---|
| 434 | self.update_boundary() |
---|
[1280] | 435 | |
---|
| 436 | #Or maybe restore from latest checkpoint |
---|
| 437 | if self.checkpoint is True: |
---|
| 438 | self.goto_latest_checkpoint() |
---|
| 439 | |
---|
| 440 | yield(self.time) #Yield initial values |
---|
| 441 | |
---|
| 442 | while True: |
---|
| 443 | |
---|
| 444 | #Compute fluxes across each element edge |
---|
| 445 | self.compute_fluxes() |
---|
| 446 | |
---|
| 447 | #Update timestep to fit yieldstep and finaltime |
---|
| 448 | self.update_timestep(yieldstep, finaltime) |
---|
| 449 | |
---|
| 450 | #Update conserved quantities |
---|
| 451 | self.update_conserved_quantities() |
---|
| 452 | |
---|
[1402] | 453 | #update ghosts |
---|
| 454 | self.update_ghosts() |
---|
| 455 | |
---|
[1280] | 456 | #Update vertex and edge values |
---|
| 457 | self.distribute_to_vertices_and_edges() |
---|
| 458 | |
---|
[1402] | 459 | #Update boundary values |
---|
| 460 | self.update_boundary() |
---|
| 461 | |
---|
[1280] | 462 | #Update time |
---|
| 463 | self.time += self.timestep |
---|
| 464 | self.yieldtime += self.timestep |
---|
| 465 | self.number_of_steps += 1 |
---|
| 466 | if self.order == 1: |
---|
| 467 | self.number_of_first_order_steps += 1 |
---|
| 468 | |
---|
| 469 | #Yield results |
---|
| 470 | if finaltime is not None and abs(self.time - finaltime) < epsilon: |
---|
| 471 | |
---|
| 472 | #FIXME: There is a rare situation where the |
---|
| 473 | #final time step is stored twice. Can we make a test? |
---|
| 474 | |
---|
| 475 | # Yield final time and stop |
---|
| 476 | yield(self.time) |
---|
| 477 | break |
---|
| 478 | |
---|
| 479 | |
---|
| 480 | if abs(self.yieldtime - yieldstep) < epsilon: |
---|
| 481 | # Yield (intermediate) time and allow inspection of domain |
---|
| 482 | |
---|
| 483 | if self.checkpoint is True: |
---|
| 484 | self.store_checkpoint() |
---|
| 485 | self.delete_old_checkpoints() |
---|
| 486 | |
---|
| 487 | #Pass control on to outer loop for more specific actions |
---|
| 488 | yield(self.time) |
---|
| 489 | |
---|
| 490 | # Reinitialise |
---|
| 491 | self.yieldtime = 0.0 |
---|
| 492 | self.min_timestep = max_timestep |
---|
| 493 | self.max_timestep = min_timestep |
---|
| 494 | self.number_of_steps = 0 |
---|
| 495 | self.number_of_first_order_steps = 0 |
---|
| 496 | |
---|
| 497 | |
---|
| 498 | def evolve_to_end(self, finaltime = 1.0): |
---|
| 499 | """Iterate evolve all the way to the end |
---|
| 500 | """ |
---|
| 501 | |
---|
| 502 | for _ in self.evolve(yieldstep=None, finaltime=finaltime): |
---|
| 503 | pass |
---|
| 504 | |
---|
| 505 | |
---|
| 506 | |
---|
| 507 | def update_boundary(self): |
---|
| 508 | """Go through list of boundary objects and update boundary values |
---|
| 509 | for all conserved quantities on boundary. |
---|
| 510 | """ |
---|
| 511 | |
---|
| 512 | #FIXME: Update only those that change (if that can be worked out) |
---|
| 513 | for i, ((vol_id, edge_id), B) in enumerate(self.boundary_objects): |
---|
| 514 | q = B.evaluate(vol_id, edge_id) |
---|
| 515 | |
---|
| 516 | for j, name in enumerate(self.conserved_quantities): |
---|
| 517 | Q = self.quantities[name] |
---|
| 518 | Q.boundary_values[i] = q[j] |
---|
| 519 | |
---|
| 520 | |
---|
| 521 | def compute_fluxes(self): |
---|
| 522 | msg = 'Method compute_fluxes must be overridden by Domain subclass' |
---|
| 523 | raise msg |
---|
| 524 | |
---|
| 525 | |
---|
| 526 | def update_timestep(self, yieldstep, finaltime): |
---|
| 527 | |
---|
| 528 | from config import min_timestep |
---|
| 529 | |
---|
| 530 | timestep = self.timestep |
---|
| 531 | |
---|
| 532 | #Record maximal and minimal values of timestep for reporting |
---|
| 533 | self.max_timestep = max(timestep, self.max_timestep) |
---|
| 534 | self.min_timestep = min(timestep, self.min_timestep) |
---|
| 535 | |
---|
| 536 | #Protect against degenerate time steps |
---|
| 537 | if timestep < min_timestep: |
---|
| 538 | |
---|
| 539 | #Number of consecutive small steps taken b4 taking action |
---|
| 540 | self.smallsteps += 1 |
---|
| 541 | |
---|
| 542 | if self.smallsteps > self.max_smallsteps: |
---|
| 543 | self.smallsteps = 0 #Reset |
---|
| 544 | |
---|
| 545 | if self.order == 1: |
---|
| 546 | msg = 'WARNING: Too small timestep %.16f reached '\ |
---|
| 547 | %timestep |
---|
| 548 | msg += 'even after %d steps of 1 order scheme'\ |
---|
| 549 | %self.max_smallsteps |
---|
| 550 | print msg |
---|
| 551 | timestep = min_timestep #Try enforcing min_step |
---|
| 552 | |
---|
| 553 | #raise msg |
---|
| 554 | else: |
---|
| 555 | #Try to overcome situation by switching to 1 order |
---|
| 556 | self.order = 1 |
---|
| 557 | |
---|
| 558 | else: |
---|
| 559 | self.smallsteps = 0 |
---|
| 560 | if self.order == 1 and self.default_order == 2: |
---|
| 561 | self.order = 2 |
---|
| 562 | |
---|
| 563 | |
---|
| 564 | #Ensure that final time is not exceeded |
---|
| 565 | if finaltime is not None and self.time + timestep > finaltime: |
---|
| 566 | timestep = finaltime-self.time |
---|
| 567 | |
---|
| 568 | #Ensure that model time is aligned with yieldsteps |
---|
| 569 | if self.yieldtime + timestep > yieldstep: |
---|
| 570 | timestep = yieldstep-self.yieldtime |
---|
| 571 | |
---|
| 572 | self.timestep = timestep |
---|
| 573 | |
---|
| 574 | |
---|
| 575 | |
---|
| 576 | def compute_forcing_terms(self): |
---|
| 577 | """If there are any forcing functions driving the system |
---|
| 578 | they should be defined in Domain subclass and appended to |
---|
| 579 | the list self.forcing_terms |
---|
| 580 | """ |
---|
| 581 | |
---|
| 582 | for f in self.forcing_terms: |
---|
| 583 | f(self) |
---|
| 584 | |
---|
| 585 | |
---|
| 586 | |
---|
| 587 | def update_conserved_quantities(self): |
---|
| 588 | """Update vectors of conserved quantities using previously |
---|
| 589 | computed fluxes specified forcing functions. |
---|
| 590 | """ |
---|
| 591 | |
---|
| 592 | from Numeric import ones, sum, equal, Float |
---|
| 593 | |
---|
| 594 | N = self.number_of_elements |
---|
| 595 | d = len(self.conserved_quantities) |
---|
| 596 | |
---|
| 597 | timestep = self.timestep |
---|
| 598 | |
---|
| 599 | #Compute forcing terms |
---|
| 600 | self.compute_forcing_terms() |
---|
| 601 | |
---|
| 602 | #Update conserved_quantities |
---|
| 603 | for name in self.conserved_quantities: |
---|
| 604 | Q = self.quantities[name] |
---|
| 605 | Q.update(timestep) |
---|
| 606 | |
---|
| 607 | #Clean up |
---|
| 608 | #Note that Q.explicit_update is reset by compute_fluxes |
---|
| 609 | Q.semi_implicit_update[:] = 0.0 |
---|
| 610 | |
---|
[1402] | 611 | def update_ghosts(self): |
---|
| 612 | pass |
---|
[1280] | 613 | |
---|
| 614 | |
---|
| 615 | def distribute_to_vertices_and_edges(self): |
---|
| 616 | """Extrapolate conserved quantities from centroid to |
---|
| 617 | vertices and edge-midpoints for each volume |
---|
| 618 | |
---|
| 619 | Default implementation is straight first order, |
---|
| 620 | i.e. constant values throughout each element and |
---|
| 621 | no reference to non-conserved quantities. |
---|
| 622 | """ |
---|
| 623 | |
---|
| 624 | for name in self.conserved_quantities: |
---|
| 625 | Q = self.quantities[name] |
---|
| 626 | if self.order == 1: |
---|
| 627 | Q.extrapolate_first_order() |
---|
| 628 | elif self.order == 2: |
---|
| 629 | Q.extrapolate_second_order() |
---|
| 630 | Q.limit() |
---|
| 631 | else: |
---|
| 632 | raise 'Unknown order' |
---|
| 633 | Q.interpolate_from_vertices_to_edges() |
---|
| 634 | |
---|
| 635 | |
---|
| 636 | |
---|
| 637 | ############################################## |
---|
| 638 | #Initialise module |
---|
| 639 | |
---|
| 640 | #Optimisation with psyco |
---|
| 641 | from config import use_psyco |
---|
| 642 | if use_psyco: |
---|
| 643 | try: |
---|
| 644 | import psyco |
---|
| 645 | except: |
---|
| 646 | import os |
---|
| 647 | if os.name == 'posix' and os.uname()[4] == 'x86_64': |
---|
| 648 | pass |
---|
| 649 | #Psyco isn't supported on 64 bit systems, but it doesn't matter |
---|
| 650 | else: |
---|
| 651 | msg = 'WARNING: psyco (speedup) could not import'+\ |
---|
| 652 | ', you may want to consider installing it' |
---|
| 653 | print msg |
---|
| 654 | else: |
---|
| 655 | psyco.bind(Domain.update_boundary) |
---|
| 656 | #psyco.bind(Domain.update_timestep) #Not worth it |
---|
| 657 | psyco.bind(Domain.update_conserved_quantities) |
---|
| 658 | psyco.bind(Domain.distribute_to_vertices_and_edges) |
---|
| 659 | |
---|
| 660 | |
---|
| 661 | if __name__ == "__main__": |
---|
| 662 | pass |
---|