| 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 | |
|---|
| 13 | class Domain(Mesh): |
|---|
| 14 | |
|---|
| 15 | def __init__(self, coordinates, vertices, boundary = None, |
|---|
| 16 | conserved_quantities = None, other_quantities = None): |
|---|
| 17 | |
|---|
| 18 | Mesh.__init__(self, coordinates, vertices, boundary) |
|---|
| 19 | |
|---|
| 20 | from Numeric import zeros, Float |
|---|
| 21 | from quantity import Quantity |
|---|
| 22 | |
|---|
| 23 | #List of quantity names entering |
|---|
| 24 | #the conservation equations |
|---|
| 25 | #(Must be a subset of quantities) |
|---|
| 26 | if conserved_quantities is None: |
|---|
| 27 | self.conserved_quantities = [] |
|---|
| 28 | else: |
|---|
| 29 | self.conserved_quantities = conserved_quantities |
|---|
| 30 | |
|---|
| 31 | if other_quantities is None: |
|---|
| 32 | other_quantities = [] |
|---|
| 33 | |
|---|
| 34 | #Build dictionary of Quantity instances keyed by quantity names |
|---|
| 35 | self.quantities = {} |
|---|
| 36 | for name in self.conserved_quantities + other_quantities: |
|---|
| 37 | self.quantities[name] = Quantity(self) |
|---|
| 38 | |
|---|
| 39 | #Defaults |
|---|
| 40 | from config import max_smallsteps, beta |
|---|
| 41 | self.beta = beta |
|---|
| 42 | self.default_order = 1 |
|---|
| 43 | self.order = self.default_order |
|---|
| 44 | self.smallsteps = 0 |
|---|
| 45 | self.max_smallsteps = max_smallsteps |
|---|
| 46 | self.number_of_steps = 0 |
|---|
| 47 | self.number_of_first_order_steps = 0 |
|---|
| 48 | |
|---|
| 49 | #Model time |
|---|
| 50 | self.time = 0.0 |
|---|
| 51 | self.finaltime = None |
|---|
| 52 | self.min_timestep = self.max_timestep = 0.0 |
|---|
| 53 | |
|---|
| 54 | #Checkpointing |
|---|
| 55 | self.filename = 'domain' |
|---|
| 56 | self.checkpoint = False |
|---|
| 57 | |
|---|
| 58 | #Realtime visualisation |
|---|
| 59 | self.visualise = False |
|---|
| 60 | |
|---|
| 61 | #Stored output |
|---|
| 62 | self.store=False |
|---|
| 63 | self.format = 'dat' |
|---|
| 64 | #self.smooth = False |
|---|
| 65 | self.smooth = True |
|---|
| 66 | |
|---|
| 67 | #Reduction operation for get_vertex_values |
|---|
| 68 | #from pytools.stats import mean |
|---|
| 69 | #self.reduction = mean |
|---|
| 70 | self.reduction = min #Looks better near steep slopes |
|---|
| 71 | |
|---|
| 72 | |
|---|
| 73 | #Public interface to Domain |
|---|
| 74 | def get_conserved_quantities(self, vol_id, vertex=None, edge=None): |
|---|
| 75 | """Get conserved quantities at volume vol_id |
|---|
| 76 | |
|---|
| 77 | If vertex is specified use it as index for vertex values |
|---|
| 78 | If edge is specified use it as index for edge values |
|---|
| 79 | If neither are specified use centroid values |
|---|
| 80 | If both are specified an exeception is raised |
|---|
| 81 | |
|---|
| 82 | Return value: Vector of length == number_of_conserved quantities |
|---|
| 83 | |
|---|
| 84 | """ |
|---|
| 85 | |
|---|
| 86 | from Numeric import zeros, Float |
|---|
| 87 | |
|---|
| 88 | if not (vertex is None or edge is None): |
|---|
| 89 | msg = 'Values for both vertex and edge was specified.' |
|---|
| 90 | msg += 'Only one (or none) is allowed.' |
|---|
| 91 | raise msg |
|---|
| 92 | |
|---|
| 93 | q = zeros( len(self.conserved_quantities), Float) |
|---|
| 94 | |
|---|
| 95 | for i, name in enumerate(self.conserved_quantities): |
|---|
| 96 | Q = self.quantities[name] |
|---|
| 97 | if vertex is not None: |
|---|
| 98 | q[i] = Q.vertex_values[vol_id, vertex] |
|---|
| 99 | elif edge is not None: |
|---|
| 100 | q[i] = Q.edge_values[vol_id, edge] |
|---|
| 101 | else: |
|---|
| 102 | q[i] = Q.centroid_values[vol_id] |
|---|
| 103 | |
|---|
| 104 | return q |
|---|
| 105 | |
|---|
| 106 | def set_quantity(self, name, values, location='vertices'): |
|---|
| 107 | """Set values for named quantity |
|---|
| 108 | |
|---|
| 109 | name: Name of quantity |
|---|
| 110 | values: Compatible list or Numeric array (see below) |
|---|
| 111 | location: Where values are to be stored. |
|---|
| 112 | Permissible options are: vertices, edges, centroid |
|---|
| 113 | |
|---|
| 114 | In case of location == 'centroid' the dimension values must |
|---|
| 115 | be a list of a Numerical array of length N, N being the number |
|---|
| 116 | of elements in the mesh. Otherwise it must be of dimension Nx3 |
|---|
| 117 | |
|---|
| 118 | The values will be stored in elements following their |
|---|
| 119 | internal ordering. |
|---|
| 120 | """ |
|---|
| 121 | |
|---|
| 122 | self.quantities[name].set_values(values, location) |
|---|
| 123 | |
|---|
| 124 | |
|---|
| 125 | def set_boundary(self, boundary_map): |
|---|
| 126 | """Associate boundary objects with tagged boundary segments. |
|---|
| 127 | |
|---|
| 128 | Input boundary_map is a dictionary of boundary objects keyed |
|---|
| 129 | by symbolic tags to matched against tags in the internal dictionary |
|---|
| 130 | self.boundary. |
|---|
| 131 | |
|---|
| 132 | As result one pointer to a boundary object is stored for each vertex |
|---|
| 133 | in the list self.boundary_objects. |
|---|
| 134 | More entries may point to the same boundary object |
|---|
| 135 | |
|---|
| 136 | Schematically the mapping is: |
|---|
| 137 | |
|---|
| 138 | self.boundary_segments: k: (vol_id, edge_id) |
|---|
| 139 | self.boundary: (vol_id, edge_id): tag |
|---|
| 140 | boundary_map (input): tag: boundary_object |
|---|
| 141 | ---------------------------------------------- |
|---|
| 142 | self.boundary_objects: k: boundary_object |
|---|
| 143 | |
|---|
| 144 | |
|---|
| 145 | Pre-condition: |
|---|
| 146 | self.boundary and self.boundary_segments have been built. |
|---|
| 147 | |
|---|
| 148 | Post-condition: |
|---|
| 149 | self.boundary_objects is built |
|---|
| 150 | |
|---|
| 151 | If a tag from the domain doesn't appear in the input dictionary an |
|---|
| 152 | exception is raised. |
|---|
| 153 | However, if a tag is not used to the domain, no error is thrown. |
|---|
| 154 | FIXME: This would lead to implementation of a |
|---|
| 155 | default boundary condition |
|---|
| 156 | """ |
|---|
| 157 | |
|---|
| 158 | self.boundary_objects = [] |
|---|
| 159 | for k, (vol_id, edge_id) in enumerate(self.boundary_segments): |
|---|
| 160 | tag = self.boundary[ (vol_id, edge_id) ] |
|---|
| 161 | |
|---|
| 162 | if boundary_map.has_key(tag): |
|---|
| 163 | B = boundary_map[tag] |
|---|
| 164 | self.boundary_objects.append(B) |
|---|
| 165 | |
|---|
| 166 | else: |
|---|
| 167 | msg = 'ERROR (domain.py): Tag "%s" has not been ' %tag |
|---|
| 168 | msg += 'bound to a boundary object.\n' |
|---|
| 169 | msg += 'All boundary tags defined in domain must appear ' |
|---|
| 170 | msg += 'in the supplied dictionary.\n' |
|---|
| 171 | msg += 'The tags are: %s' %self.get_boundary_tags() |
|---|
| 172 | raise msg |
|---|
| 173 | |
|---|
| 174 | |
|---|
| 175 | |
|---|
| 176 | #MISC |
|---|
| 177 | def check_integrity(self): |
|---|
| 178 | Mesh.check_integrity(self) |
|---|
| 179 | |
|---|
| 180 | for quantity in self.conserved_quantities: |
|---|
| 181 | msg = 'Conserved quantities must be a subset of all quantities' |
|---|
| 182 | assert quantity in self.quantities, msg |
|---|
| 183 | |
|---|
| 184 | def write_time(self): |
|---|
| 185 | if self.min_timestep == self.max_timestep: |
|---|
| 186 | print 'Time = %.4f, delta t = %.8f, steps=%d (%d)'\ |
|---|
| 187 | %(self.time, self.min_timestep, self.number_of_steps, |
|---|
| 188 | self.number_of_first_order_steps) |
|---|
| 189 | else: |
|---|
| 190 | print 'Time = %.4f, delta t in [%.8f, %.8f], steps=%d (%d)'\ |
|---|
| 191 | %(self.time, self.min_timestep, |
|---|
| 192 | self.max_timestep, self.number_of_steps, |
|---|
| 193 | self.number_of_first_order_steps) |
|---|
| 194 | |
|---|
| 195 | |
|---|
| 196 | |
|---|
| 197 | |
|---|
| 198 | ########################### |
|---|
| 199 | #Main components of evolve |
|---|
| 200 | |
|---|
| 201 | def evolve(self, yieldstep = None, finaltime = None): |
|---|
| 202 | """Evolve model from time=0.0 to finaltime yielding results |
|---|
| 203 | every yieldstep. |
|---|
| 204 | |
|---|
| 205 | Internally, smaller timesteps may be taken. |
|---|
| 206 | |
|---|
| 207 | Evolve is implemented as a generator and is to be called as such, e.g. |
|---|
| 208 | |
|---|
| 209 | for t in domain.evolve(timestep, yieldstep, finaltime): |
|---|
| 210 | <Do something with domain and t> |
|---|
| 211 | |
|---|
| 212 | """ |
|---|
| 213 | |
|---|
| 214 | #import data_manager |
|---|
| 215 | from config import min_timestep, max_timestep, epsilon |
|---|
| 216 | |
|---|
| 217 | if yieldstep is None: |
|---|
| 218 | yieldstep = max_timestep |
|---|
| 219 | |
|---|
| 220 | |
|---|
| 221 | self.order = self.default_order |
|---|
| 222 | |
|---|
| 223 | |
|---|
| 224 | self.yieldtime = 0.0 #Time between 'yields' |
|---|
| 225 | |
|---|
| 226 | #Initialise interval of timestep sizes (for reporting only) |
|---|
| 227 | self.min_timestep = max_timestep |
|---|
| 228 | self.max_timestep = min_timestep |
|---|
| 229 | self.finaltime = finaltime |
|---|
| 230 | self.number_of_steps = 0 |
|---|
| 231 | self.number_of_first_order_steps = 0 |
|---|
| 232 | |
|---|
| 233 | #Initial update of vertex and edge values |
|---|
| 234 | self.distribute_to_vertices_and_edges() |
|---|
| 235 | |
|---|
| 236 | |
|---|
| 237 | #Or maybe restore from latest checkpoint |
|---|
| 238 | if self.checkpoint is True: |
|---|
| 239 | self.goto_latest_checkpoint() |
|---|
| 240 | |
|---|
| 241 | |
|---|
| 242 | #Store model data, e.g. for visualisation |
|---|
| 243 | if self.store is True and self.time == 0.0: |
|---|
| 244 | self.store_bathymetry() |
|---|
| 245 | self.store_conserved_quantities() |
|---|
| 246 | |
|---|
| 247 | if self.visualise is True and self.time == 0.0: |
|---|
| 248 | import realtime_visualisation as visualise |
|---|
| 249 | visualise.create_surface(self) |
|---|
| 250 | |
|---|
| 251 | |
|---|
| 252 | yield(self.time) #Yield initial values |
|---|
| 253 | |
|---|
| 254 | while True: |
|---|
| 255 | #Update boundary values |
|---|
| 256 | self.update_boundary() |
|---|
| 257 | |
|---|
| 258 | #Compute all fluxes and timestep suitable for all volumes |
|---|
| 259 | self.compute_fluxes() |
|---|
| 260 | |
|---|
| 261 | #Update timestep to fit yieldstep and finaltime |
|---|
| 262 | self.update_timestep(yieldstep, finaltime) |
|---|
| 263 | |
|---|
| 264 | #Update conserved quantities |
|---|
| 265 | self.update_conserved_quantities() |
|---|
| 266 | |
|---|
| 267 | #Update vertex and edge values |
|---|
| 268 | self.distribute_to_vertices_and_edges() |
|---|
| 269 | |
|---|
| 270 | #Update time |
|---|
| 271 | self.time += self.timestep |
|---|
| 272 | self.yieldtime += self.timestep |
|---|
| 273 | self.number_of_steps += 1 |
|---|
| 274 | if self.order == 1: |
|---|
| 275 | self.number_of_first_order_steps += 1 |
|---|
| 276 | |
|---|
| 277 | #Yield results |
|---|
| 278 | if finaltime is not None and abs(self.time - finaltime) < epsilon: |
|---|
| 279 | # Yield final time and allow inspection of domain |
|---|
| 280 | yield(self.time) |
|---|
| 281 | break |
|---|
| 282 | |
|---|
| 283 | if abs(self.yieldtime - yieldstep) < epsilon: |
|---|
| 284 | # Yield (intermediate) time and allow inspection of domain |
|---|
| 285 | |
|---|
| 286 | if self.checkpoint is True: |
|---|
| 287 | self.store_checkpoint() |
|---|
| 288 | self.delete_old_checkpoints() |
|---|
| 289 | |
|---|
| 290 | #Store model data, e.g. for subsequent visualisation |
|---|
| 291 | if self.store is True: |
|---|
| 292 | self.store_conserved_quantities() |
|---|
| 293 | |
|---|
| 294 | #Real time viz |
|---|
| 295 | if self.visualise is True: |
|---|
| 296 | visualise.update(self) |
|---|
| 297 | |
|---|
| 298 | yield(self.time) |
|---|
| 299 | |
|---|
| 300 | # Reinitialise |
|---|
| 301 | self.yieldtime = 0.0 |
|---|
| 302 | self.min_timestep = max_timestep |
|---|
| 303 | self.max_timestep = min_timestep |
|---|
| 304 | self.number_of_steps = 0 |
|---|
| 305 | self.number_of_first_order_steps = 0 |
|---|
| 306 | |
|---|
| 307 | |
|---|
| 308 | def evolve_to_end(self, finaltime = 1.0): |
|---|
| 309 | """Iterate evolve generator all the way to the end |
|---|
| 310 | """ |
|---|
| 311 | |
|---|
| 312 | for _ in self.evolve(yieldstep=None, finaltime=finaltime): |
|---|
| 313 | pass |
|---|
| 314 | |
|---|
| 315 | |
|---|
| 316 | |
|---|
| 317 | def update_boundary(self): |
|---|
| 318 | """Go through list of boundary objects and update boundary values |
|---|
| 319 | for all conserved quantities on boundary. |
|---|
| 320 | """ |
|---|
| 321 | |
|---|
| 322 | #FIXME: Update only those that change (if that can be worked out) |
|---|
| 323 | for i, B in enumerate(self.boundary_objects): |
|---|
| 324 | vol_id, edge_id = self.boundary_segments[i] |
|---|
| 325 | q = B.evaluate(vol_id, edge_id) |
|---|
| 326 | |
|---|
| 327 | for j, name in enumerate(self.conserved_quantities): |
|---|
| 328 | Q = self.quantities[name] |
|---|
| 329 | Q.boundary_values[i] = q[j] |
|---|
| 330 | |
|---|
| 331 | |
|---|
| 332 | def compute_fluxes(self): |
|---|
| 333 | msg = 'Method compute_fluxes must be overridden by Domain subclass' |
|---|
| 334 | raise msg |
|---|
| 335 | |
|---|
| 336 | |
|---|
| 337 | def update_timestep(self, yieldstep, finaltime): |
|---|
| 338 | |
|---|
| 339 | from config import min_timestep |
|---|
| 340 | |
|---|
| 341 | timestep = self.timestep |
|---|
| 342 | |
|---|
| 343 | #Record maximal and minimal values of timestep for reporting |
|---|
| 344 | self.max_timestep = max(timestep, self.max_timestep) |
|---|
| 345 | self.min_timestep = min(timestep, self.min_timestep) |
|---|
| 346 | |
|---|
| 347 | #Protect against degenerate time steps |
|---|
| 348 | if timestep < min_timestep: |
|---|
| 349 | |
|---|
| 350 | #Number of consecutive small steps taken b4 taking action |
|---|
| 351 | self.smallsteps += 1 |
|---|
| 352 | |
|---|
| 353 | if self.smallsteps > self.max_smallsteps: |
|---|
| 354 | self.smallsteps = 0 #Reset |
|---|
| 355 | |
|---|
| 356 | if self.order == 1: |
|---|
| 357 | msg = 'Minimal timestep %.16f reached ' %self.min_timestep |
|---|
| 358 | msg += 'using 1 order scheme' |
|---|
| 359 | |
|---|
| 360 | raise msg |
|---|
| 361 | else: |
|---|
| 362 | #Try to overcome situation by switching to 1 order |
|---|
| 363 | self.order = 1 |
|---|
| 364 | |
|---|
| 365 | else: |
|---|
| 366 | self.smallsteps = 0 |
|---|
| 367 | if self.order == 1: |
|---|
| 368 | if self.order != self.default_order: |
|---|
| 369 | self.order = 2 |
|---|
| 370 | |
|---|
| 371 | |
|---|
| 372 | #Ensure that final time is not exceeded |
|---|
| 373 | if finaltime is not None and self.time + timestep > finaltime: |
|---|
| 374 | timestep = finaltime-self.time |
|---|
| 375 | |
|---|
| 376 | #Ensure that model time is aligned with yieldsteps |
|---|
| 377 | if self.yieldtime + timestep > yieldstep: |
|---|
| 378 | timestep = yieldstep-self.yieldtime |
|---|
| 379 | |
|---|
| 380 | self.timestep = timestep |
|---|
| 381 | |
|---|
| 382 | |
|---|
| 383 | |
|---|
| 384 | def compute_forcing_terms(self): |
|---|
| 385 | """If there are any forcing functions driving the system |
|---|
| 386 | they should be defined in Domain subclass |
|---|
| 387 | """ |
|---|
| 388 | pass |
|---|
| 389 | |
|---|
| 390 | |
|---|
| 391 | def update_conserved_quantities(self): |
|---|
| 392 | """Update vectors of conserved quantities using previously |
|---|
| 393 | computed fluxes specified forcing functions. |
|---|
| 394 | """ |
|---|
| 395 | |
|---|
| 396 | from Numeric import ones, sum, equal, Float |
|---|
| 397 | |
|---|
| 398 | |
|---|
| 399 | N = self.number_of_elements |
|---|
| 400 | d = len(self.conserved_quantities) |
|---|
| 401 | |
|---|
| 402 | timestep = self.timestep |
|---|
| 403 | |
|---|
| 404 | self.compute_forcing_terms() |
|---|
| 405 | |
|---|
| 406 | #Update conserved_quantities from explicit updates |
|---|
| 407 | for i, name in enumerate(self.conserved_quantities): |
|---|
| 408 | Q = self.quantities[name] |
|---|
| 409 | Q.update(timestep) |
|---|
| 410 | |
|---|
| 411 | |
|---|
| 412 | def distribute_to_vertices_and_edges(self): |
|---|
| 413 | """Extrapolate conserved quantities from centroid to |
|---|
| 414 | vertices and edge-midpoints for each volume |
|---|
| 415 | |
|---|
| 416 | Default implementation is straight first order, |
|---|
| 417 | i.e. constant values throughout each element and |
|---|
| 418 | no reference to non-conserved quantities. |
|---|
| 419 | """ |
|---|
| 420 | |
|---|
| 421 | for name in self.conserved_quantities: |
|---|
| 422 | Q = self.quantities[name] |
|---|
| 423 | Q.distribute_to_vertices_and_edges() |
|---|
| 424 | |
|---|
| 425 | |
|---|
| 426 | |
|---|
| 427 | |
|---|
| 428 | def compute_fluxes_c(domain, max_timestep): |
|---|
| 429 | """Compute all fluxes and the timestep suitable for all volumes |
|---|
| 430 | This version works directly with consecutive data structure |
|---|
| 431 | and calls C-extension. |
|---|
| 432 | Note that flux function is hardwired into C-extension. |
|---|
| 433 | """ |
|---|
| 434 | |
|---|
| 435 | import sys |
|---|
| 436 | |
|---|
| 437 | neighbours = Volume.neighbours |
|---|
| 438 | neighbour_faces = Volume.neighbour_faces |
|---|
| 439 | normals = Volume.normals |
|---|
| 440 | flux = Volume.explicit_update |
|---|
| 441 | |
|---|
| 442 | area = Volume.geometric[:,0] |
|---|
| 443 | radius = Volume.geometric[:,1] |
|---|
| 444 | edgelengths = Volume.geometric[:,2:] |
|---|
| 445 | |
|---|
| 446 | |
|---|
| 447 | flux[:] = 0.0 #Reset stored fluxes to zero |
|---|
| 448 | from domain_ext import compute_fluxes |
|---|
| 449 | timestep = compute_fluxes(domain.flux_function, |
|---|
| 450 | Volume.conserved_quantities_face0, |
|---|
| 451 | Volume.conserved_quantities_face1, |
|---|
| 452 | Volume.conserved_quantities_face2, |
|---|
| 453 | Volume.field_values_face0, |
|---|
| 454 | Volume.field_values_face1, |
|---|
| 455 | Volume.field_values_face2, |
|---|
| 456 | Boundary_value.conserved_quantities, |
|---|
| 457 | Boundary_value.field_values, |
|---|
| 458 | Vector.coordinates, |
|---|
| 459 | neighbours, |
|---|
| 460 | neighbour_faces, |
|---|
| 461 | normals, |
|---|
| 462 | flux, area, radius, edgelengths, |
|---|
| 463 | max_timestep) |
|---|
| 464 | |
|---|
| 465 | domain.timestep = timestep |
|---|
| 466 | |
|---|
| 467 | |
|---|
| 468 | |
|---|
| 469 | |
|---|
| 470 | |
|---|
| 471 | ############################################## |
|---|
| 472 | #Initialise module |
|---|
| 473 | |
|---|
| 474 | #C-extensions |
|---|
| 475 | #import util |
|---|
| 476 | #if util.can_use_C_extension('domain_ext.c'): |
|---|
| 477 | # compute_fluxes = compute_fluxes_c |
|---|
| 478 | #distribute_to_vertices_and_edges = distribute_to_vertices_and_edges_c |
|---|
| 479 | #update_conserved_quantities = update_conserved_quantities_c |
|---|
| 480 | #else: |
|---|
| 481 | # from shallow_water import compute_fluxes |
|---|
| 482 | #from python_versions import distribute_to_vertices_and_edges |
|---|
| 483 | #from python_versions import update_conserved_quantities |
|---|
