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